Timothée Lecomte

David Taj

Aristide Lemaitre

Jacqueline Bloch

Claude Delalande

Jerome Tignon

Philippe Roussignol

Polariton-polariton interaction potentials determination by pump-probe degenerate scattering in a multiple microcavity

Keywords: numbers: 68.65.-k, 73.21.-b, 78.67.Pt, 71.36.+c

We study the polarisation-dependent polariton-polariton interaction through its effect on a parametric scattering process in a microcavity (MC). The ratio of the anti-circular interaction strength V2 over its co-circular counterpart V1 is involved in defining the regime in which many nonlinear processes arise in MCs, such as parametric conversion or condensation. We measure the ratio V2/V1 using a stimulated energy-degenerate parametric scattering process in a multiple MC. The sample is pumped at normal incidence, probed with a non-zero angle, and the phase-matched idler is observed at the opposite angle. The idler behaviour, both in power and polarisation, is compared to a hamiltonian interaction model that takes into account the two polarisation-dependent parametric scattering channels characterized by V1 and V2. The proposed method to measure the ratio V2/V1 is convenient and precise. The flexibility of the triple MC allows to observe the process and measure this ratio over a large range of detunings, where we find it to be highly dependent on the detuning. These measurements complement the previous study of Vladimirova et al. [Phys. Rev. B 82, 075301 (2010)] with an original approach and for detunings that were unexplored up to now.

I. INTRODUCTION

Semiconductor microcavities (MCs) are a model system where numerous exotic optical effects can be observed: parametric amplification 1 and oscillation 2,3 , polariton condensation 4 , vortices inside condensates 5 , polariton superfluidity 6 , analogues of magnetic monopoles 7 , etc. These phenomena involve interactions between polaritons, mixed states of confined photons and quantum well (QW) excitons. Those quasi-particles inherit interesting properties from each of their two strongly-coupled components. Their excitonic part brings Coulombian interaction, and the photonic part brings cavity amplification and steep dispersion. They have two allowed spin states, and in turn the interaction effects are strongly spin-dependent. In particular, we distinguish the interaction potential V 1 between polaritons of the same circular polarisation, and the potential V 2 between polaritons of opposite circular polarisation. The behaviours of the non-linear phenomena heavily depend 8,9 on the sign of V 1 and V 2 and on the ratio V 2 /V 1 . For example, a gas of polaritons in a planar cavity at zero temperature will either condensate in the real-space (i.e. collapse spatially) or in the reciprocal space depending 9 on V 1 and V 2 . Similarly, the existence of a polarization bistability 10,11 and the feasibility of exciton polariton spin-switching 12 depend on the sign and value of V 2 /V 1 . The value of V 2 /V 1 also determines directly how the gain of a parametric amplifier depends on the polarisation configuration, as illustrated in this work.

Parametric amplification was first observed in a single MC 1 in a co-circular configuration, whereas no gain was visible in the anti-circular setup. This indicated that the potential V 2 was noticeably smaller than V 1 , as other studies confirmed 13,14 . Direct experimental studies of the relative value of V 2 followed. Based on the spin-resolved temporal evolution of the polariton populations 15,16 , it was found that V 2 was approximately ten times smaller than V 1 and of the opposite sign. The latter was confirmed by the rotation of the linear polarisation in spontaneously-formed continuous-wave optical parametric oscillation in a MC 11,17 . Studying the polarisation degree of MC luminescence 18 , the value of V 1 + V 2 was evaluated. The measure of the blueshift and of the broadening of the polariton lines in a MC with two coupled QWs under an electric field was proposed 19 as a way to measure V 1 and V 2 . Nevertheless the dependence of the values of these potentials as a function of energy has kept being neglected when interpreting MC experiments. Meanwhile, various theoretical studies focused on V 2 . In Ref. 20 the spin-dependent interaction between Wannier excitons in a crystal is estimated. In Ref. 21 a model of the spin-dependant exciton-exciton elastic scattering is proposed. In Ref. 22, the consequences of the spin-dependent interactions on the dark excitons and on the polariton condensation are studied. In Refs. 23-25 the interaction mechanisms are studied using scattering matrices made of polarisation states. The recent work of Ref. 9 focused on measuring experimentally V 2 /V 1 through the spin-resolved and energy-resolved nonlinear transmission of a MC. The ratio is found to strongly depend on the energy detuning and is compared to a theoretical model of the underlying interaction mechanisms. Since the ratio V 2 /V 1 has a major effect on the nature of processes found in MCs, we aim at using an original approach to confirm that V 2 /V 1 is not constant and to explore further the effect of the detuning.

Here, we propose to use a multiple MC as a tool to investigate the properties of the potentials V 2 and V 1 . A multiple MC is a layered structure designed to produce a rich set of resonant parametric processes between polaritons. Although parametric amplification and oscillation have been observed in single planar MCs [START_REF] Savvidis | [END_REF]26,[START_REF] Ciuti | Quattropani[END_REF] , those processes are limited to "magic-angle" configurations where the pump has to excite the sample at a very precise angle so that energy and momentum are conserved in the nonlinear process. This severely limits the range of accessible configurations. In a multiple planar MC, several cavities are stacked on top of each other and coupled by their intermediate mirrors. Appropriate coupling results in several delocalized modes of polaritons and a fan of accessible branches, enabling interbranch parametric processes to occur. The set of four-wave mixing configurations that conserve both energy and wave-vectors is consequently much richer. For example, parametric amplification and oscillation have been observed in a triple MC in an angle-degenerate configuration 3,28,29 .

In the present paper, we show that parametric scattering in a multiple MC provides a direct and precise measure of V 2 /V 1 over a large range of unexplored detunings, thus providing a reference to predict and interpret spin-dependent phenomena in MC 8 . The flexibility of the multiple MC allows us to study the potentials V 2 and V 1 over a very wide range of detunings compared to previous works. We specifically focus on a pump-probe parametric scattering that is degenerate in energy so that we can precisely describe the energy dependence of the ratio V 2 /V 1 . In this geometry, the sample is pumped at normal incidence and probed at the same energy in the lower polariton branch with a non-zero angle. The parametric scattering produces an idler at the opposite angle, phase-matched with the probe. We first verify that the process is unambiguously of parametric origin by studying the power behaviour of the ider, compared to an hamiltonian model of the parametric interaction between polaritons. Then we study the behaviour of the pumpprobe process in relation to polarisation, compared to the hamiltonian model augmented to take into account the two polarisation-dependent scattering channels described by the potentials V 1 and V 2 . We then deduce the value of the ratio of the interaction potentials V 2 /V 1 . The energydegenerate configuration makes the interpretation of our results straightforward and very slightly dependent on fit parameters. We show that the value of V 2 /V 1 strongly depends on the laser-exciton detuning.

II. PARAMETRIC DIFFUSION IN RESONANT EXCITATION

The sample used for this study is made of three coupled cavities grown using molecular beam epitaxy (Fig. 1). The three cavities are made of a total of four distributed Bragg reflectors (DBRs) with 13, 13, 13 and 25 Al 0.05 Ga 0.95 As and AlAs pairs. The structure is meant to be studied in reflection, so it was designed with a thick back mirror. The two outer cavities both contain three In 0.07 Ga 0.93 As QWs, with a 1s excitonic transition at 1.4754 eV. The inner cavity is empty. Although two cavities are enough to produce the necessary modes for a triply-resonant energy-degenerate scattering process, here a third cavity is interleaved to serve as an adjusting knob for the cavity coupling. The outer cavities are wedged in the same direction, so that they are always resonant and the position of the cavity modes can be adjusted by studying an appropriate position on the surface of the sample. The resonance of the inner cavity is placed at an energy high above the two others, and above the excitonic transition, so that the inner cavity just behaves as a mirror for the two outer cavities. It is wedged perpendicularly to the outer cavity wedges, providing an effective mirror of adjustable reflectivity, which allows to vary the optical coupling between the two outer cavities. The sample is cooled down to 5 K in a cold-finger, liquidhelium cryostat. Consequently, the structure is made of three optical modes and two excitons (one set for each outer cavity), which gives five polariton branches. The Fig. 2 shows a typical dispersion obtained by nonresonant photoluminescence, where four polariton branches are visible. This dispersion was obtained in the case where the two outer cavities and the excitons are degenerate. The fifth polariton branch is essentially composed of the inner cavity mode, and lies at a higher energy. The four visible polariton modes are mixed states of excitons and photons delocalized in the two outer cavities. The Rabi splitting, which describes the coupling between excitons and photons and corresponds to the distance between the two lower polariton and their higher counterparts, is 6 meV. The two degenerate outer cavities gives delocalized modes separated by 1 meV, as can be seen by the distance between the first two lower polariton branches.

The sample is excited in a pump-probe configuration as shown if Fig. 3. The excitation is realized with a tunable cw Ti:Sapphire laser. The laser beam is splitted in two beams both focused on the sample surface. The pump spot diameter is 50 µm. The probe spot is made appropriately smaller. Attention is paid to proper beam superposition, and correct incidence angles. The beams that are produced or reflected by the sample are observed with two CCD cameras: one is focused on the near-field and used to control the spots sizes and their superposition on the sample, and the second is focused on the far-field (i.e. the back focal plane of the excitation lens), behind a 50 cm spectrometer, to measure energy dispersions (Fig. 2) and resonant emission diagrams in the Fourier plane (Fig. 4). The pump and the probe beams are degenerate in energy. While the pump beam excites the second lower polariton branch at normal incidence, the probe beam excites the first lower polariton branch with a non-zero angle (Fig. 2). When these conditions are met, and when the pump and probe energy is chosen close enough to the exciton, an additional beam is observed coming from the sample, as shown on Fig. 4. It has the same energy as the two input beams and its wave-vector is the symmetric of the probe wave-vector with respect to the pump wavevector. This beam is identified as the idler produced by the parametric scattering of pump polaritons stimulated by the presence of probe polaritons. In the Fourier plane, the conservation of momentum imposes that the idler is always the symmetric of the probe with respect to the pump. This is verified by moving the probe beam on a ring as shown on Fig. 5. It is also possible to tilt the pump slightly, as long as the pump and the idler remain close to the modes of the structure. We then observe that the idler follows the tilt of the pump by staying the symmetric of the probe.

III. MODEL OF INTERACTING POLARITONS

We will compare the experimental results to a model describing the polaritonic excitations in second quantization with the hamiltonian H = H 0 + H int obtained by truncating the coupled excitons-photon hamiltonian to the second order in exciton density, and developing it in the basis of polaritons 21,30 . H 0 gives the linear behaviour of the polaritons and H int describes the interaction between polaritons, that comes from their excitonic part:

H 0 = σ=± ω p † pσ p pσ + p † sσ p sσ + p † iσ p iσ (1)
H int = σ=± V 1 (q = 0)p † pσ p † pσ p pσ p pσ + V 2 (q = 0)p † pσ p † p-σ p pσ p p-σ + V 1 (q = κ)(p † pσ p † pσ p sσ p iσ + h.c.) + V 2 (q = κ)(p † pσ p † p-σ p sσ p i-σ + h.c.) (2)
p p/s/i,σ are the annihilation operators for the pump, probe (or signal) and idler polaritons respectively, with circular polarisation σ = ±. In our degenerate configuration, ω = ω p = ω s = ω i is the single particle energy of the modes. V 1/2 are the interaction potentials for the co-circular and anti-circular polarisation scattering channels, respectively. The first two terms in H int correspond to the Kerr effect (or blueshift ) caused by the interaction of pump polaritons between each other. The last two terms correspond to the parametric scattering of two pump polaritons into a signal polariton and an idler polariton, respectively conserving or inverting the circular polarisation. q = κ is the distance in momentum between the pump and the signal. H int contains only the four interaction channels that are of importance in our configuration. We have neglected the terms for the cross Kerr effect ( p † pσ p † pσ p sσ p sσ and p † pσ p † pσ p iσ p iσ ) and all other terms where the signal and the idler operators appear at a superior order. This hamiltonian gives a set of coupled differential equations for the mean fields which boils down to the linear system M v = S, for the time-independent complexvalued vector v = (p s+ , p * i+ , p s-, p * i-), where M is:

M =     Ω V 1 (p p+ ) 2 0 V 2 p p+ p p- V 1 (p * p+ ) 2 Ω * V 2 p * p+ p * p- 0 0 V 2 p p+ p p- Ω V 1 (p p-) 2 V 2 p * p+ p * p- 0 V 1 (p * p-) 2 Ω *    
(3) Here we use Ω = (ω 0 +δω 0 -ω)-iγ, where ω 0 is the bare frequency of the signal and idler modes, δω 0 takes into account possible blueshifts of the resonance condition, γ is the signal and idler inverse lifetime and ω is the laser frequency. So (ω 0 + δω 0 -ω) corresponds to the detuning between the laser and the polaritonic modes.

We study the behaviour of the total idler intensity I i = |p i+ | 2 +|p i-| 2 by inverting M , using a source vector S describing a vertically polarised probe : S = (1, 0, 1, 0) and a pump linearly polarised with an angle θ p : p p+ = e iθp , p p-= e -iθp . We assume that the resonance condition is met: Ω = -iγ. Solving this system for pump powers far below the oscillation threshold, we obtain an expression of the idler intensity I i that shows how it depends on the FIG. 6. Idler power Ii plotted against the probe power Is. The dots are experimental data points, the solid curve is a linear fit to the low power data. See the text for the value of the fixed parameters. pump intensity I p , signal (i.e. probe) intensity I s , on the pump-probe angle θ p , on the interaction potential V 1 and on the ratio V 2 /V 1 :

I i ∝ I s I 2 p V 2 1 1 + V 2 V 1 2 + 2 V 2 V 1 cos(2θ p ) (4)

IV. POWER BEHAVIOUR

We now study the experimental dependence of the idler intensity with respect to the excitation parameters. First, the idler intensity dependence on the probe power is shown in Fig. 6. The other parameters are fixed: the detuning δ = E laser -E exciton is -2.5 meV, the pump power is 50 mW and the pump and the probe beams are cross-polarised. The Eq. 4 anticipates a linear dependence for that measurement: I i ∝ I s . This is what we observe up to a probe power of 8 mW. At higher probe powers, the idler intensity saturates. This indicates that the process efficiency is limited: an additional stimulation above 10 mW does not increase the scattering rate. Competing processes can explain this saturation, including mode shifts, pump depletion and crossed two-photon absorption of the probe with pump polaritons. [START_REF] Shen | The Principles of Nonlinear Optics[END_REF] Second, the idler intensity dependence on the pump power is shown in Fig. 7. In this set of data, the probe power is 2 mW, the detuning δ = E laser -E exciton is -2.8 meV, and the pump and the probe beams are crosspolarised. The Eq. 4 anticipates a quadratic dependance for that measurement:

I i ∝ I 2
p . The data is correctly fitted with a square power law up to a pump power of 70 mW. For a pump power above 80 mW, the idler intensity decreases dramatically, showing that the parametric process stops. This is due to non-linear effects induced by the pump on the modes of the structure: the polariton branches widen and shift, the weak-coupling regime The dots are experimental data points, the solid curve is a square-law fit to the low-power data.

is approached, and consequently the probe beam is no longer in resonance and cannot enter the structure.

Third, the idler intensity dependence on the detuning δ = E laser -E exciton is shown in Fig. 8. Here, the pump power is fixed at 60 mW. The probe power is 2 mW. The probe and pump beams are cross-polarised. The enveloppe of this curve is bell-shaped with a maximum around zero detuning. For detunings lower than -3 meV, the idler was not visible. Variations below the enveloppe of Fig. 8 are attributed to situations where the spot superposition or the angles of incidence were sub-optimal. We observed that the process works at detunings which are close to zero in a weak-coupling configuration, where the beams can enter the sample through the cavity modes only weakly-coupled to the excitons. Still, the excitonic nonlinearity allows the parametric process to occur efficiently.

In the Eq. 4, the terms that directly depend 32 on the detuning δ are V 1 and V 2 . The experimental setup corresponds to θ p = π/2, which gives

I i ∝ (V 1 -V 2 ) 2 .
Now, each of these potentials depends on the detuning δ through the Hopfield coefficient [START_REF] Hopfield | [END_REF] for the excitonic part of the polaritons: V i ∝ X(δ) 4 . This gives:

I i ∝ X(δ) 8 = Ω ( Ω) 2 + 4(E -E X ) 2 8 ( 5
)
On Fig. 8, the fit with X(δ) 8 correctly reproduces the envelope of the data points.

All these observations validate our theoretical description of the scattering process.

V. POLARISATION BEHAVIOUR AND V2/V1 MEASUREMENTS

Now that it is clear that the process that we observe is parametric and that the Eq. 4 satisfyingly describes the power dependences, we focus on the polarisation behaviour that will lead to measurements of the ratio

V 2 /V 1 .
The fact that two polarisation-dependent scattering channels are present induces a rotation effect on the polarisation of the idler. At least at the lowest order in I p , the interaction model anticipates that the idler is the symmetric of the signal with respect to the pump. This translates on the linear polarisation angles as:

θ i -θ p = θ p -θ s (6) ⇔ θ i = 2θ p -θ s (7)
The Fig. 9 presents a serie of measurements of the idler polarisation angle when the pump polarisation is rotated. The pump power is 60 mW, the probe power is 4 mW, and the laser-exciton detuning is -2.6 meV. The data shows that the idler angle θ i varies with 2θ p -θ s .

Now, looking at Eq. 4, it is clear that measuring how I i depends on θ p should provide enough data to estimate

V 2 /V 1 , since I i ∝ 1 + (V 2 /V 1 ) 2 + 2(V 2 /V 1 )
cos(2θ p ). The Fig. 10 shows one serie of raw measurements of I i as a function of θ p . The pump power is 60 mW, the probe power is 4 mW and the laser-exciton detuning is -2.6 meV. Although we see an obvious dependence on θ p , the period of this dependence is twice as small as the one of cos(2θ p ). We see two local maxima for parallel and perpendicular polarisations and a minimum around 40 degrees. This pronounced minimum at 40 degrees comes from the polarisation response of our visualisation setup, where the beams reach the CCD sensor used for intensity measurements by first going through various waveplates and being reflected on the spectrometer grating. This grating and the waveplates are sensible to the angle of linear polarisation of the idler. The waveplates also introduce a small elliptical component to the idler polarisation. The effective transfer function for these effects can be described with three separately-calibrated parameters: first, the ratio of the setup response with a horizontal polarisation with respect to a vertical polarisation ρ = horizontal /vertical = 0.068; second, a constant angle shift θ 0 = -3.2 degrees which corresponds to the difference between the vertical as measured by the reference polariser and as defined by the slit and the grating of the spectrometer; third, an elliptical component represented by an angle ψ = 15 degrees. The polarisation response of the setup brings an additional multiplicative term f (θ i ) in the expression of I i :

f (θ i ) ∝ cos 2 (ψ sin(2(θ i -θ 0 ))) cos 2 (2(θ i -θ 0 )) + sin 2 (ψ sin(2(θ i -θ 0 ))) sin 2 (2(θ i -θ 0 )) + ρ cos 2 (ψ sin(2(θ i -θ 0 ))) cos 2 2(θ i -θ 0 ) - π 2 + ρ sin 2 (ψ sin(2(θ i -θ 0 ))) sin 2 2(θ i -θ 0 ) - π 2 (8)
The shape of f (θ i ) correctly describes the minimum at 40 degrees. Now the assymetry between 0 and 90 degrees is not explained at all by the polarisation response of the setup but by the Eq. 4 when V 2 is non-zero. We finally fit the experimental points in Fig. 10 with:

af (θ i ) 1 + V 2 V 1 2 + 2 V 2 V 1 cos(2θ p ) ( 9
)
where a is a multiplicative constant. The fit leads us to a relative measure of V 2 for this particular detuning:

V 2 ≈ -0.19V 1 (10)
Alternatively we show on Fig. 11 the experimental measurements normalised by the setup response f (θ i ) to better highlight the variation that solely comes from the parametric interaction model. The fit is directly done with The dots are the experimental data. The probe polarisation is fixed vertically, while the pump polarisation is rotated incrementally between 0 degree (parallel) and 90 degrees (orthogonal).The solid curve is a fit including the effect of the two polarisation-dependent scattering channels.

I i ∝ 1 + (V 2 /V 1 ) 2 + 2(V 2 /V 1 ) cos(2θ p ), so that the unique adjustable parameter is V 2 /V 1 .

VI. V2/V1 AS A FUNCTION OF THE DETUNING

Since the value of V 2 /V 1 has been shown to depend on the detuning in previous studies 9 , we now use our setup to measure that value over a range of detunings that have not been explored up to now. We conducted this study over a range of laser-exciton detunings going from -2 meV to +5 meV. The results are reproduced in Fig. 12. We 12. Ratio V2/V1 of the polarisation-dependent scattering strengths measured over a range of laser-exciton detunings δ = Ep -EX. The square data point corresponds to the measurements presented in Fig. 10. The solid line is a linear regression. The shaded area is a 95% confidence region for this linear regression, bounded by the 2.5% and 97.5% quantiles.

find that V 2 /V 1 increases progressively from -0.2 to 0 on this range of detunings.

This range of values is in good agreement with the previous estimation in Ref. 34, which gives V 2 ≈ -0.04V 1 in a single MC with one QW for a fixed cavity-exciton detuning of 0 meV, which approximately corresponds to a laser-exciton detuning δ of -2 meV in our sample. Ref. 16 gives V 2 ≈ -0.08V 1 for a fixed cavity-exciton detuning of 0 meV, which is in even better agreement with our measurements.

In the range of detunings studied here, our measurements show that the ratio V 2 /V 1 is increasing with the detuning. Assuming that the dependency of the ratio V 2 /V 1 as a function of detuning is linear in this range, we have run a regression analysis on the measurements, finding an average increase from -0.175 to -0.01. Performing a Monte-Carlo Bayesian inference regression [START_REF] Gelman | Bayesian Data Analysis, Second Edition, Chapman & Hall/CRC Texts in Statistical Science[END_REF][START_REF] Development | Stan: A C++ Library for Probability and Sampling[END_REF] , we compute the probability distribution for the parameters of the fit to quantify how confident we can be in the increase of the ratio V 2 /V 1 . The shaded area in Fig. 12 represents the 95% confidence region for the regression, bounded by the 2.5 % and 97.5 % quantiles obtained by the Bayesian method. This statistically confirms that the interaction ratio is actually increasing with the detuning.

The measurements in this paper shall now be compared to those obtained in Ref. 9. This work indicates that V 2 /V 1 decreases progressively from 0 to -1 for cavityexciton detunings going from -3.0 meV to 0.0 meV. Given the Rabi splitting of 3.5 meV in their sample, this range of cavity-exciton detunings corresponds to laser-exciton detunings δ going from -4.0 to -1.2 meV. For larger detunings, the measurements of Ref. 9 are highly dispersed and no tendency stands out. In contrast, our method provides data for detunings δ above -1.5 meV and up to 4.5 meV and indicates that V 2 /V 1 increases progressively from -0.2 to 0.0. Our study confirms the major result of Ref. 9: the ratio V 2 /V 1 is not a constant. We note three differences in the results of the two studies : we have studied a different range of detunings, we see a lower amplitude of variation of the ratio V 2 /V 1 (0.2 meV versus 1 meV), and we see this ratio increasing with detuning whereas it was found to be decreasing in Ref. 9. It is worth noting that these two studies are done in very different conditions. Not only the accessible detuning range is different, but there are also differences in the sample structure: we study a multiple MC instead of a single MC and our structure contains six QWs instead of one. These choices provide us more flexibility and more efficiency in the nonlinear processes. It is also worth noting that we study the sample using a resonant continuous incident light source instead of a pulsed source with a wide spectrum. The single-vs-multiple QWs situation is the main candidate for the difference in behaviour. We conclude this comparison by emphasizing that the method we used is experimentally convenient and provides straightforward measurements of V 2 /V 1 over a wide range of parameters.

The values of V 2 /V 1 measured in this work have interesting consequences. We mention here what these consequences are in the context of parametric oscillation. In the parametric scattering process which can lead to parametric oscillation, the gain is proportional to V 1 + V 2 when the linear polarisation is conserved, and to V 1 -V 2 when the linear polarisation is reversed. Since we have measured that V 2 /V 1 is negative and its absolute value is smaller than 1 in the range of detunings studied, it is in the configuration that reverses the linear polarisation that the parametric oscillation is expected to happen. This corresponds well to experimental observations 3,[START_REF] Dasbach | [END_REF]38 . This polarisation selectivity is all the more pronounced when the pump energy is well below the excitonic transition. On the other hand, if we extrapolate the data points for laser-exciton detuning above 5 meV, we can expect V 2 to become positive. It means that for these large detunings the regime of parametric oscillation would be expected to change dramatically and the process would conserve the linear polarisation.

CONCLUSION

We have reported the study of the polarisationdependent polariton-polariton interaction potentials in a triple-MC by observing the energy-degenerate stimulated parametric scattering of polaritons. In this process, a probe beam stimulates the scattering of two pump polaritons into a signal and an idler polaritons of the same energy, conserving energy and momentum. The idler behaviour follows an interaction model that takes into account two polarisation-dependent scattering channels. This has allowed us to measure the relative value of the anti-circular interaction potential V 2 compared to its co-circular counterpart V 1 with a convenient and pre-cise method. We analysed this ratio over a large range of detuning and found it to be highly detuning-dependant: it goes from -0.2 to 0 for laser-exciton detunings going from -2 meV to +5 meV. This study complements earlier results 9 in a range of detunings unexplored to-date. This shows the flexibility of a multiple MC, and we believe these results are precious for the study of spin-dependent nonlinear phenomena in MCs.

4

  FIG. 1. (Color online) Structure of the asymmetric triple MC.

FIG. 2 .

 2 FIG. 2. Dispersion of the polaritonic modes obtained by photoluminescence under nonresonant excitation. Four polariton branches are visible. The energy-degenerate scattering process at stakes is schematically represented on the dispersion.

FIG. 3 .

 3 FIG. 3. (Color online) Excitation and visualisation setup for the reflective pump-probe experiment.

  FIG. 4. (a) Image of the parametric scattering process in the Fourier plane of the excitation lens. The polarisations of the pump and the probe are perpendicular. The idler has the same polarisation as the probe. A polariser blocks most of the reflected pump. (b) Schematic view of the energy dispersion of the polaritonic modes against energy and angle and of the parametric process.

FIG. 5 .

 5 FIG. 5. Images of the parametric scattering process in the Fourier plane, in conditions similar to those of Fig. 4. The probe is moved over a ring centred on the pump, while the spots of the probe and the pump remain superposed on the surface of the sample. The idler follows the probe position in the Fourier plane, showing the momentum conservation.

20 FIG. 7 .

 207 FIG. 7. Idler intensityIi plotted against the pump power Ip. The dots are experimental data points, the solid curve is a square-law fit to the low-power data.

FIG. 8 .

 8 FIG. 8. Idler intensity Ii normalized by the probe intensityIs and plotted against the laser detuning δ = Ep -EX. The solid curve is a fit with X(δ)8 .

FIG. 9 .

 9 FIG.9. Linear polarisation angle θi of the idler as a function of the linear polarisation angle of the pump θp. The linear polarisation of the probe is kept vertical while the pump angle is rotated from 0 degree (parallel to the probe) to 90 degrees (perpendicular). The dots are the experimental data points. The solid line is a fit with θi = 2θp -θs.

FIG. 10 .

 10 FIG.10. Idler intensity Ii plotted against the relative polarisation θp of the pump and the probe. The dots are the raw experimental data. The probe polarisation is fixed vertically, while the pump polarisation is rotated incrementally between 0 degree (parallel) and 90 degrees (orthogonal). The solid curve is a fit including the polarisation response of the visualization setup and the existence of the two polarisation channels. The dashed curve illustrates the same fit for a flat polarisation response. The dotted curve illustrates the same fit when neglecting V2, i.e. when ignoring the scattering channel that inverses circular polarisation.

  FIG.12. Ratio V2/V1 of the polarisation-dependent scattering strengths measured over a range of laser-exciton detunings δ = Ep -EX. The square data point corresponds to the measurements presented in Fig.10. The solid line is a linear regression. The shaded area is a 95% confidence region for this linear regression, bounded by the 2.5% and 97.5% quantiles.


Pierre-Jean Meyer

Hosein Nazarpour

Antoine Girard

Emmanuel Witrant

Experimental Implementation of UFAD Regulation based on Robust Controlled Invariance*

In this paper, we implement a temperature regulation strategy on a small-scale building equipped with Under-Floor Air Distribution (UFAD). Using a 0-dimensional model of the temperature variations in the building identified from experimental measurements, this strategy is based on the criteria of Robust Controlled Invariance, which give a characterization of the limits for a robust controller to keep the state of the system in an interval. These criteria are defined for monotone systems with bounded disturbances and thus are independent of the chosen control strategy.

I. INTRODUCTION

Energy consumption in buildings (both commercial and residential) has reached up to 40% of the total energy use in developed countries and is rapidly increasing with the growing demand in comfort [START_REF] Perez-Lombard | A review on buildings energy consumption information[END_REF]. This fact motivates many countries to make energy savings in buildings a priority, as mentioned by the European Energy Performance of Buildings Directives [2], [3]. As a result, a substantial amount of work has been done in the past decades toward intelligent buildings [START_REF] Wong | Intelligent building research: a review[END_REF]. This has been particularly the case for climate regulation in buildings, with research on modeling, simulation [START_REF] Trčka | Overview of HVAC system simulation[END_REF] and control [START_REF] Mirinejad | Control techniques in heating, ventilating and air conditioning systems[END_REF] of Heating, Ventilating and Air Conditioning (HVAC) systems, to improve comfort and energy efficiency. On these matters, the UnderFloor Air Distribution (UFAD) solution has shown some interesting results compared to the more traditional ceiling ventilation [START_REF] Bauman | Underfloor air distribution (UFAD) design guide[END_REF].

Various paths have already been explored for the control of HVAC systems in intelligent buildings. When the focus is mainly on control, numerous feedback strategies have been devised, based on simple PID or On/Off control, more robust controllers with the H ∞ approach [START_REF] Witrant | Limitations and performances of robust control over WSN: UFAD control in intelligent buildings[END_REF], or non-linear approaches [START_REF] Arguello-Serrano | Nonlinear control of a heating, ventilating, and air conditioning system with thermal load estimation[END_REF]. For more energy-efficient controllers, we can look for the optimal tradeoff between comfort and energy saving [START_REF] House | Optimal control of building and HVAC systems[END_REF], a model-predictive strategy [START_REF] Oldewurtel | Energy efficient building climate control using stochastic model predictive control and weather predictions[END_REF], or a fuzzy logic controller [START_REF] Hamdi | A fuzzy control system based on the human sensation of thermal comfort[END_REF].

In this paper we consider the notion of Robust Controlled Invariance which characterizes whether a controller can maintain the state of the system in an interval for any value of the disturbances. With the monotonicity property [START_REF] Angeli | Monotone control systems[END_REF], intuitive for thermal systems, we obtain simple new criteria on Robust Controlled Invariance. The goal of this paper is to use these criteria to implement a temperature regulation *This work was partly supported by a PhD scholarship and the research project COHYBA funded by Région Rhône-Alpes.

† Univ. Grenoble Alpes, LJK, F-38000 Grenoble, France:

{Pierre-Jean.Meyer, Antoine.Girard}@imag.fr nazarpur@gmail.com ‡ Univ. Grenoble Alpes, GIPSA-Lab, F-38000 Grenoble, France:

Emmanuel.Witrant@gipsa-lab.grenoble-inp.fr on a small-scale UFAD flat. We should note that our focus is not on the efficiency of a controller but simply on checking if the invariance can be achieved in the chosen interval. For a deeper discussion on the monotonicity of our system and the results involving the Robust Controlled Invariance, the reader is referred to the theoretical counterpart of this paper in [START_REF] Meyer | Controllability and invariance of monotone systems for robust ventilation automation in buildings[END_REF] 1 . The paper is organized as follows. In Section II, we describe the 0-dimensional UFAD model from [START_REF] Witrant | Limitations and performances of robust control over WSN: UFAD control in intelligent buildings[END_REF], which is then identified with a small-scale experiment of a flat to match its behavior. In Section III, the monotonicity property is presented and analyzed for our 0D model, followed by the introduction of the Robust Controlled Invariance. The final Section IV gives an application of the obtained control criteria on the experiment and a comparison between these results and the identified model.

II. MODELING

The system considered is a small-scale flat equipped with UnderFloor Air Distribution (UFAD), as sketched in Fig. 1 and pictured in Fig. 2. The cooling of the room temperatures is done as follows. Peltier coolers are used to decrease the temperature of the air in the underfloor plenum, which is sent into each room of the building using the controlled fans. The excess of air is pushed into the ceiling plenum through exhausts in the fake ceiling and then sent back to the underfloor through a return pipe outside of the flat.

There are two components in the control problem. At the building level, we control the air recirculation with the speed of the fan in the return pipe and the temperature in the underfloor plenum by acting on the Peltier coolers. At the room level, we use the fans in the underfloor separately to Fig. 2. Small-scale experiment of a flat with UFAD regulation2 obtain a decentralized control of each room temperature. In this paper, the controls at the building level are assumed to be set and our focus is on achieving separate climate regulation using decentralized control.

A. Theoretical model

We consider the following 0-dimensional model:

∀i, dT i dt = j∈Ni a i,j (T j -T i ) + b i u i (T u -T i ) (1)
+ j∈Ni δ dij c i,j * h(T j -T i ) + δ si d i (T 4 si -T 4 i ),
with all constants a, b, c, d being positive. A detailed description on establishing this model is given in [START_REF] Witrant | Limitations and performances of robust control over WSN: UFAD control in intelligent buildings[END_REF] and general information on the heat transfers involved in (1) can be found in the following books on the design of HVAC systems [START_REF] Levermore | Building Energy Management Systems: Applications to Low-Energy HVAC and Natural Ventilation Control[END_REF], [START_REF] Mcquiston | Heating, ventilating, and air conditioning: analysis and design[END_REF]. This model is obtained by combining energy conservation (the first law of thermodynamics) and mass conservation in each room i of the building. It is assumed that: the temperature of air is uniform in each room; the air is incompressible; the potential and kinetic energies of air can be neglected due to its reduced mass and speed.

Here, we review some of the hypotheses used to derive (1) in [START_REF] Witrant | Limitations and performances of robust control over WSN: UFAD control in intelligent buildings[END_REF] in order to capture the heat transfers observed in the experiment.

The first heat transfer in [START_REF] Perez-Lombard | A review on buildings energy consumption information[END_REF] corresponds to the conduction between the considered room i and the neighbor rooms (or the outside) j ∈ N i , through the PVC walls. The second term is linked to the mass flow rate u i forced by the fan from the underfloor (of temperature T u ) into room i. The command over the fan speed, and by extension the induced air flows u i , are the control inputs of our system. On the second line of (1) are the disturbances: the flow going through open doors and the radiation from heat sources (human body, computer, . . . ). The booleans δ dij and δ si describe the discrete state of these disturbances: they are equal to 1 when the door between rooms i and j is open and when a heat source is active in room i, respectively. The function h introduced in the term for the doors is defined by:

h(x) = 0 if x ≤ 0 x 3/2 if x > 0,
and it comes from the removal, using the mass conservation, of the flow going into the ceiling. As a result, the door heat transfer only appears in the equation of the coldest of the rooms connected by the door. Finally, our experiment has lamps in each room included as the radiating heat sources of temperature T si .

In this paper, we write the system dynamics as

dT dt = f (T, u, w, δ), (2)
where the four-dimensional vector field f depends on the state T = [T 1 , T 2 , T 3 , T 4 ], the control input u = [u 1 , u 2 , u 3 , u 4 ], the boolean disturbances (doors and sources) δ, and the exogenous input vector w corresponding to the other temperatures (T u for the underfloor, T c for the ceiling, T o for the outside). Each component f i of f follows (1).

B. Process identification

The theoretical model ( 1) is adapted to the measured behavior of our experiment with the following identification strategy. We run several experiments on the flat to capture the main behaviors modeled in [START_REF] Perez-Lombard | A review on buildings energy consumption information[END_REF]. In these tests, the outside temperature T o is varying around 30 • C and the underfloor temperature T u is regulated at 17 • C using a PID controller. Our experiments are quickly summarized as follows, with the desired physical phenomena.

• To observe the heat radiation, in each room, we turn on a lamp, wait for an equilibrium to be reached, then turn the lamp off; • for the air flow coming from the underfloor, in each room, we turn on a fan, wait for an equilibrium to be reached, then turn the fan off; • for the air flow through a door, we heat one of the neighbor room with a lamp to create a temperature gradient, then open the door; • similarly to the previous test, we create a temperature gradient at the door by cooling down a neighbor room using a fan, then opening the door. This test is also added to check the influence of the fans on the direction of the air flow through the door; • in each room, we alternatively turn on and off both fans and lamps to generate a data set that is sufficiently representative of the different operating conditions.

All these tests aim to quantify heat transfers due to heat radiation and the exchange of air flows, as well as the conduction in the walls. The control of the experiment (fans, lamps, doors) and the measurements are done using LabVIEW TM . The data are then processed in MATLAB . In our problem we have 40 unknown parameters, with 10 for each room: 5 for the conduction (underfloor, ceiling, outside, and 2 neighbor rooms), 1 for the fan, 1 for each door, and 2 for the heat source (emissivity and temperature). We solve this optimization problem with a least-squares algorithm initialized with a set of values based on known physical parameters and observations. We run this optimization on 45652 algebraic relationships obtained from the discretized system, representing the 57079 data points from the above experiments. The resulting values of the identified model ( 1) are given in Table I in the appendix.

C. Model evaluation

The identified model is evaluated on experimental scenarios not included in the identification data set. An example of such scenario is described below, where we start with all lamps and fans off and all doors closed: The data measured during this experiment (noisy, blue) is displayed and compared to the model (dashed, red) in Fig. 3. The vertical lines correspond to the transitions described in the previous paragraph: solid when the switched element has a direct link to the corresponding room and dashed otherwise. For the most part, we can see in Fig. 3 that the model fits reasonably well the data. For this data set of 1211 points, the mean squared error between the data and the model is 0.18 • C, with a standard deviation of 0.42.

t = 150s, lamp 1
Even though this data set is representative of the typical phenomena observed in the building, it is only a sample of all the possible combinations between doors, fans and lamps to be verified. In a more thorough verification, we found several limitations to our model, some of which are visible in Fig. 3. First, if we look at the behavior in room 1 (top left of Fig. 3), the temperature of the model keeps increasing without any influence (i.e. cold air inflow) from the open doors with rooms 2 or 4, since room 1 has the highest temperature (0D hypothesis that the air flows from the hot volume toward the cold one, due to pressure difference at constant density). In reality, we can see on the measured data that the door opening has some circulation effect on the hottest room. This is particularly visible between 810 and 1050 seconds, when the door between rooms 1 and 4 is open and the fan is running in room 4. This fan generates an air flow going in the opposite direction of the flow created by the pressure difference at the door. This behavior is not modeled in (1), and may become a problem when using the fans too fast. An extra term representing the turbulent diffusive mixing could be introduced if a refined model is required.

Another difference between the model and the data can be seen on the graph for room 3 (bottom left of Fig. 3). In this room, while the fan is running, the lamp is turned on at 570 seconds and turned off at 930 seconds. The difference between the model and the data in that case can be explained by the fact that the fan forces convection on the lamp. Indeed, the moving air creates a new convective heat transfer between the hot lamp (even when it has been turned off) and the air.

Also, the identified model is based on the ON/OFF use of the fans. The air flow induced by the maximal fan voltage is thus part of the parameters to identify and is included in the obtained values for b i in (1). In the next sections, when we apply a voltage between the extremal admissible values, we consider in the model that the air flow is linear to the voltage. To model the friction in the fans, we add a lower saturation of the air flow (u i = 0) when the voltage is below a threshold.

III. MONOTONICITY AND ROBUST CONTROLLED INVARIANCE

For rigorous definitions and proofs of the notions presented in this section, see the extended version of [START_REF] Meyer | Controllability and invariance of monotone systems for robust ventilation automation in buildings[END_REF] available at the link given in page 1.

A. Monotonicity

To obtain the controlled invariance proposed in the next section, the system needs to satisfy the monotonicity property. In [START_REF] Angeli | Monotone control systems[END_REF], Angeli and Sontag define a monotone system as a system preserving partial orderings chosen for the state and the inputs. The cooperative systems, a subclass of the monotone systems but with a very similar definition, take their names from an interpretation which gives a clearer understanding of how such systems behave. A system is called cooperative because its variables reinforce each other, or more generally, "each pair of variables may affect each other in either positive or negative forms" [START_REF] Angeli | Interconnections of monotone systems with steady-state characteristics[END_REF].

This property is particularly easy to observe on thermal systems: if we take the example of the conduction a i,j (T j -T i ) (when considering room i), we naturally see that the hotter the neighbor T j , the bigger the heat transfer toward room i. This comment can be reproduced for any of the heat transfers modeled in (1), whether T j represents a state or a temperature from w. In a similar way, for the discrete disturbances δ: if more δ dij and δ si are equal to 1, we have more heat transfers going toward the considered room i, reinforcing the state T i . This is natural for heat sources, since it is assumed that the source temperature is always hotter than the room temperature. In the case of the doors, remember that the heat transfer in (1) only appears in the coldest room. So having δ dij = 1 with T j > T i indeed reinforces T i , while having no effect on T j .

The variables T , w and δ thus all affect the state T in a positive way. In the case of our flat, we aim to cool down the four rooms, so the underfloor temperature is necessarily set to a value below any room temperature: T u < min i (T i ). Then, the heat transfer in room i corresponding to the air flow forced by the fan is always negative (b i u i (T u -T i ) < 0), resulting in a decrease of T i . The control input u thus affects the state T in a negative way: increasing the air flow u i would decrease the room temperature T i faster.

The fact that the variables affect the state in different ways is not a problem: what matters is that each variable always has the same effect, which is indeed our case. According to the interpretation given in the first paragraph, our system is cooperative, and therefore monotone.

We can now address the initial definition of a monotone system, involving preservation of the partial orderings chosen for the state and the inputs. Let's consider two sets of variables (T, u, w, δ) and (T , u , w , δ ) where T is an initial state and u, w, δ are functions of time. We can choose partial orderings consistent with the above observations:

         T > T ∀t > 0, u(t) < u(t) ∀t > 0, w(t) > w(t) ∀t > 0, δ(t) > δ(t) .
In that case, knowing that the system is monotone means that the trajectory of the system with (T, u, w, δ) stays at all time above the trajectory of the same system with (T , u , w , δ ).

B. Robust Controlled Invariance

The Robust Controlled Invariance is the main notion used in the remaining of this paper. We aim to use the control input u to keep the state in a chosen target interval [T , T ] for any external conditions (w, δ).

We consider that all the inputs of our system are bounded: either by physical constraints for the control inputs (u i ∈ [0, u i ]) and the disturbances (δ ∈ [0, 1]); or by a range obtained from observations: w ∈ [w, w]. With that assumption, we can use the monotonicity to simplify the conditions for Robust Controlled Invariance by considering the worst case for the external conditions (w and δ) while taking the best control input u to counteract the situation induced by the disturbances. If we take the example of the upper bound T of the interval, the worst situation of the disturbances is the hottest one (w = w and δ = 1) which tries to push the state above the interval. To compensate this effect, we choose the best ventilation u = u to cool down as much as possible. If in that situation we are able to decrease the temperature, it means that even in the hottest conditions, we can prevent the state from going over the upper bound T .

We can now characterize the Robust Controlled Invariance using componentwise inequalities on the vector field.

Theorem 1: The system is Robust Controlled Invariant in [T , T ] if and only if f (T , u, w, 1) ≤ 0 f (T , 0, w, 0) ≥ 0.

To simplify the notations in what follows, we define the subspaces C and W as the sets of admissible upper bound T and lower bound T respectively, so that the conditions in Theorem 1 can simply be written T ∈ C and T ∈ W.

IV. APPLICATION TO THE UFAD EXPERIMENT

A. Controller implementation

When the state space is of dimension 2, we can easily obtain a graphical representation of the subspaces C and W, which facilitates the choice of the target interval [START_REF] Meyer | Controllability and invariance of monotone systems for robust ventilation automation in buildings[END_REF]. In our case, such a graphical representation is more complicated to obtain due to the 4 dimensions. Instead, we simply compute the extremal points of these subspaces (lower bound of C and upper bound of W). These values would represent the minimal allowed T (lower bound of C), and the maximal allowed T (upper bound of W) in order to obtain the Robust Controlled Invariance.

The characterization of the subspaces C and W (and therefore of the Robust Controlled Invariance) requires to consider that all the inputs of the system are bounded. The control input u is clearly bounded between 0 and the air flow induced by the maximal voltage applied to the fans. The discrete disturbances in δ can only take two values: 0 or 1. For the exogenous temperatures w, we take the intervals corresponding to observations in similar situations: both the ceiling and the outside temperatures are assumed to be varying between 27 • C and 30 • C. For the underfloor temperature, we use a PID controller to maintain it at 17 • C. However, it has been observed that when all the fans are active, the Peltier coolers are not powerful enough to cool down the hot air coming from the ceiling. That is why we consider the following boundaries for the underfloor temperature: T u ∈ [START_REF] Angeli | Interconnections of monotone systems with steady-state characteristics[END_REF]21].

Knowing these boundaries, we can now obtain the lower bound of C: min(T ) = [26.1, 26.2, 29.0, 28.1]; and the upper bound of W: max(T ) = [26. 2, 26.3, 26.4, 26.5]. The smallest Robust Controlled Invariant interval is wider for rooms 3 and 4 because room 3 is the smallest room and the heat radiated by the lamp in room 4 is more important than in other rooms. When we consider the hottest situation (C), the heat radiation from the lamps thus has more impact on rooms 3 and 4. For this control experiment, we choose a 3 • C-wide target interval for each room, apart from room 3 for which we take it 4 • C-wide:

T =     26 25 26 26     T =     29 28 30 29     . (3)
In Theorem 1, only one condition is implied on the controller: to use the maximal ventilation u i when a room temperature reaches its upper bound T i and to stop the fan when T i = T i . This leaves the choice of the control strategy relatively free to meet some performance specifications. In the scope of this paper, we do not aim to discuss the efficiency of a controller but simply to check the Robust Controlled Invariance with any controller. This is why we choose a simple linear controller, with saturations to satisfy the condition above. In the model described in (1), we consider that the control input is the mass flow rate u sent by the fans. In reality, we do not have a direct control on the air flow, but only on the voltage V applied to the fans. We consequently control the fan voltages in this Decentralized Linear Saturated (DLS) strategy.

Definition 1 (Decentralized Linear Saturated Controller):

∀i,

       T i ≤ T i ⇒ V i = 0 T i ∈ [T i , T i ] ⇒ V i (T i ) = V i * Ti-Ti Ti-Ti T i ≥ T i ⇒ V i = V i
As explained at the end of Section II-C, we only used the control input in ON/OFF mode until now. Indeed, since intuitively the air flow is monotone with respect to the voltage applied to the fan, our only concerns to ensure the invariance are the extremal values. However the controller from Definition 1 can assign to the fan any value in the voltage range [0, 6], so if we want to compare the measurements from the controlled experiment with the corresponding simulation we need to add to the model a function giving the relation between voltage and air flow of the fan. First, to model the frictions preventing the fan from working with a voltage below a threshold, we add a lower saturation (u i = 0) for these small voltages. The observed voltages for the threshold in each room are: [2, 3.5, 2.5, 2.5]. To keep the model as simple as possible for the remaining voltage range, we consider that the air flow is linear in the voltage, with the maximal air flow already obtained in the identification. The reality is that these two variables are not directly proportional, which can explain some of the variations between the measurements and the simulation in the results displayed in the next section (Fig. 4).

Using the DLS control strategy from Definition 1, we aim to regulate each temperature of the experiment in its target interval defined in (3). To ensure that the system encounters several disturbances, including both extremal conditions, we apply the following disturbance schedule, starting when all doors are closed and all lamps are off:

t = 330s, lamps 2

B. Experimental results

The results from that scenario are displayed (in blue) in Fig. 4. The horizontal black lines are the boundaries of the target interval given in (3). As in Fig. 3, the vertical lines represent the switching times in the schedule described above: solid lines when the switched component (door or lamp) is link to the room; dashed otherwise. In addition, the red dashed curves in Fig. 4 correspond to the simulation of this control using the identified model from Section II-B and the same exogenous inputs w and disturbances δ.

In room 3 (bottom left of Fig. 4), when the temperature increases the model is significantly higher than the measurements. This can be explained by the observed limitation of the model, described in Section II-C, when open doors and fans are associated. In the experiment of Fig. 4, room 3 has a higher temperature than its neighbors (2 and 4). With the DLS controller from Definition 1, the fans are running most of the time. When the doors in room 3 are open (between 930 and 2100 seconds and between 3510 and 4350 seconds), the flows from the fans 2 and 4 will go in the opposite direction of the flows at the doors (supposed to go from the hot room 3 to the colder 2 and 4). This phenomenon, which is not modeled in [START_REF] Perez-Lombard | A review on buildings energy consumption information[END_REF], explains the differences between the model and the data in room 3. The last observation we can do on the results in Fig. 4 is that the system is indeed invariant in the target interval chosen in (3). As previously listed in the disturbance schedule, both the hottest and the coldest situations are met for the discrete disturbances δ. However, we can neither ensure that the extremal conditions of the exogenous inputs w were met, nor that they were met at the same time as the extremal values of δ. But we can note that these temperatures mostly intervene in the conduction which has a reduced effect on the system compared to the other heat transfers associated with the heat radiation of the air flows. We can thus almost claim (due to the model uncertainties) that the system is indeed Robust Controlled Invariant in the chosen target interval.

V. CONCLUSION

In this paper, we use a 0-dimensional model to describe the temperature dynamics in a building equipped with UnderFloor Air Distribution. Identification and verification methods are then used to fit the behavior of a smallscale experiment of an UFAD flat. The Robust Controlled Invariance is presented as the ability to maintain the state in a given target interval for any value of the disturbances and this result is greatly simplified by the monotonicity property satisfied by our model. Finally, the invariance criteria are computed using the identified model and confirmed in a control experiment on the small-scale UFAD flat.

The results on Robust Controlled Invariance presented in this paper can be generalized to any system provided it satisfies the monotonicity property and a local control property (each controlled input only affects one state variable). While we do not know a general category of systems satisfying the latter, we can find the monotonicity property in various fields such as molecular biology, chemical reactions, vehicle platooning, . . . Furthermore, this approach does not depend on the feedback control strategy used in the Robust Controlled Invariant interval, which leaves a large degree of freedom for the performance specifications.

Fig. 1 .

 1 Fig. 1. 4-room flat with UFAD

Fig. 3 .

 3 Fig. 3. Comparison between the identified model (dashed, red) and the verification data set (noisy, blue), with the transitions (vertical lines: solid when linked to the room, dashed otherwise)

  and fan 3 on; t = 570s, lamp 3 on, door 1 -4 open; t = 810s, fan 4 on; t = 930s, lamp 3 off, door 1 -2 open; t = 1050s, door 1 -4 closed. This data set gathers the main situations that can happen in the building: conduction alone (rooms 2 and 4); a lamp (room 1); a fan (room 3); an open door (rooms 2 and 4); a lamp and a fan (room 3); a lamp and an open door (room 1); a fan and an open door (room 4).

  and 3 on; t = 930s, doors 1 -2 and 2 -3 open; t = 1530s, lamp 4 on, door 3 -4 open; t = 2010s, lamp 3 off, doors 2 -3 and 3 -4 closed; t = 2500s, lamp 4 off, door 4 -1 open; t = 2910s, lamp 1 on, door 4 -1 closed; t = 3510s, all lamps on, all doors open; t = 4350s, all lamps off, all doors closed.

Fig. 4 .

 4 Fig. 4. Controlled state, using DLS controller, for the experiment (blue, noisy) and the simulation (red, dashed), with the transitions (vertical lines: solid when linked to the room, dashed otherwise)

An extended version of this paper containing the proofs is available at http://hal.archives-ouvertes.fr/docs/00/85/87/15/ PDF/Controllability.pdf

Built at the physics department (UFR PhITEM) of University Joseph Fourier, Grenoble, France.


Benoît Gaüzère
email: benoit.gauzere@ensicaen.fr

Luc Brun
email: luc.brun@ensicaen.fr

Didier Villemin
email: didier.villemin@ensicaen.fr

Graph kernel encoding substituents' relative positioning

Chemoinformatics aims to predict molecular properties using informational methods. Computer science's research fields concerned by this domain are machine learning and graph theory. An interesting approach consists in using graph kernels which allow to combine graph theory and machine learning frameworks. Graph kernels allow to define a similarity measure between molecular graphs corresponding to a scalar product in some Hilbert space. Most of existing graph kernels proposed in chemoinformatics do not allow to explicitly encode cyclic information, hence limiting the efficiency of these approaches. In this paper, we propose to define a cyclic representation encoding the relative positioning of substituents around a cycle. We also propose a graph kernel taking into account this information. This contribution has been tested on three classification problems proposed in chemoinformatics.

I. INTRODUCTION

Chemoinformatics aims to predict molecular properties from their structural similarities. Most of existing methods are based on fingerprints defined as collections of descriptors such as the boiling point, logP, molar refractivity, etc. An alternative strategy consists in computing a set of descriptors directly from the molecular graph G of a molecule. A molecular graph G = (V, E, µ, ν) consists of an unlabeled graph (V, E) encoding the structure of a molecule and two labeling functions µ and ν mapping respectively each node to a label encoding atom's chemical element and each edge to a label identifying a type of atomic bond (single, double, triple or aromatic). Considering this molecular representation, the similarity between molecules can be directly deduced from the similarity of their molecular graphs instead of using a set of descriptors chosen a priori.

Graph kernels can be understood as symmetric graph similarity measures. Using a semi definite positive kernel, the value k(G, G ′ ), where G and G ′ encode two graphs, corresponds to a scalar product between two vectors ψ(G) and ψ(G ′ ) in an Hilbert space. This latter property allows us to combine graph kernels with widely-used machine learning algorithms, such as SVM, thanks to the kernel trick which consists in replacing the scalar product ψ(G), ψ(G ′ ) by the value of k(G, G ′ ) in these algorithms. Graph kernels provide thus a natural connection between molecular graphs and machine learning frameworks. A large family of kernels is based on bags of patterns. These methods consist in extracting a bag of patterns from a graph and deducing the similarity between two graphs from the similarity between their respective bags. The Hilbert space associated to the computed kernel is directly induced by the space associated to the set of extracted patterns. Therefore, the similarity measure encoded by a kernel mostly depends on the patterns extracted from molecular graphs. Most of existing methods are defined on linear patterns [START_REF] Kashima | Kernels for graphs[END_REF]. Such methods are generally associated to a low complexity but they are limited by the lack of expressivity of linear patterns. Nevertheless, in order to encode more structural information, some kernels are based on non-linear acyclic patterns [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF], [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF]. However, kernels defined directly on molecular graphs do not explicitly take into account the cyclic similarities of molecules. For example, treelet kernel [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF] deduces the similarity of two molecules from the similarity of two bags of patterns corresponding to the set of all labeled subtrees composed of at most 6 nodes (figure 4). Obviously, considering a set of acyclic structures extracted from the molecular graph does not allow to encode the cyclic information. Nonetheless, cycles have an important impact on molecules' properties since they reduce atoms' degrees of freedom. Therefore, this information must be taken into account in order to define accurate similarity measures.

Considering this, Horváth proposed in [START_REF] Horváth | Cyclic pattern kernels for predictive graph mining[END_REF] to combine a kernel encoding an acyclic similarity with an intersection kernel defined on the set of simple cycles of a graph. Despite the high complexity of the enumeration of all simple cycles of a graph, this method can be efficiently used when graphs have a low number of cycles. In order to tackle the complexity induced by the enumeration of all simple cycles, Horváth proposed in [START_REF] Horváth | Cyclic pattern kernels revisited[END_REF] to use a subset of simple cycles. This set is first initialized using the set of relevant cycles [START_REF] Vismara | Union of all the minimum cycle bases of a graph[END_REF] of a graph. Then, additional simple cycles are iteratively enumerated by combining relevant cycles and newly discovered cycles. Horváth shown that a low number of iterations is sufficient to obtain results similar to the ones obtained using all simple cycles. In order to encode a finer cyclic information, Gaüzère et al. proposed in [START_REF] Gaüzére | Graph kernels based on relevant patterns and cycle information for chemoinformatics[END_REF] to deduce the cyclic similarity of two molecules from the similarity of their relevant cycle graphs [START_REF] Vismara | Reconnaissance et reprsentation d'lments structuraux pour la description d'objets complexes. application l'laboration de stratgies de synthse en chimie organique[END_REF]. This approach allows to encode adjacency relationships between cycles in addition to isolated cycles as proposed by Horváth. However, this molecular representation does not include any acyclic information. Such information may be encoded by the relevant cycle hypergraph as defined in [START_REF] Gaüzère | Relevant cycle hypergraph representation for molecules[END_REF]. This molecular representation defines a global representation of a molecule while explicitly encoding cyclic information. Treelet kernel can then be adapted in order to define a similarity measure based on this molecular representation.

However, these approaches can not encode the relative positioning of the set of substituents around a given cycle. This information is particularly important in chemoinformatics since it characterizes the spatial configuration of a molecule which in turn strongly influences many biological properties. Therefore, we propose in this paper to define a graph kernel encoding this information. To this end, we first review in section II the different approaches to encode cyclic information into graph kernels. Then, we propose in section III-A a representation, called augmented cycles, which allows to encode the relative positioning of substituents for each relevant cycle of a molecule. Then, we propose in section III-B a kernel based on the treelet kernel which includes the information encoded by augmented cycles. This approach relies on increasing the information carried by a treelet. Finally, section IV demonstrates the insight of our contribution.

II. MOLECULAR CYCLIC INFORMATION

A. Encoding cyclic information

Considering a molecular graph G = (V, E, µ, ν), a simple cycle C is defined as a connected subgraph C = (V ′ , E ′ , µ, ν) of G where each vertex v ∈ V ′ has a degree equal to 2. In order to explicitly encode the cyclic information of a molecule, Horváth proposed to consider all simple cycles [START_REF] Horváth | Cyclic pattern kernels for predictive graph mining[END_REF]. However, this first approach has an high complexity which restricts its use to molecules having a low number of cycles. A less complex approach consists in considering a subset of simple cycles such as the set of relevant cycles [START_REF] Vismara | Union of all the minimum cycle bases of a graph[END_REF]. In order to define the notion of relevant cycle, let us first introduce the notion of cyclic vector space. Each cycle C ⊆ G can be represented as a vector C ∈ {0, 1} |E| where C i equals 1 if i is an edge of C and 0 otherwise. The set of vectors encoding the cycles of G defines a vector where the addition of two cycles C and C ′ corresponds to a XOR bitwise [START_REF] Vismara | Union of all the minimum cycle bases of a graph[END_REF]. The set of relevant cycles of a graph G, denoted C R (G), is defined by the union of all bases of the vector space of minimum length, the length of a base being defined as the sum of lengths of its cycles. Therefore, a relevant cycle can not be encoded by a combination of cycles having a lower size. Intuitively, the set of relevant cycles corresponds to a canonical set of elementary cycles which allows to retrieve all cycles included in a molecule. This set of cycles can be enumerated in polynomial time with the number of nodes of the graph. Therefore, enumerating the set of relevant cycles instead of the set of all simple cycles allows to obtain a reduced complexity.

The first approach which consists in encoding a molecule by a set of cycles allows to explicitly characterize an important part of the cyclic information. However, this approach is only based on isolated cycles which do not allow to identify the cyclic system of a molecule. This system is defined by the relevant cycle graph (figure 1(c)), initially proposed by [START_REF] Vismara | Reconnaissance et reprsentation d'lments structuraux pour la description d'objets complexes. application l'laboration de stratgies de synthse en chimie organique[END_REF] and modified by [START_REF] Gaüzére | Graph kernels based on relevant patterns and cycle information for chemoinformatics[END_REF]. The relevant cycle graph defines a graph representation which encodes the set of relevant cycles and their adjacency relationships:

Definition 1 (Relevant cycle graph). Let us consider a graph G = (V, E, µ, ν). The relevant cycle graph associated to G is defined by a graph G C (G) = (V C , E C , µ C , ν C ) such that : • Each vertex v ∈ V C corresponds to a relevant cycle c v ∈ C R (G); • two nodes (u, v) ∈ V 2
C are connected by an edge (u, v) ∈ E C if the corresponding cycles c u and c v share at least one vertex in G.

The relevant cycle graph allows to encode the cyclic system of a given molecule by mapping each relevant cycle to a node.

Each node v ∈ V C is then labeled by a label's sequence µ(v) identifying the set of atoms and atomic bonds included within the cycle c v encoded by v. Similarly, each edge e ∈ E C is labeled by a label's sequence identifying the set of nodes and edges common to the two cycles incident to e. Including cycles' adjacency relationships allows to encode a finer cyclic information than the one encoded by a set of isolated cycles. However, we can note that this representation encodes the cyclic system of a molecule by removing all acyclic parts. Therefore, an important part of the molecule, such as adjacency relationships between cycles and acyclic parts, is not encoded.

In order to encode a finer cyclic information, Gaüzère et al. proposed in [START_REF] Gaüzère | Relevant cycle hypergraph representation for molecules[END_REF] a molecular representation, called relevant cycle hypergraph (figure 1(d)), which allows to add acyclic parts to the representation defined by the relevant cyclic graph. A relevant cycle hypergraph

H CH (G) = (V CH , E CH , µ CH , ν CH ) can be defined from a molecular graph G = (V, E, µ, ν) and its relevant cycle graph G C (G) = (V C , E C , µ C , ν C ).
The set of nodes V C can be associated to the set V C R ⊆ V which corresponds to the set of nodes of G included within a cycle. Similarly, E C R ⊆ E encodes the set of edges included in at least one relevant cycle. The complement of V C R and E C R in V and E corresponds thus to the acyclic parts of the molecular graph G which are not encoded by its relevant cycle graph. In order to define a complete molecular representation including both cyclic and acyclic parts, the set of nodes V CH of the relevant cycle hypergraph is defined by the union of two subsets. A first subset V C corresponding to the set of relevant cycles and thus encoding an explicit cyclic information, and a second subset V -V C R corresponding to the set of atoms not included within a cycle. Considering the set of vertices V CH , [START_REF] Gaüzère | Relevant cycle hypergraph representation for molecules[END_REF] defines a function p :

V → P(V CH ) as p(u) = {u} if u / ∈ V C R and p(u) = {c ∈ V C | u ∈ V (c)} if not
, where V (c) corresponds to the set of nodes included within the relevant cycle encoded by c. In other words, if u ∈ V C R , p(u) is defined as all nodes v ∈ V C whose associated relevant cycle includes u. This function p encodes thus the print of vertex v ∈ V on V CH . In the same way as for vertices, the set of hyperedges E CH is composed of two subsets:

1) A set of edges E e CH composed of: • edges between relevant cycle vertices, corresponding to the set of edges E C included within the relevant cycle graph,

• edges e = (p(u), p(v)) such that (u, v) ∈ E -E C R , |p(u)| = 1 and |p(v)| = 1.
This set of edges corresponds to edges of the molecular graph G connecting two acyclic atoms or two relevant cycles or a relevant cycle to an acyclic part of G, 2) and a set

E h CH composed of oriented hyperedges e = (p(u), p(v)) such that (u, v) ∈ E -E C R , |p(u)| > 1 or |p(v)| > 1.
This set of hyperedges corresponds to special cases where an edge connects at least two distinct relevant cycles to another part of the molecule. Such an edge connects two sets of vertices s 1 = p(u) and s 2 = p(v) and is thus encoded by an oriented hyperedge e = (s 1 , s 2 ) ∈ E h CH . Thanks to the use of hyperedges, the relevant cycle hypergraph allows to encode adjacency relationships between cycles and acyclic parts. These relationships may be useful to identify

B. Kernels computing a cyclic similarity

The different representations presented in the previous section allow to encode different levels of cyclic information. These different representations can thus be used to define different similarity measures. First graph kernels introduced in chemoinformatics consists in extracting a bag of patterns composed of acyclic sub-structures directly from the original molecular graph [START_REF] Kashima | Marginalized Kernels Between Labeled Graphs[END_REF], [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF], [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF]. The similarity of the two molecules being compared is then deduced from the similarity of the two extracted bags of patterns. For example, the treelet kernel [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF] is based on a bag of patterns T (G), called treelets, defined as all labeled subtrees composed of at most 6 nodes (figure 4). Given two bags of treelets T (G) and T (G ′ ) extracted from two molecular graphs G and G ′ , the similarity between G and G ′ is encoded by a kernel k(G, G ′ ) defined as:

k(G, G ′ ) = t∈T (G)∩T (G ′ ) k ′ (f t (G), (f t (G ′ )) (1)
where each term of the sum corresponds to a subkernel k ′ encoding a similarity measure between the number of occurrences of treelet t in G (encoded by f t (G)) and G ′ (encoded by f t (G ′ )). This kernel allows to take into account most of the information encoded within a molecular graph by considering labeled and nonlinear patterns. However, since the similarity is computed from acyclic structures extracted from the molecular graph, such kernels can only partially and implicitly encode the cyclic information included in molecular graphs.

In order to compute a graph kernel taking into account the cyclic similarity of molecules, Horváth proposed in [START_REF] Horváth | Cyclic pattern kernels revisited[END_REF] to define a kernel on a set of cycles extracted from both graphs to be compared. This kernel is defined as an intersection kernel on the set of relevant cycles, hence computing the number of common relevant cycles of two molecular graphs. Acyclic similarity is then included by a second kernel computing the number of common subtrees extracted from the bridges of the two molecules to be compared. While providing a first kernel including cyclic similarity, this approach is only based on a set of isolated cycles and can not take into account adjacency relationships between cycles.

An alternative graph kernel, defined in [START_REF] Gaüzére | Graph kernels based on relevant patterns and cycle information for chemoinformatics[END_REF], consists in applying the treelet kernel on relevant cycle graphs instead of original molecular graphs in order to compute a cyclic system similarity measure. This kernel consists thus in extracting the set of treelets from the relevant cycle graphs associated to the molecules. This set of patterns allows us to encode most of the adjacency relationships encoded within relevant cycle graphs and thus the cyclic system of molecules. Considering this bag of patterns, similarity between relevant cycle graphs is defined analogously as the case of the original treelet kernel:

k C (G, G ′ ) = t∈T (G C (G))∩ T (G C (G ′ )) k ′ (f t (G C (G)), (f t (G C (G ′ ))) (2)
The main difference with the kernel defined in equation 1 is that treelets are extracted from the relevant cycle graph, hence encoding an explicit cyclic information and cycle's adjacency relationships.

Similarly to the method proposed by Horváth, the cyclic system similarity measure can be combined with an acyclic similarity measure, such as the treelet kernel applied on the molecular graph, in order to define a global similarity measure between molecules including both cyclic and acyclic similarities. Despite the fact that this approach leads to good results on experiments involving cyclic molecules, this representation, as the one proposed by Horváth [START_REF] Horváth | Cyclic pattern kernels revisited[END_REF], separates cyclic and acyclic information by using two different molecular representations. Then, global similarity between molecules is computed by the combination of two distinct similarity measures, each of them being applied on one representation. This separation induces thus a loss of adjacency relationships between cyclic and acyclic parts in the similarity measure.

Adjacency relationships between a cycle and its substituents can be encoded within the relevant cycle hypergraph (section II-A). In the same way as the adaptation of the treelet

R 2 R 1 R 3 R 2 R 1 R 3
Fig. 2. Two cycles and their respective sets of substituents having a different relative positioning.

kernel to relevant cycle graph comparisons, [START_REF] Gaüzère | Relevant cycle hypergraph representation for molecules[END_REF] proposed to adapt the treelet kernel to the comparison of relevant cycle hypergraphs. The set of treelets extracted from the relevant cycle hypergraph then encodes adjacency relationships between a cycle and its substituents in addition to the information encoded by the relevant cycle graph. This additional information allows to define a finer similarity measure which leads to an increasing accuracy on several molecular property prediction problems [START_REF] Gaüzère | Relevant cycle hypergraph representation for molecules[END_REF]. This last contribution defines thus a global similarity measure based on an unique molecular representation encoding the whole molecule and an explicit cyclic information.

III. KERNEL BASED ON AUGMENTED CYCLES

A. Definition of augmented cycles

Previous approaches allow to encode an important part of the cyclic information by defining different molecular representations encoding different levels of cyclic information. While these representations encode most of the cyclic information, they do not allow to distinguish the relative positioning of substituents around a given cycle. Indeed, the two patterns displayed in figure 2 are composed of a same cycle, denoted C 1 , surrounded by a same set of substituents R 1 , R 2 and R 3 . However, we can note that these substituents are not positioned in the same way around C 1 . Using a molecular representation encoding cycles adjacency relationships such as the relevant cycle hypergraph, these two patterns can not be distinguished since their representations are isomorphic. However, a molecular cycle is a planar structure and the relative positioning of its substituents characterizes the spatial configuration of the whole pattern. Since different spatial configurations may lead to different molecular properties, the relative positioning of the substituents should be taken into account in order to define an accurate similarity measure.

Let us consider a molecular graph G and its associated relevant cycle hypergraph H CH (G). Given a relevant cycle C

of G encoded by c in H CH (G), the set of substituents of C is defined as the neighbourhood Γ(c) of c. For example, N and C 2 are two substituents of C 3 in figure 1(a). Since a cycle is a planar structure, we can define a coherent orientation on its boundary and hence a cyclic order on the neighbourhood Γ(c) of each vertex c ∈ V CH encoding a cycle.

Our basic idea to encode the relative positioning of the substituents of a cycle C consists in associating each substituent to a position on the boundary of C.

Given the node c ∈ V CH R 1 R 3 R 3 R 2 R 1
Fig. 3. Encoding of positions around a cycle. associated to C, edges encoding the adjacency between c and nodes belonging to Γ(c) may either correspond to an edge of the original molecular graph G or to an edge of the relevant cycle graph encoding a cycle's adjacency relationship. The position of substituents incident to edges of the original molecular graph may be easily defined since these edges are incident to a single atom of C (e.g. (C 3 , N ) in figure 1(a)). However, edges of the relevant cycle graph correspond to a common path between two cycles (e.g. (C 3 , C 2 ) in figure 1(a)). In such cases we fix the position of the substituents incident to such edges to the middle of the shared path. Since this last position may not correspond to a cycle's atom but to the middle of an edge, the number of possible positions is equal to twice the number of edges of the cycle (figure 3). Using an arbitrary origin and orientation on C, these positions may be efficiently encoded by a function pos C : Γ(c) → N counting the number of half-edges between our origin and any substituent.

The relative positioning of two substituents according to C may be encoded by their angle. The angle between two substituents of C is defined as the number of half-edges separating them using the orientation fixed on C. This number of half-edges is directly retrieved from the substituents' positions encoded by the function pos C : Definition 2 (Angle between substituents). Using the above notations and hypothesis, given two nodes s 1 and s 2 of Γ(c), the angle between s 1 and s 2 according to c is defined by:

ang c (s 1 , s 2 ) = (pos C (s 2 ) -pos C (s 1 )) mod (2|C|) (3)
where |C| denotes the number of edges of the cycle C of G.

We can note that equation 3 does not depend on the origin but only on the orientation defined on C. The combination of a cycle c and a function is defined by

α i c = ang c (s i , s i+1 ), ∀i ∈ {1, . . . , n -1}. (4)
For example in figure 3, the ordering σ = (R 1 , R 3 , C 1 ) of Γ(C 2 ) corresponds to the sequence of angle's [START_REF] Vismara | Union of all the minimum cycle bases of a graph[END_REF][START_REF] Horváth | Cyclic pattern kernels revisited[END_REF]. Choosing the alternative orientation on C 2 or a different origin leads to another ordering associated to a different angle's sequence. Let us consider a sequence of substituents σ = (s 1 , . . . , s n ) and its associated sequence angles(σ). We can define a coherent function pos ′ c from both sequences by setting pos ′ c (s 1 ) = 0 and pos ′ c (s i ) =

i-1 j=1 α j c mod 2|C| for each i ∈ {2, . . . , n}. Such a function corresponds thus simply to a shift of the original function pos c where pos ′ c (s 1 ) is set to the origin.

B. Augmented treelets

A treelet is defined by a subtree T = (V T , E T , µ T , ν T ) extracted from a graph G = (V, E, µ, ν). Each treelet extracted from a graph is canonically identified by a function key(T ) defined as the combination of an index encoding its structure (V T , E T ) (figure 4) and a code encoding the labels associated to each node and edge of T given a particular treelet's structure. This function is defined such as two isomorphic treelets T and T ′ are associated to the same key, i.e. T ≃ T ′ ⇔ key(T ) = key(T ′ ) [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF]. The code identifying the labeling information of a treelet is defined as a node and edges' label sequence. When applied to relevant cycle hypergraph, each node of an extracted treelet may encode either an atom or a cycle and this node is labeled by the labeling function µ CH defined on the relevant cycle hypergraph (section II-A). Two treelets encoding the two configurations shown in figure 2 are associated to a same structure and the same set of node labels. Therefore, these two treelets are isomorphic since the label of a node encoding a relevant cycle does not include the relative positioning of its adjacent nodes.

In order to encode this information, we propose to introduce the notion of augmented treelets by encoding the relative positioning of substituents into the notion of treelet. Given a treelet T = (V T , E T , µ T , ν T ) let us consider an order {v 1 , . . . , v n } on V T induced by our key construction scheme [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF]. This order is unique up to the automorphisms of T and induces an order on the neighborhood Γ(c) of each vertex c ∈ T encoding a cycle. Let us denote by Σ c (T ) the sequence of subtrees incident to c, the order being deduced from the one defined on the neighborhood Γ(c). Each automorphism τ i ∈ Aut(T ) of T induces a new numbering of the vertices of T and hence a new order on the neighborhood of vertices c corresponding to a relevant cycle. Let us denote by Σ c (τ i (T )) the new ordering of Σ c (T ) induced by τ i . By definition of automorphism if Σ c (T ) = t 1 . . . .t p and Σ c (τ i (T )) = t ′ 1 . . . .t ′ p we have t i ≃ t ′ i for all i ∈ {1, . . . , p}. In other words, if a tree automorphism maps two branches of a tree one onto the other, these branches should be isomorphic. However, these trees may be attached differently on the cycle encoded by c and thus correspond to different sequences of angles. We thus define the canonical sequence of angles angles ⋆ (c) for each node c ∈ V T encoding a relevant cycle as:

angles ⋆ (c) = min τi∈Aut(T ) angles(σ c (τ i (T ))) (5)
where σ c (τ i (T )) corresponds to the order on Γ(c) induced by Σ c (τ i (T )). Let us consider a canonical sequence angles ⋆ (c) = α 1 c . . . α p c and the initial order of subtrees Σ c (T ) = t 1 . . . .t p defined on Γ(c) in T . Each angle α i c encodes an angle between a subtree isomorphic to t i and a subtree isomorphic to t i+1 . Hence, the sequence angles ⋆ (c) encodes the angles between the different subtrees of T incident to c. Using Section II-A, the position of each subtree on the cycle associated to c may be retrieved up to a shift.

In order to include the relative positioning of substituents into our augmented treelets, we define a new labeling function µ ′ T (v) as the concatenation of µ T (v) and angles ⋆ (v). The key associated to an augmented treelet is then computed in the same way as in [START_REF] Gaüzère | Two new graphs kernels in chemoinformatics[END_REF] according to the augmented treelet T a = (V T , E T , µ ′ T , ν T ). Two augmented treelets with a same key are thus isomorphic. Hence two vertices c and c ′ corresponding to relevant cycles and mapped one onto the other have the same label and hence the same set of connected subtrees and the same sequence of angles. According to our previous considerations, the positions of each subtree may be retrieved without ambiguities from the sequence of angles in both graphs. Both vertices encode thus cycles with the same connected subtrees separated by the same angles. Considering such a labeling function, the two configurations depicted in figure 2 may thus be distinguished. The set of augmented treelets allows us to extract a set of patterns encoding cyclic, acyclic and cyclic/acyclic adjacency relationships. In addition, our bag of treelets allows to distinguish configurations having a different subtituents' relative positioning, hence encoding a finer cyclic information than the set of treelets extracted from the relevant cycle hypergraph.

IV. EXPERIMENTS

Our extension has been tested on three classification experiments involving molecular graphs composed of cycles. Results obtained by different kernels using SVM are displayed in table I. The first experiment, denoted PAH, consists in 94 molecules only composed of carbons. This experiment aims to predict if a molecule is carcinogenic or not. Results displayed in table I correspond to the best results obtained by using a 10-fold cross classification over a grid of parameters including both machine learning and kernel parameters. The second experiment, denoted Mutag, is taken from [START_REF] Riesen | Iam graph database repository for graph based pattern recognition and machine learning[END_REF] and is composed of 4 337 molecules identified as mutagen or nonmutagen. This set of molecules is divided into 3 distinct subsets: 1 500 molecules are used to train our model, 500 as a validation set used to tune the parameters and the accuracy is evaluated on the remaining 2 337 molecules, hence composing the testset. Finally, the last experiment corresponds to the PTC [START_REF] Toivonen | Statistical evaluation of the predictive toxicology challenge 20002001[END_REF] dataset which consists in predicting the toxicity of different molecules on four classes of animals: male mice (MM), female mice (FM), male rats (MR) and female rats (FR). Each class of animal is composed of ten trainsets and ten testsets and, similarly to the first experiment, parameters are tuned over the testset for each trainset for all methods. Table I displays the mean accuracy over the ten testsets of each class of animals.

In order to evaluate our contribution, we compare it with different kernels. Treelet kernel (Line 1) is a kernel based on a bag of patterns (section II-B) which do not explicitly encode cyclic information. The second kernel (Line 2) is computed by applying a gaussian kernel on an approximate edit distance [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] which thus do not explicitly encode cyclic information. Note that this kernel may need some regularization step in order to define a positive definite kernel. Next kernels explicitly include different levels of cyclic information (sections II-A and II-B). Line 3 corresponds to the kernel defined by Horváth in [START_REF] Horváth | Cyclic pattern kernels revisited[END_REF] which only compares the set of relevant cycles extracted from molecules to be compared. The kernel based on the relevant cycle graph (Line 4) allows to encode cycle adjacency relationships and the one based on the relevant cycle hypergraph (Line 5) encodes in addition the acyclic parts of molecules. Finally, line 6 corresponds to our contribution which encodes the position of each substituent along the boundary of each relevant cycle.

First, we can observe that the similarity measures which do not include explicitly the cyclic information are generally outperformed by the treelet kernel on relevant cycle hypergraph (line 5) which explicitly includes the cyclic information. This observation allows to highlight the relevancy of explicitly encoding cyclic information. However, we can note that considering only the cyclic system (Horváth and relevant cycle graph kernels) do not allow to obtain a good accuracy. Second, we can note that finer the cyclic information is encoded, better are the results. This observation is systematically observed on PTC and Mutag dataset whereas we observe a lower accuracy for the relevant cycle hypergraph on PAH dataset. This can be explained by the fact that the molecules included within PAH dataset are mainly cyclic. In this case, the addition of the acyclic information without including the relative positioning of substituents induces some irrelevant features. Conversely, adding acyclic parts together with the relative positioning of substituents lead to the best results on this dataset. In addition, we can note that our contribution allows to systematically outperform the results obtained by other kernels, hence showing the usefulness of our contribution. From a computational point of view, computing a Gram matrix on 150 molecules of kernel including the relative positioning of substituents requires about 1 s whereas it requires about 0.5 s for the kernel on relevant cycle hypergraph. These experiments have been performed thanks to the resources of computer center CRIHAN.

V. CONCLUSION

In this paper, we proposed a new cycle representation, called augmented cycles, which allows to encode cyclic information of molecules in a finer way. Our contribution consists in including the relative positioning of substituents around a cycle. This new representation can be used in molecular representations where the cyclic information is explicitly encoded, such as the relevant cycle hypergraph. In order to use this additional information, we also proposed to adapt the treelet kernel to encode the relative positioning of substituents by defining the augmented treelets. The different experiments show the relevancy of this additional information to solve classification problems in chemoinformatics. Our contribution defines only locally this orientation for a given cycle. Further works will aim to study the insight of using a coherent orientation, when it exists, between adjacent cycles. To that end, we will have to determine for each cycle its maximal planar connected component.

Fig. 1 .

 1 Fig. 1. Different molecular representations encoding explicitly the cyclic information.

  ang c : Γ(c) 2 → N defines the notion of augmented cycle. Let us denote by π(Γ(c)) the set of possible orderings of Γ(c). Using equation 3, any ordering σ ∈ π(Γ(c)) may be associated to a sequence of angles, each angle being defined between two successive elements of σ. Such a sequence is encoded by the function angles(σ): Definition 3 (Angle's sequence). Given a cycle encoded by c ∈ V CH and an ordering σ = (s 1 , . . . , s n ) on Γ(c), the angle's sequence angles(σ) = α 1 c . . . . α n-1 c

Fig. 4 .

 4 Fig.[START_REF] Horváth | Cyclic pattern kernels for predictive graph mining[END_REF]. Augmented treelets and the set of angles (α 1 , . . . , αn) associated to each potential relevant cycle having at least one pair of substituents.

  some particular substructures which have an important influence on some molecular properties. Note that, by adding the acyclic parts to the cycles encoded by the relevant cycle graph, each node and each edge of the original molecular graph is encoded within the relevant cycle hypergraph, either by the node or the edge itself or by a node encoding a relevant cycle. Therefore, this representation corresponds to a complete molecular representation and acyclic, cyclic and acyclic/cyclic adjacency relationships are included into the relevant cycle hypergraph. In addition, since this representation is based on the relevant cycle graph, the cyclic information is explicitly encoded.

TABLE I .

 I ACCURACY ON MUTAG, PAH AND PTC DATASETS.

	Method	Mutag	PAH		PTC
				MM	FM	MR	FR
	Treelet Kernel [3]	77.1	71.3	61.9	58.7	60.8	60.4
	Graph edit distance [14]	71.5	72	66.4	60.7	56.4	66.7
	Cycles [5]	64.5	63.0	62.2	59.3	58.7	65.0
	Relevant cycle graph [7]	67.8	77.7	62.8	60.2	59.0	66.0
	Relevant cycle hypergraph [9]	78.3	76.3	64.6	64.2	60.2	66.4
	Augmented cycle	80.2	80.7	67.9	64.8	61.3	68.7