Khaled Metwally

Serge Mensah

Guillaume Baffou

Isosbestic Thermoplasmonic Nanostructures

Keywords: plasmonics, nanoparticles, photothermal, thermoplasmonics

This article introduces the concept of photothermal isosbesticity in plasmonics. In analogy with absorbance spectros-copy, this concept designates nanostructures that feature an invariance of their temperature increase upon varying the illumination polarization angle. We show that nontrivial (i.e., non-centrosymmetric) isosbestic nanostructures exist and prove valuable when the optical near-field intensity remains, on the contrary, highly dependent on the illumination polarization. The concept is introduced with the case of a sphere-dimer, where the conditions for isosbesticity can be derived analytically. The cases of a spheroid and a disc-dimer are also studied in order to draw a general theory and explain how isosbesticity conditions can be obtained from the visible to the infrared range. Nontrivial isosbestic plasmonic nanostructures represent powerful systems to elucidate the origin (thermal or optical) of mechanisms involved in plasmonics-assisted nanochemistry, liquid-gas phase transition, or heat-assisted magnetic recording.

H eat generation and near-field enhancement have been the main functionalities of metal nanoparticles for applications in nanoplasmonics. On the one hand, the gigantic optical nearfield enhancements (NFEs) of metal nanotips or nanogaps prove valuable for applications in nanochemistry, [START_REF] Baffou | Nanoplasmonics for Chemistry[END_REF] SERS, or heatassisted magnetic recording (HAMR). [START_REF] Stipe | Magnetic Recording at 1.5 Pb m -2 Using an Integrated Plasmonic Antenna[END_REF][START_REF] Pan | Heat-Assisted Magnetic Recording[END_REF][START_REF] Challener | Heat-Assisted Magnetic Recording by a Near-Field Transducer with Efficient Optical Energy Transfer[END_REF][START_REF] Kryder | Heat Assisted Magnetic Recording[END_REF][START_REF] Zhang | The Use of Plasmonic Resonances in Thermally Assisted Magnetic Recording[END_REF][START_REF] Matsumoto | Thermally Assisted Magnetic Recording on a Bit-Patterned Medium by Using a Near-Field Optical Head with a Beaked Metallic Plate[END_REF] On the other hand, heat generation (HG) from gold nanoparticles is at the basis of numerous promising biomedical applications, [START_REF] Baffou | Thermo-plasmonics: Using Metallic Nanostructures as Nano-Sources of Heat[END_REF] namely, photothermal cancer therapy, [START_REF] Jain | Au Nanoparticles Target Cancer[END_REF][START_REF] Huang | Gold Nanoparticles: Interesting Optical Properties and Recent Applications in Cancer Diagnostic and Therapy[END_REF][START_REF] Gobin | Near-Infrared Resonant Nanoshells for Combined Optical Imaging and Photothermal Cancer Therapy[END_REF][START_REF] Gu | Targeted Nanoparticles for Cancer Therapy[END_REF] drug and gene delivery, [START_REF] Doane | The Unique Role of Nanoparticles in Nanomedicine: Imaging, Drug Delivery and Therapy[END_REF][START_REF] Levy | Gold Nanoparticles Delivery in Mammalian Live Cells: a Critical Review[END_REF][START_REF] Ghosh | Gold Nanoparticles in Delivery Applications[END_REF][START_REF] Han | Drug and Gene Delivery Using Gold Nanoparticles[END_REF] or remote surgery. [START_REF] Kharlamov | Silica-gold Nanoparticles for Atheroprotective Management of Plaques: Results of the NANOM-FIM Trial[END_REF] In this article, the concept of photothermal isosbesticity in nanoplasmonics is introduced. Isosbestic nanostructures are defined by an HG that is independent of the illumination polarization. The interest in isosbestic nanostructures becomes manifest when considering nanoparticles with no rotational invariance along the direction of light propagation (such as rods, dimers, spheroids, etc.), the near-field of which remains highly polarization dependent. The illustrative textbook case of a plasmonic sphere-dimer is treated in the first part, which is intended to give a physical picture of the concept with an analytical description of the dimer as interacting dipoles. Then, the case of a spheroid is also treated in order to propose a unified description of isosbesticity in plasmonics. These results are also aimed at giving the recipe for adjusting the isosbestic wavelength at will, from the visible range to the infrared. The final part of the article is devoted to explain how isosbestic structures can help identify the mechanisms (NFE or HG) in various plasmonics applications.

■ DEFINITION OF PHOTOTHERMAL ISOSBESTICITY

Let us consider an optical plane wave traveling in a dielectric medium of refractive index n s (n s = 1.33 in this work) and characterized by a wave vector k. An orthonormal frame (O, x, y, z) is oriented so that k is aligned with the (Oz) axis:

= π λ k u n z 2 s 0 .
In this case, the complex amplitude E 0 of the electric field of an incoming linearly polarized light reads E 0 (r)e -iωt = E 0 (cos θu x + sin θu y )e ik•r-iωt where ω = k/n s c and θ defines the light polarization. In this article, all complex number quantities are underlined (such as E 0 ). The irradiance (power per unit area) I of this plane wave reads thus I = n s cε 0 E 0 2 /2. Let us consider now an absorbing nanostructure within this medium, characterized by an absorption cross section σ abs such that the power absorbed by the nanostructure reads P = σ abs I. In this work, an absorbing structure is said to be isosbestic relative to the wave vector k if its absorption cross section is independent of θ. This adjective is borrowed from the field of molecular absorbance spectroscopy, where an isosbestic point refers to the wavelength for which the absorbance of a sample under chemical transformation is independent of a parameter, usually time (because reactants and products have the same molar absorptivity). Here, polarization rather than time is the parameter of interest, acting on the absorbance of the system. In a good approximation, the direct consequence of an invariance of the absorption cross section is an invariance of the temperature increase as well. A sphere is obviously isosbestic for any wave vector due to its point symmetry, but this case is not of particular interest. The response of any system endowed with a 3-fold symmetry (e.g., equilateral triangle) is also polarizationindependent by nature. We coin these systems trivial isosbestic structures. Nontrivial isosbestic systems, which rather consist of nanoparticles that have no axial symmetry along the direction of light propagation (like rods or dimers), are much more relevant for several applications (as explained in the last section of this article).

■ THE CASE OF A SPHERE DIMER

As an introductory nontrivial isosbestic structure, a dimer consisting of two neighboring nanospheres is considered herein. This example is ideal to introduce the concept of isosbesticity because it can be treated analytically in the dipolar approximation and it easily gives a physical picture of the underlying physics of isosbesticity in plasmonics. In this section, the sphere-dimer will be composed of two spheres aligned along the (Ox) axis, in positions r 1 = (0, 0, 0) and r 2 =(d, 0, 0). They have a radius a, they are separated by a gap g, and their center-to-center distance

is d = |r 1 -r 2 | =2a + g.
Numerical Simulations. Prior to presenting the analytical formalism, Figure 1 plots the results of numerical simulations for a sphere-dimer, in order to illustrate the concept. In this whole study, optical simulations have been conducted using the boundary element method (BEM) [START_REF] Hohenester | MNPBEM -A Matlab Toolbox For the Simulation of Plasmonic Nanoparticles[END_REF][START_REF] García De Abajo | Retarded Field Calculation of Electron Energy Loss in Inhomogeneous Dielectrics[END_REF] and thermal simulations using a Green's function approach, [START_REF] Baffou | Thermoplasmonics Modeling: a Green's Function Approach[END_REF] considering a water thermal conductivity of 0.6 W•m -1 •K -1 . The thermal conductivity of the nanoparticles was considered as infinite, which naturally results in uniform temperature distributions within the nanoparticles. This assumption is justified by the much larger thermal conductivity of gold compared to water (500 times larger). This assumption is usually valid in nanoplasmonics, no matter the size of the nanoparticles, provided the nanoparticles are not partially illuminated as shown in ref 21. Figure 1a plots the absorption cross section of a sphere-dimer structure in the case a = 30 nm and g = 10 nm, for perpendicular and parallel polarizations. An invariance of the absorption cross section is observed for λ 0 = 531 nm. This is the isosbestic wavelength according to the above-mentioned definition. The fact that all the cross section line shapes at different angles meet in a single point for some wavelength is not surprising. Any incoming linearly polarized light can be decomposed into two light beams crossedpolarized along the two main axes of the structure. So if the absorption cross sections for the transverse and longitudinal modes are equal, the absorption cross section naturally does not vary by changing the polarization angle. In other words, it suffices to find crossing points of the transverse and longitudinal spectra to make sure that the response of the structure is polarizationangle invariant. Note that this result illustrates a simple procedure to determine the isosbestic wavelength(s) of a plasmonic structure experimentally or numerically: It suffices to measure or compute the absorption cross section spectra of a given structure for the two orthogonal polarizations and to look for possible crossing points. At this isosbestic wavelength of 531 nm, the evolution of the absorption cross section and of the maximum NFE as a function of the illumination angle θ have been plotted in Figure 1b. This figure evidences that for this specific wavelength the energy absorption (and thus the temperature increase) is not dependent on the incident light polarization, while the NFE is highly polarization-dependent. The distributions of the temperature increase and near-field intensity for the two orthogonal polarizations are represented in Figure 1c,d. While the temperature maps remain identical, the optical near-field distributions are different in amplitude and shape.

Let us give some explanation about the physical origin of isosbesticity in the case of a sphere-dimer. The near-field enhancement of a dimer is bound to be highly polarizationdependent, no matter the wavelength, because of the presence of a gap. Within a gap, because of the proximity of the accumulated charges on both sides of the gap, the electric field is necessarily gigantic when the polarization is along the dimer axis, even out of resonance, which naturally implies a strong polarization dependence of the NFE. However, since the heat generation has nothing to do with the electric field within the gap, but rather with the electric field inside the metal, the presence of a gap does not prevent having a constant heat generation as a function of the light polarization, which opens the path for isosbesticity.

Analytical Formalism in the Dipolar Approximation. Determining the optical response of a dimer of nanospheres considered as two interacting dipoles is a self-consistent problem, governed by this set of equations:

α α ̲ =+ ̲ ̲ ̲ ̲ =+ ̲ ̲ ̲ ω ω ⎪ ⎪ ⎧ ⎨ ⎩ EE r r r E EE r r r E () (, ) () (, ) 1 ext 01 1 2 2 ext 2 ext 02 2 1 1 ext (1)
where E j ext is the electric field experienced by sphere j ∈ {1, 2}, generated by its environment. It has two origins, the incoming plane wave E 0 and the electric field generated by the neighboring nanoparticle. αE j ext = p j is the polarization vector of sphere j where α is the complex polarizability of the sphere. ̲ ω rr (, ) mn is the field propagator; it has a closed-form expression in a uniform medium given in the general case by [START_REF] Baffou | Luminescence induite par microscopie àeffet tunnel et etude des proprieteśeĺectroniques, chimiques et optiques de la surface de carbure de silicium 6H-SiC[END_REF]

πε

̲ =- - + || k rr R R R (, ) ( () i k() () ) e 4 mn k R 2 12 3 i 0 where =- = •- = •- = •- R R R R R R Rr r R RR R RR R RR () () 3
() 3 mn 1 T2 3 2 T2 4 3 T2 5
In this series of equations, R T means the transpose of the column vector R, R is the norm of R, is the identity matrix, R ()

1
represents the propagative field, while R ()

2
and R ()

3
represent the optical near-field. In the following, 1 will be neglected, although closed-form expressions of a dimer response can also be obtained without this approximation. [START_REF] Rolly | Metallic Dimers: When Bonding Transverse modes Shine Light[END_REF] This approximation is valid if kd ≪ 1. ̲ rr (, ) mn is dependent only on the column vector R and is invariant by the inversion of m and n. Thanks to these relations, system (1) can be recast into a self-consistent matrix equation involving only the unknowns E j ext . This system can be simplified. First, as the two nanoparticles are located in the (Oxy) plane, the z coordinate can be skipped and the dimension of the matrix equation reduced accordingly. Second, because of the symmetry of the sphere-dimer and of the illumination, the electric field experienced by each sphere must be identical:

E 1 ext = E 2
ext =̂E ext , as well as E 0 (r 1 )=E 0 (r 2 )=̂E 0 , which becomes a real vector. Thus, system (1) can be recast into a single equation with a single unknown,

E ext = E x ext u x + E y ext u y : α ̲ =+ ̲ ̲ ̲ E Er r E (, ) ext 01 2 ext (2) α ̲ =-̲ ̲ - E rr E [1 ( , )] ext 12 1 0 (3) θ θ ̲ ̲ = -̲ + ̲ • - ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ E E s s E E
where

α πε ̲ =̲̂ - s kd d 1i 4 e kd 0 3 i (5)
Hence, calculating E ext just implies the inversion of a simple 2 × 2 diagonal matrix. It finally yields

θθ ̲ = -̲ ̲ = + ̲ EE s EE s cos 12 , sin 1
xy ext 0 ext 0 (6)
The condition of isosbesticity reads "|E ext | 2 independent of θ" where

θθ | ̲ |= -̲ + + ̲ E s E s E cos 12 sin 1 ext 2 0 2 0 2 (7)
This implies |1 -2s| = |1+s|, or equivalently

| ̲ -| = s 11 (8)
The expression of s can be expressed from eq 5 by calculating α using the following closed-form expression, valid for small values of ka: [START_REF] Maier | Fundamentals and Applications[END_REF] αε ε εε εε ε . Refined estimations of α can be obtained numerically using the Mie theory of higher orders. [START_REF] Metwally | Threshold for Photothermal Bubble Generation Using Plasmonic Nanoparticles[END_REF] Hence, writing s = s′ +i s″, the condition of isosbesticity considered in the (s′, s″) complex plane is elegantly described by a circle of radius 1 centered in (s′, s″) = (1, 0). The derivation of this condition is provided in the SI (as a Mathematica notebook). This condition has to be compared with the actual expression of s, which can be expressed from eq 5 using, for instance, a closedform expression of α 24 or Mie theory. [START_REF] Metwally | Threshold for Photothermal Bubble Generation Using Plasmonic Nanoparticles[END_REF] In Figure 2, representing the (s′, s″) complex plane, the condition of isosbesticity (eq 8) corresponds to a circle represented by a solid line, and the actual values of s defined by eq 5 are represented by dashed lines, for three different dimer geometries. The isosbestic conditions correspond to crossings of these two lines. In the case a = 30 nm and g = 10 nm, an isosbestic condition is fulfilled at λ 0 = 527 nm. Noteworthily, this is in good agreement with the numerical simulations reproduced in Figure 1 for the same morphology that predicted an isosbestic wavelength of 531 nm.

̲ = - ̲ + +- ̲ +- ε εε ̲ - ( ) () V x xx 1( ) ( 1
Making the spheres bigger and reducing the gap make the isosbestic wavelength larger, up to 550 nm in the case of a = 40 nm and g = 5 nm. The effect of enlarging an asymmetric structure (no rotational invariance along the direction of propagation) such as a sphere-dimer is to damp the stronger (longitudinal) resonance with respect to the weaker (transverse). As a consequence, the two spectra get closer, which increases the chance to observe crossing points. However, for smaller spheres and/or larger gaps, no isosbesticity condition occurs, as in the case of a = 20 nm and g = 10 nm. Thus, the existence of an isosbestic wavelength for a given nanoparticle geometry is not systematic.

Isosbesticity with Materials Other than Gold. As shown above, for gold, no solution of eq 8 exists for small enough nanoparticles, i.e., for small values of |s| (lim a→0 s = 0). One could wonder why and one could wonder if this happens for any material. Equation 8 can be fulfilled at the first order for arbitrarily small values of s if s is an imaginary number. For small nanoparticles, one can consider a first-order approximation of s using eqs 5 and 9, which yields

ε εε εε ̲ = ̲ - ̲ + s a d 2 s 3 3 s s (10)
Consequently, eq 8 can be fulfilled for arbitrarily small nanoparticles only if the permittivity of the material can fulfill

εε εε ℜ ̲ - ̲ + = ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 2 0 s s (11)
or equivalently |ε| 2 + ε′ = 2, or equivalently

εε ′+ + ″ = ⎜ ⎟ ⎜ ⎟ ⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎞ ⎠ 1 2 () 3 2 2 2 2 (12)
which is the equation of a circle in the complex plane of ε. This condition never happens for gold irrespectively of the wavelength. It is actually even very rare in plasmonics. According to our investigations, it only occurs for silver (Ag) and zirconium nitride (ZrN). Interestingly, this result makes Ag and ZrN singular materials in plasmonics. Figure 3 displays polar plots of the complex permittivity of several metals used in plasmonics, along with the circle corresponding to eq 12. The materials we have considered are Ag, Al, Au, Co, Cr, Cu, Fe, Mn, Mo, Ni, Pd, Pt, Ta, TiN, and ZrN, and we have only plotted the permittivity values down to λ = 300 nm. Other materials can exhibit isosbestic points such as Hg or Al, but at very low wavelengths close to 100 nm. They have not been represented in Figure 3.

These results are corroborated by numerical simulations presented in Figure 4a-f. Both small and large nanoparticles composed of Au, ZrN, and Ag have been modeled. While a small sphere-dimer made of gold does not exhibit isosbestic points (Figure 4a), small sphere-dimers made of ZrN and Ag feature isosbestic conditions for specific wavelengths (Figure 4b,c). For larger structures, isosbestic conditions can be obtained for any material (Figure 4d-f).

■ THE CASE OF SPHEROIDS

Spheroids represent another textbook case that deserves attention. For small spheroids in the quasistatic regime, the condition of isosbesticity can be derived analytically as well. Using Mie-Gans theory [START_REF] Myroshnychenko | Modeling the optical response of gold nanoparticles[END_REF] and the expression of the Joule number defined in ref 28, the isosbestic condition reads

ξξ | -| = | - | LL 13 13 11 3 3 (13)
where

ξ εε εε ε = - +-L 1 3( ) j j s ss ( 14
)
and where L j are the depolarization factors. For a sphere L x = L z = 1/3. For a prolate spheroid along the z direction, L z < L x , but 2L x + L z = 1/3 remains constant. Replacing ξ j in eq 13 by its expression in 14, one gets Let us now consider a small elongation from a sphere along the z axis. In this case, one has

εε εε ε εε εε ε - - +- =- - +- L L L L 1 () () 1 () ()
δ =+ L L 1 3 x (18) δ =- L L 1 3 2 z ( 19
)
where δL is arbitrarily small. Equation 17 now gives

ε ε εε δ ε ε εε δ ξδ ξδ ξδ ξδ ++- = +- - |+ | = |- | + - = LL LL L L 2 33
where

ξ εε εε = - + 2 s s ( 21
)
At the first order of δL, eq 20 gives

ξδ ξδ ξδ |+ + | ≈ |+ | ≈ LL L (1 3 )(1 6 ) 1 19 1 (22)
Just like eq 8, this equation has a solution for δL arbitrarily small if ξ is imaginary, i.e., if

ξ εε εε ℜ= ℜ ̲ - ̲ + = ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ () 2 0 s s ( 23
)
which is exactly the same condition as in the case of a spheredimer. Consequently, eq 23 seems to be a general condition to achieve isosbesticity in plasmonics for small nanoparticles. Since we could derive this condition only for dimers and elongated nanoparticles, this rule of thumb is only conjecture. These results are corroborated by numerical simulations presented in Figure 5a-f. Both small and large nanoparticles composed of Au, ZrN, and Ag have been modeled. While a small spheroid made of gold does not exhibit isosbestic points (Figure 5a), small spheroids made of ZrN and Ag feature isosbestic conditions for specific wavelengths (Figure 5b,c). For larger structures, isosbestic conditions can be obtained for any material (Figure 5d-f).

■ HOW TO ADJUST THE ISOSBESTIC WAVELENGTH

Just like the case of a sphere-dimer again, the isosbesticity wavelength for a spheroid remains in the green region of the spectrum, close to the resonance wavelength of a dipolar gold sphere. For structures with 2-fold symmetries, the isosbestic wavelength is actually supposed to be close to the resonance wavelength corresponding to the transverse mode, i.e., when the polarization of the incoming light is along the shorter dimension of the structure. This is why sphere-dimers and spheroids are isosbestic around λ = 530 nm. In order to shift the isosbestic wavelength to higher values, it suffices therefore to consider structures with a red-shifted transverse resonance. Structures made by e-beam lithography typically behave this way. We conducted numerical simulations of the cross-polarization spectra of disc-dimers of various thicknesses (see Figure 6). In each case, the discs are 60 nm in diameter and separated by a gap of 10 nm, but their thickness varies from 10 to 60 nm. Varying the thickness makes it possible to vary the transverse resonance. For a thickness of 60 nm, corresponding to the diameter of a disc, the resonance is close to that of a sphere and the isosbestic wavelength is 542 nm, close to the resonance of a sphere. For much smaller thicknesses, the transverse aspect ratio is increased, which naturally increases the transverse resonance wavelength and the isosbestic wavelength. Results of this section evidence that engineering the transverse resonance of plasmonic nanostructures enables the selection of isosbestic conditions over the entire visible-near-IR range using gold. Disc-dimers, which can be easily designed using e-beam lithography, appear to be efficient structures for this purpose.

The absorption spectra of a disc-dimer for longitudinal and transverse illuminations have been calculated using BEM (Figure 6). Both polarizations feature a red-shifted resonance, which naturally drives the isosbestic wavelengths to higher values. Thus, the isosbestic wavelength can be a priori adjusted at will, which is a strong benefit to address a wide range of applications.

■ DISCUSSION

The influence of the light polarization on plasmonic structures and related invariances have been reported in the literature. First, the presence of isosbestic points in plasmonics were first mentioned by Thompson et al. in 2011, 29 but the definition was different. This work rather dealt with the transmission intensity through plasmonic arrays, not photothermal effects. Second, in 2016, Geraci et al. investigated invariances with plasmonic trimers. [START_REF] Geraci | Polarisation-Independent Enhanced Scattering by Tailoring Asymmetric Plasmonic Systems[END_REF] The authors explained that the usual strong dependence of the optical response of plasmonic nanoparticles on the incident light polarization can be problematic for some applications. To address this problem, disc-trimer structures have been shown to be possible candidates. However, the authors focused on a different objective compared to our work: a broadband invariance of the optical response while we are seeking invariance of the heat generation while maintaining substantial dependence of the optical response. Finally, Lachaine et al. evidenced an invariance of heat generation from gold spheres from linear to circular polarizations. [START_REF] Lachaine | From Thermo-to Plasma-Mediated Ultrafast Laser-Induced Plasmonic Nanobubbles[END_REF] Isosbesticity is a concept that can be beneficial for several applications in plasmonics. In nanochemistry, plasmonic nanostructures can favor chemical reactions. [START_REF] Baffou | Nanoplasmonics for Chemistry[END_REF] However, the origin of this enhancement is not always straightforward. Both an NFE and an HG can be at the origin of an increase of the yield of a chemical reaction, and both these quantities are difficult to discriminate at the nanoscale in samples used in nanochemistry. In 2009, this problem was honestly pointed out by the group of El Sayed. [START_REF] Yen | Plasmonic Field Effect on the Hexacyanoferrate (III)-Thiosulfate Electron Transfer Catalytic Reaction on Gold Nanoparticles: Electromagnetic or Thermal?[END_REF] Very often, thermal effects are experimentally [START_REF] Christopher | Visible-Light-Enhanced Catalytic Oxidation Reactions on Plasmonic Silver Nanostructures[END_REF] or numerically [START_REF] Neumann | Solar Vapor Generation enabled by nanoparticles[END_REF] underestimated because of inappropriate approaches. A plasmonic structure behaving just like in Figure 1 could be used to straightforwardly discriminate one mechanism or another (HG or NFE), experimentally. If the chemical yield produced by nontrivial isosbestic nanostructures proves independent of incident light polarization, it will evidence a thermal-induced effect, while if it is strongly polarizationdependent, it will rather highlight a photochemical process. Other applications of plasmonics also suffer from this blurred duality, such as in HAMR: [START_REF] Stipe | Magnetic Recording at 1.5 Pb m -2 Using an Integrated Plasmonic Antenna[END_REF][START_REF] Pan | Heat-Assisted Magnetic Recording[END_REF] in this application, the near-field intensity of a plasmonic nanoantenna localized a few nanometers away from a hard-disk drive is supposed to transfer energy via its optical near-field toward the disk to locally increase the disk temperature. However, a temperature increase due to thermal radiation (not optical radiation) is also possible, although it was never demonstrated. The use of isosbestic plasmonic structures could decipher between both transduction pathways: an optical near-field transfer or a thermal radiation near-field transfer. Finally, the origin of bubble generation around plasmonic nanoparticles under pulsed illumination is also unclear and suffers from the same dual issue: besides a simple heating mechanism, the generation of a plasma due to a strong electric near-field enhancement could also be responsible for bubble generation, especially under femtosecond pulsed illumination. [START_REF] Lachaine | From Thermo-to Plasma-Mediated Ultrafast Laser-Induced Plasmonic Nanobubbles[END_REF][START_REF] Boulais | Plasma Mediated off-Resonance Plasmonic Enhanced Ultrafast Laser-Induced Nanocavitation[END_REF] Generating a bubble around metal nanoparticles under pulsed illumination at the isosbestic wavelength of the nanoparticles (typically rods immobilized on a substrate) could confidently reveal the nature of the underlying mechanism of bubble formation.

Note that the concept of photothermal isosbesticity in plasmonics could be extended to the case of circular polarization and chiral nanostructures. [START_REF] Fan | Chiral Nanocrystals: Plasmonic Spectra and Circular Dichroism[END_REF][START_REF] Valev | Chirality and Chiroptical Effects in Plasmonic Nanostructures: Fundamentals, Recent Progress, and Outlook[END_REF][START_REF] Schaferling | Chiral Nanophotonics[END_REF] The optical properties of chiral plasmonic structures depend on the left-or right-handed nature of an incoming circular illumination. Thus, possible isosbesticity could be a priori obtained when observing crossing points in the absorption spectra corresponding to left-and right-handed polarization, i.e., for zero circular dichroism (CD), as defined in ref 36. This idea will possibly be further investigated in a future study.

In summary, we introduced the concept of isosbestic plasmonic structures in nanoplasmonics and give all the recipes to make an efficient use of it: (i) Isosbestic conditions are not always possible for a given nanoparticle morphology, (ii) Ag and ZrN are singular isosbestic materials, (iii) the isosbestic wavelength can be found for a given morphology by considering the intercrossing of two cross-polarized absorption spectra, (iv) enlarging nanoparticles favor the occurrence of isosbesticity, (v) the isosbestic wavelength can be arbitrarily chosen by engineering the transverse resonance wavelength of the nanoparticle. Finally, we have explained how such a class of plasmonic structures can have powerful applications for elucidating mechanisms in many applications such as nanochemistry, heatassisted magnetic recording, or microbubble formation.

Figure 1 .

 1 Figure 1. (a) Absorption cross section of a dimer structure composed of two 60 nm gold spheres in water separated by a 10 nm gap, for longitudinal and perpendicular polarizations of the incident light. Gray lines correspond to intermediate polarization angles (15°, 30°, 45°, 60°, and 75°). A crossing is observed for the isosbestic wavelengths λ = 496 and 531 nm. (b) Absorption cross section and normalized electric near-field as functions of the illumination angle, at the isosbestic wavelength λ = 531 nm, identified in (a). (c) Temperature map for a longitudinal polarization of the incident light and an arbitrary light irradiance of 1 mW/μm 2 . (d) Corresponding normalized electric near-field. (e, f) Same as (c) and (d) for a perpendicular incident light polarization.

ε s = n s 2 ,

 2 ε is the permittivity of the nanoparticle material (taken from Johnson and Christy's[START_REF] Johnson | Optical Constants of the Noble Metals[END_REF] data set in our simulations using gold), =

Figure 2 .

 2 Figure 2. Space of s values for a gold sphere-dimer in water in the complex space. The black solid line represents the values that s must equal to obtain an isosbestic sphere-dimer. Dashed lines represent the actual values of s as a function of the wavelength from λ = 400 to λ = 1000 nm, for three different gold sphere-dimer morphologies (a, g). Particular wavelength values have been indicated, in particular the isosbestic wavelength (boxed numbers).

Figure 3 .

 3 Figure 3. Polar plots of the complex permittivity of common plasmonic materials. The condition given by eq 11 is represented as a solid line circle. Only materials crossing this circle can exhibit isosbestic conditions in the quasistatic regime.

Figure 4 .

 4 Figure 4. Calculations of cross-polarization absorption spectra using the BEM for different spheroid sizes and materials. (a-c) Small sphere-dimer composed of two 20 nm spheres separated by a gap of 2 nm. (d-f) Large sphere-dimer composed of two 200 nm spheres separated by a gap of 20 nm. (a, d) Case of Au. (b, e) Case of ZrN. (c, f) Case of Ag. Dashed lines correspond to calculations with light polarization parallel to the axis of the dimer. Solid lines correspond to light polarization perpendicular to the dimer axis. Dotted lines correspond to an intermediate light polarization (45°).

Figure 5 .

 5 Figure 5. Calculations of cross-polarization absorption spectra using the BEM for different spheroid sizes and materials. (a-c) Small prolate spheroids of 10 nm × 20 nm. (d-f) Large prolate spheroids of 100 nm × 200 nm. (a, d) Case of Au. (b, e) Case of ZrN. (c, f) Case of Ag. Dashed lines correspond to calculations with a light polarization parallel to the long axis of the spheroid. Solid lines correspond to a light polarization along the short axis of the spheroid. Dotted lines correspond to an intermediate light polarization (45°).

Figure 6 .

 6 Figure 6. Calculations of cross-polarization absorption spectra using the BEM for different gold disc-dimer thicknesses e: (a) e = 60 nm, (b) e = 30 nm, and (c) e = 10 nm. Insets represent the morphology of the structures. The illumination axis is vertical on these insets. Dashed lines correspond to calculations with light polarization parallel to the axis of the dimer. Solid lines correspond to light polarization perpendicular to the dimer axis. Dotted lines correspond to an intermediate light polarization (45°).

■ AUTHOR INFORMATION

Corresponding Author *E-mail: guillaume.baffou@fresnel.fr.

ORCID

Guillaume Baffou: 0000-0003-0488-1362

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

The authors wish to acknowledge financial support from ITMO Cancer AVIESAN (Alliance Nationale pour les Sciences de la Vie et de la Sante, National Alliance for Life Sciences & Health) within the framework of the Cancer Plan (GRAVITY project) and from the Agence National de la Recherche (PlasmoTherm project, ANR-15-CE24-0025-01).


G Prabhudesai

S Perrard

F Pétrélis

S Fauve

Experimental observation of spontaneous temperature fluctuations in turbulent flows

Even in the absence of externally applied temperature gradients, spontaneously generated temperature fluctuations arise in turbulent flows. We experimentally study these fluctuations in a closed von Kármán swirling flow of air at Mach number of order 10 -3 , whose boundaries are maintained at a constant temperature. We observe intermittent peaks of low temperature correlated with pressure drops within the flow and show that they are caused by vorticity filaments. The measured ratio of temperature to pressure fluctuation agrees with the prediction based on adiabatic cooling within vortex cores. This experimental study shows that although the Mach number of the flow is small, there exist regions within the flow where compressible effects cannot be discarded in the equation for temperature and locally dominate the effect of viscous dissipation.

Introduction. -Temperature fluctuations spontaneously generated in turbulent flows have not been studied experimentally so far, in contrast to advection of temperature in the presence of an externally applied temperature gradient [1,[START_REF] Obukhov | Tech. Rep. Air force systems command Wright-Patterson[END_REF]. Two processes are responsible for spontaneously generated temperature fluctuations. First, viscous dissipation of the flow kinetic energy generates heat, as recently studied in detail using numerical simulations [START_REF] Marinis | [END_REF][4][5]. Fluctuations of temperature result from bursts of dissipation that occur at small scales and in a strongly intermittent way in space and time. These fluctuations could therefore provide additional information about small scale intermittency of turbulent flows. There exists a second mechanism other than dissipative processes, which subsists in the limit of a perfect fluid with conserved entropy. It is related to heating or cooling induced by pressure fluctuations generated by the flow. The equation for temperature in the absence of any externally applied large scale gradients is given by [START_REF] Monin | Statistical fluid mechanics. 1 1st Edition[END_REF],

ρc p DT Dt = ρν 2 Σ 2 + σ Dp Dt + α △T (1)
where D /Dt = ( ∂ /∂t + u • ∇) and Σ = (∇u) + (∇u) T respectively denote the material derivative and strain rate tensor, and, Σ2 = Σ ij Σ ji . The quantities of density, specific heat capacity at constant pressure, kinematic viscos-ity, thermal conductivity and coefficient of thermal expansion are denoted by ρ, c p , ν, α and β respectively. We denote the non-dimensionalized coefficient of thermal expansion by σ = βT . The first two terms on the right hand side of eqn. 1 correspond to two sources of spontaneously generated temperature fluctuations; viscous dissipation and compressible heating/cooling driven by pressure fluctuations. As demonstrated by Bayly et al. [START_REF] Bayly | [END_REF], both the source terms in eqn. 1 are retained in the limit of incompressible flow with the quantities replaced by ρ * , c * p , ν * , α * , β * and σ * = β * T * . The superscripted asterisk denotes the values of the corresponding quantities measured in the absence of flow. We define the Reynolds number as Re = U L/ν * and Peclet number as P e = ρ * c * p U L/α * where U , L are the characteristic velocity and length scales of the flow. In the limit of incompressible flow, eqn. 1 then implies that when σ * Re ≫ 1, temperature fluctuations are mostly driven by pressure fluctuations 1 . Additionally, in this case, the fluid will undergo adiabatic process if the condition of P e = P rRe ≫ 1 is also met, where the Prandtl number is defined as P r = ρ * c * p ν * /α * . When σ * Re ≪ 1, they are mostly generated by viscous dissipation. The above analysis, though, is not as straightforward for turbulent flows where a wide range of scales coexist and interact, with the flow exhibiting a range of Reynolds numbers across scales. We can nonetheless ascertain that the temperature fluctuations at small scales (Re ≪ 1) would be governed by viscous dissipation and at large scales (Re ≫ 1) by pressure fluctuations. We note that this second mechanism has been discarded in numerical simulations of temperature equation, even though it should certainly be taken into account at large Reynolds numbers, for instance in astrophysical or geophysical flows. At the laboratory scale, adiabatic cooling is observed within strong coherent vortices such as in the Ranque tube [8] or in von Kármán swirling flows between co-rotating disks [9,[START_REF] Morales | Ecoulements tourbillonnaires de von Kármán Ph[END_REF] where a significant temperature drop can be observed in the vortex core [START_REF] Labbé | [END_REF]. Generating a von Kármán turbulent swirling flow in air (σ * ≈ 1, Pr = 0.71) between counter-rotating disks, we show that even in the absence of an externally driven vortex and at moderate values of the Reynolds number, these compressible events contribute significantly to the temperature fluctuations that are generated spontaneously by the flow.

Experimental set-up and results. -We generate a von Kármán swirling flow in air between two counterrotating disks in a closed cylinder of diameter 130 mm made out of copper of 2 mm thickness (see fig. 1

top).

Copper has been chosen for its high thermal conductivity to minimize temperature gradients along the boundaries of the experiment. The top and the bottom faces of the cylinder are closed with circular copper plates of the same thickness as the cylinder. In order to maintain fixed the temperature at the boundaries, a copper tubing of outer diameter 10 mm is welded to the cylindrical part of the container. The inlet and the outlet of the copper tubing are connected to a circulating water bath with thermostat control which maintains the temperature of the water at a given fixed value within ±0.01 K.

We use two loop-controlled brushless DC motors with the angular velocity Ω ranging from of 0 -2000 rotations per minute (rpm). Each of the motor drives a disk with four curved blades which generate a strong flow inside the cylindrical cavity. The thickness of the disks and the height of the blades are 7.5 mm. Holes drilled on the surface of the cylinder provide access to the probes: one cold-wire temperature probe (Dantec 55P31)2 , one 1D hot-wire velocity probe (Dantec 55P16) and one acceleration compensated piezoelectric pressure probe (PCB 103B02) 3 . All probes are placed in the midplane, with the pressure probe flushed to the wall as shown in the top-view sketch of fig. 1 bottom. The mean power consumption per unit mass by the turbulent flow ǫ is measured from the power required by the motors to maintain the flow, ranging from 0 to 500 m 2 /s 3 . This gives the Kolmogorov length scale η = ν 3 / ǫ 1/4 ≥ 50 µm and the corresponding Kolmogorov timescale τ η = ǫ -1/3 η 2/3 ≥ 0.2 ms with ν = 1.5 × 10 -5 m 2 /s the kinematic viscosity of air at room temperature. We evaluate the Taylor microscale λ = 15νu 2 rms / ǫ 1/2 ≥ 0.8 mm, the corresponding Taylor timescale τ λ = ǫ -1/3 λ 2/3 ≥ 1 ms and the Taylor microscale based Reynolds number at the highest rotation rate Re λ ∼ O(10 2 ) using standard estimates for homogeneous and isotropic turbulent flows [START_REF] Tennekes | A first course in turbulence nachdruck des originals von 1972 Edition[END_REF]. The cold-wire temperature probe should thus be able to resolve turbulent scales as small as λ. Data is acquired from the probes at a sampling frequency of 5 kHz using a NI DAQ system. Two point spatial correlation of velocity is evaluated in a different set of experiments using two 1D hot-wire probes moved radially about the center axis in the midplane. An exponential decay is observed giving an integral length scale of l 0 = 5 mm. The Mach number defined as M a = u rms /c * reaches a maximum value of 4 × 10 -3 for the current setup where u rms is the RMS of velocity fluctuations and c * is the isentropic speed of sound in the fluid at rest. The intermittent pressure drops in a turbulent flow have been fairly well-studied and so is the physical phenomenon associated with it. These have been attributed to vorticity filaments [START_REF] Douady | [END_REF][14][15][16][17] and were shown to result in the behaviour at low frequency of the pressure power spectrum [18]. The observed amplitude of pressure drops is of O(ρ * R 2 Ω 2 ) which is in agreement with previous experimental observations [15]. Similar to the pressure signal, we observe intermittent temperature drops much larger than the standard deviation of the temperature fluctuations. They correspond to the low exponential tail of the temperature PDFs (fig. 3a). Previous studies [19][20][21] have shown that the process of viscous dissipation is intermittent, with the PDFs of ǫ being positively skewed. As observed from eqn. 1, this would result in peaks in temperature and the process of viscous dissipation thus cannot explain the observed drops.

One immediate question that arises is whether the drops observed in pressure and temperature signals are correlated and result from the same type of structures in the turbulent flow. To check that, we evaluate the joint PDFs between pressure and temperature signals. To do so, we place the temperature probe close to the pressure probe (d ≈ 6 mm apart) near the boundary. The joint PDF between pressure and temperature is defined as,

Π(p 0 , T 0 , ∆t)dp 0 dT 0 = (2) prob p(t) ∈ [p 0 , p 0 + dp 0 ], T (t + ∆t) ∈ [T 0 , T 0 + dT 0 ]
for a time lag ∆t between the pressure and temperature signals. The correlation between pressure and temperature shows a maximum for a non zero time lag ∆t max ∼ 2 ms which can be attributed to the advection of fluctuations between the two probes by a mean flow. This gives a characteristic advection velocity of u c = d/∆t max = 3 m/s which is consistent with the velocity associated with disk rotation of O(RΩ). Fig. 3b shows the joint PDF Π(p 0 , T 0 , ∆t = ∆t max ) obtained for a rotation rate of Ω = 2000 rpm. We observe that the joint PDF is skewed towards negative values of both pressure and temperature fluctuations. We conclude that the vorticity filaments are responsible not only for the distribution of pressure drops but also for the temperature drops4 . From eqn. 1, in the limit of σ * Re ≫ 1 and P e ≫ 1, and for the fluid being an ideal gas, the amplitudes of pressure (∆p) and temperature (∆T ) satisfy the adiabatic relation,

∆T ∆p = γ -1 γ T * p * ≈ 0.84 mK/Pa, (3)
where γ ≈ 1.4 for air at normal temperature and pressure. For a given value of p, we can define the temperatures T max (p) and T avg (p). The value of temperature that maximizes the conditional PDF Π(T |p, ∆t max ) is denoted by T max (p), whereas, T avg (p) = +∞ -∞ T Π(T |p, ∆t max )dT denotes the conditional average of T on p. Fig. 4a shows the values of these two temperatures extracted from fig. 3b. We observe that for negative values of pressure fluctuations, T max and T avg depend linearly on p with a slope of 0.84 mK/Pa as predicted by eqn. 3. From this observation we conclude that concurrent negative fluctuations in pressure and temperature are associated with adiabatic cooling. It is worth noting that the temporal width of the drops fluctuates following an exponential distribution for large values, with a characteristic time of 5 ms for pressure drops and 1 ms for temperature drops. The difference results from the finite size of pressure probe, as studied by Abry et al. [18] using different probe sizes. In our experiments, the dimensions of the pressure and temperature probes are 2.5 mm and 1 µm respectively. Since the physical dimension of the measuring element of the temperature probe is smaller than any length scale associated with the turbulent flow, it turns out to be a better probe to access the internal structure of vorticity filaments.

From previous studies [15,22], we can estimate the mean value of σ * Re f and P e f in the core of vorticity filaments. The typical velocity scale of the vorticity filaments is u rms and their width lies between η and λ. Thus for the current experimental setup, we obtain 5 < σ * Re f < 100 and 3.6 < P e f < 71 corresponding to σ * Re f > 1 and P e f > 1, for which the adiabatic cooling/heating dominates over viscous dissipation, as we observe experimentally.

Conclusion. -

We have provided experimental evidence that the vorticity filaments in turbulent flows contribute to temperature fluctuations even for small Mach numbers of order 10 -3 through adiabatic effects. This is evidenced by the correlation between the drops observed in both pressure and temperature fluctuations. It should be noted that our observations are not limited to gas. If one considers the case of water at 20 • for which σ * ≈ 0.06 and P r ≈ 7, temperature fluctuations would be driven by adiabatic heating/cooling when Re ≫ 16, a value easily attained in water flows. The magnitude of temperature fluctuations thus generated would be 70 times smaller than those in air for same velocity as deduced from eq. 1 but should be measurable at larger velocities since the pressure drops scale as the square of velocity.

Vorticity filaments being characteristic of all turbulent flows, this process should be considered in a variety of contexts. Not only when studying temperature fluctuations, but also when studying a field that is coupled to temperature, possible examples of which are quite numerous. In most cases, the field under consideration, say n give an example, will display relative fluctuations equal to the ones of temperature, so that dn/n = dT /T . We measure here temperature drops of relative amplitude of O(10 -2 ) percent. As they are quadratic in velocity, large relative fluctuations can easily be observed with a moderate increase in velocity. One example is the scattering of acoustic or electromagnetic waves due to temperature fluctuations in turbulent flows [START_REF] Tatarskii | Wave propagation in a turbulent medium dover edition Edition Dover books on physics[END_REF], in which case n denotes the refractive index of the medium.

Another example relevant at the astrophysical scale is that of interstellar medium which is known to be turbulent [START_REF] Larson | [END_REF]25] and displays intermittency [26]. Later studies (see for example [27][28][29]) showed the effects of intermittent turbulent structures, in particular dissipative structures on the chemistry of interstellar medium. Our work shows that pressure-driven adiabatic temperature fluctuations can dominate over the effect of viscous dissipation and thus calls for studying the effect of these fluctuations on the chemistry of interstellar medium.

In the domain of cloud physics, an open question is how rain forms in warm cumulus on a short time scale (less than one hour) [START_REF] Rogers | A short course in cloud physics 3rd Edition International series in natural philosophy[END_REF]. Recent investigations have highlighted the effect of turbulent fluctuations on the rapid growth of cloud droplets across the size gap from 10 to 50 µm. Several mechanisms have been considered: enhanced geometric collision rates by air turbulence [START_REF] Falkovich | [END_REF], effect of turbulence on collision efficiencies, effect of stochastic fluctuations [32,33]. The present study shows that condensation induced by adiabatic cooling which is essentially the working principle of cloud chambers, should also be taken into account. How adiabatic cooling observed in vorticity filaments affects the formation of cloud droplets and whether this contributes to the efficiency of precipitation in warm rain process deserves further studies. * * *

  Fig. 2a and 2b show the time series of temperature and pressure fluctuations about their mean values, obtained for Ω = 2000 rpm. Sharp drops are observed in both signals, examples of which are shown in the respective inset figures.

Fig. 1 :

 1 Fig. 1: Sketch of the von Kármán flow configuration driven by two counter-rotating disks. Temperature of the boundaries is maintained constant using a circulating water bath. The pressure, velocity and temperature probes are all placed in the midplane.

Fig. 2 :

 2 Fig. 2: (a) Time series of temperature fluctuations measured in the bulk. (b) Time series of pressure fluctuations measured at the boundary. Inset figures show the zoom of one drop observed in temperature and pressure fluctuations. The drops are intermittent, i.e., they are observed to occur randomly in time.

Fig. 3 :

 3 Fig. 3: (a) PDFs of temperature fluctuations normalized by their respective RMS values for Ω = 1200 rpm (red), 1600 rpm (blue) and 2000 rpm (green). (b) The joint PDF of pressure and temperature fluctuations when the two probes are 6 mm apart for Ω = 2000 rpm. The isocontours of probability are plotted on a natural logarithmic scale. The black dashed line has a slope of 0.84 mK/P a.

Fig. 4 :

 4 Fig. 4: Tmax (the temperature at which the joint PDF between p and T is maximized, ) and Tavg (the conditional average of T on p, △) vs p for Ω = 2000 rpm. The black dashed line shows the prediction for adiabatic process (eqn. 3).

This condition is obtained on comparing the magnitudes of the first two terms on the right hand side of eqn. 1 which gives (σ * Dp/Dt) (ρ * ν * Σ

) ∼ (σ * ρ * U

/L) (ρ * ν * U 2 /L 2 ) = σ * Re, using p ∼ ρ * U 2 and Σ ∼ U/L.

Diameter of the filament dc = 1 µm, response time τ c r = 0.5 ms, resolution of 0.2 mK, uncertainty of ±1 mK.

Diaphragm diameter dp = 2.5 mm, frequency range at ±3 dB is 1 Hz -7 kHz, resolution of 0.15 Pa, uncertainty of ±1 Pa.

If the pressure and temperature signals were independent, the isocontours of joint PDFs would be skewed along negative p-axis and negative T -axis. This would be due to the PDFs of p and T individually being negatively skewed but without correlation between them.

This work has been supported by the Agence nationale de la recherche (Grant No. ANR-17-CE30-0004), CE-FIPRA (Project 6104-1) and by CNES (action 6291).


Ana P Estrada-Vargas

Member IEEE E López-Mellado

Jean-Jacques Lesage Member
email: jean-jacques.lesage@lurpa.ens-cachan.fr

Off-line Identification of Concurrent Discrete Event Systems Exhibiting Cyclic Behaviour

Keywords: Discrete event systems, Automated model identification, Interpreted Petri nets I

This paper presents a method for the identification of concurrent Discrete Event Systems (DES) from input-output sequences representing the observed behavior. The proposed off-line technique yields an input-output model expressed as an Interpreted Petri net (IPN), which represents exactly the language than that generated by the observed system, which may include cyclic sequences. First, a sample of inputoutput vectors words are processed for obtaining sequences of output changes called events; then sequences of -length event traces are built and represented by an IPN model composed by non measurable places. Then measurable places are added; they are related to transitions representing pertinent output changes. Finally inputs are associated to transitions and implicit non measurable places are removed.

INTRODUCTION

Automated building of discrete event models from external observation of system behavior has been addressed as a problem of system identification. The problem has interesting applications such as reverse engineering for (partially) unknown systems, fault diagnosis or system verification.

Previous results on the matter have been proposed as a solution to the problem of grammatical inference. In [START_REF] Gold | Language Identification in the Limit[END_REF] a finite automaton (FA) is built from positive samples of accepted words. Later several methods for obtaining Mealy [START_REF] Kella | Sequential Machine Identification[END_REF] [16] and Moore [START_REF] Biermann | On the Sysnthesis of Finite-State Machines from Samples of Their Behavior[END_REF] [START_REF] Veelenturf | An Automata theoretical approach to developing learning neural networks[END_REF] machines have been proposed. Also building of context free grammars has been studied [START_REF] Levy | Skeletal structural descriptions[END_REF], [START_REF] Takada | Grammatical Inference for Even Linear Languages Based on Control Sets[END_REF], [START_REF] Ishizaka | Polynomial Time Learnability of Simple Deterministic Languages[END_REF].

The identification of Petri net (PN) models has been proposed for coping with more complex systems exhibiting concurrent behavior. In [START_REF] Hiraishi | Construction of Safe Petri Nets by Presenting Firing Sequences[END_REF] an algorithm for constructing Petri net models is presented. First, the language of the target system is identified in the form of deterministic FA (DFA). Then, the algorithm obtains from the DFA the structure of a PN that accepts the obtained language. Three different approaches have been adopted in recent publications addressing the problem of identification of DES. The incremental synthesis approach proposed by Meda et al., deals with unknown partially measurable DES exhibiting cyclic behavior; in [START_REF] Meda | DES identification using Interpreted Petri nets[END_REF] an identification method based on the least square estimator is presented. Then in [START_REF] Meda-Campaña | On-line Identification of Discrete Event Systems: Fundamentals and Algorithms for the Synthesis of Petri Net Models[END_REF] [START_REF] Meda-Campaña | Identification of Concurrent Discrete Event Systems using Petri Nets[END_REF] several algorithms have been proposed allowing the on-line identification of concurrent DES. Although the techniques are efficient, the obtained models may represent more sequences than those observed.

Other technique oriented to fault diagnosis is presented in [START_REF] Klein | Identification of Discrete Event Systems for Fault Detection Purposes[END_REF]; it obtains non deterministic FA representing the same observed behavior, which is represented as a set of input/output sequences. The algorithms operate off-line efficiently.

The off-line techniques based on integer linear programming (ILP) approach lead to free-labeled PN models representing exactly the observed behavior [START_REF] Giua | Identification of free-labeled Petri nets via integer programming[END_REF]. However both the ILP problem statement from samples and their processing have exponential complexity. This approach is being applied to other PN classes in [START_REF] Cabasino | Identification of deterministic Petri nets[END_REF] [START_REF] Fanti | Fault diagnosis and identification of discrete event systems using Petri nets[END_REF].

In this paper a method for the identification of Petri net models from observed input-output sequences representing behavior of concurrent DES is presented. It gathers and extends some strategies form [START_REF] Meda-Campaña | On-line Identification of Discrete Event Systems: Fundamentals and Algorithms for the Synthesis of Petri Net Models[END_REF] and [START_REF] Klein | Identification of Discrete Event Systems for Fault Detection Purposes[END_REF] and focuses on offline case. Because of the inherent capacity of Interpreted Petri Nets (IPN) models to represent inputs and outputs of DES systems, internal states, non-controllable transitions and parallelism, the proposed method yields an IPN model which represents exactly the same language of length +1 than that generated by the system without taking into account information a-priori about the system other than its input and output signals. In the first stage, a sample of output vectors words are processed for obtaining sequences of output changes called events; then sequences of -length event traces are considered allowing building an underlying IPN of non measurable places. The IPN orders the events by relating the measurable places according to pertinent output changes. Finally, the inputs are added and the model is simplified by eliminating implicit non-measurable places.

The remainder of this paper is organized as follows. Section 2 overviews basic definitions on PN and IPN. Section 3 presents the proposed identification technique. Finally concluding remarks are given.

II. PETRI NET MODELS

This section presents the basic concepts and notation of PN and IPN used in this paper. ; and are the pre-incidence and post-incidence matrices respectively.

A marking function M: P Z represents the number of tokens (depicted as dots) residing inside each place. The marking of a PN with n places is usually expressed as an nentry vector. Z + is the set of nonnegative integers. Definition 2: A Petri Net system or Petri Net (PN) is the pair N = (G,M 0 ), where G is a PN structure and M 0 is an initial marking. In a PN system, a transition t j is enabled at marking M k if p i P, M k (p i ) I(p i ,t j ); an enabled transition t j can be fired reaching a new marking M k+1 which can be computed as M k+1 = M k + Cv k , where v k (i)=0, i j, v k (j)=1, this equation is called the PN state equation. The reachability set of a PN is the set of all possible reachable marking from M 0 firing only enabled transitions; this set is denoted by R(G,M 0 ). This work uses IPN, an extension to PN that allows associating input and output signals to PN models.

Definition 3: An IPN (Q, M 0 ) is a net structure Q = (G, , , , ) with an initial marking M 0 .
G is a PN structure = { 1 , 2 , ..., r } is the alphabet of input symbols i . = { 1 , 2 ,..., q } is the alphabet of output symbols i . : T { } is a labelling function of transitions with the following constraint:

t j ,t k T, j k, if p i I(p i ,t j ) = I(p i ,t k )
0 and both (t j ) , (t k ) , then (t j ) (t k ); represents a system internal event externally uncontrollable : R(Q,M 0 ) ( Z + ) q is an output function, that associates to each marking in R(Q,M 0 ) a q-entry output vector; q is the number of outputs. is represented by a q×n matrix, in which if the output symbol i is present (turned on) every time that M(p j ) 1, then (i,j)=1, otherwise (i,j)=0.

If a transition t j T of an IPN is enabled at marking M k , there are two possibilities: a) If (t j ) = i is provided to the system then t j can be fired. b) If (t j ) = and t j is enabled then t j must be fired. When an enabled transition t j is fired in a marking M k , then a new marking M k+1 is reached. This behaviour is represented as

1 k t k M M j
; M k+1 can be computed using the state equation:

M k+1 = M k + Cv k y k = M k
where C and v k are defined as in PN and y k (Z + ) q is the k-th output vector (or marking projection) of the IPN.

According to functions and , transitions and places of an IPN (Q,M 0 ) can be classified as follows.

Definition 4: If (t i )
the transition t i is said to be controllable. Otherwise it is uncontrollable. A place p i P is said to be measurable if the i-th column vector of is not null, i.e. ( ,i) 0. Otherwise it is non-measurable. P = P m P u where P m is the set of measurable places and P u is the set of non-measurable places.

Definition 5: A firing transition sequence of an IPN (Q, M 0 ) is a sequence = t i t j ...t k such that
, where is the length of the sequence . The set of all firing sequences is called the language of (Q, M 0 ) and it is denoted as .

Definition 6:

The input language of an IPN (Q, M 0 ) is defined as:

The output language of an IPN (Q, M 0 ) is defined as: Definition 7: The output language of length l of an IPN (Q, M 0 ) is defined as: .

III. IDENTIFICATION METHOD

A. Problem statement

An identification procedure of a DES builds a mathematical model that represents its behaviour from the measured evolution of the DES' inputs and outputs [START_REF] Meda | DES identification using Interpreted Petri nets[END_REF][10] [START_REF] Fanti | Fault diagnosis and identification of discrete event systems using Petri nets[END_REF]. Below the addressed identification problem is stated.

Consider a DES S with binary outputs revealing part of its internal states in which the outputs may evolve concurrently. The input data to the identification procedure is a set of output words describing behaviour and the sequences of input symbols in that yield the output sequences. Definition 8: The set of observed output words of a DES S is such that , where is the j-th observed vector in sequence and is the length of the output word .

Definition 9:

The observed output language of a DES S is defined as .

Definition 10: The output language of length l of a DES S is defined as .

Given a set of output words , the input signals of a DES S, and an accuracy parameter , the aim of the identification process is to obtain a safe IPN model (Q, M 0 ) such that . The parameter is used to adjust the accuracy of the identified model, similarly as proposed in [START_REF] Klein | Fault detection of Discrete Event Systems using an identification approach[END_REF]. It is assumed that for every output vector change an input signal i is recorded; otherwise is taken.

B. General approach

The method for identification consists of several stages that build systematically a safe IPN which represents exactly the sampled output language of length +1 of the DES.

From the output vector words, the event sequences are computed and then sequences of event substrings of length are built. Then every substring is associated to a transition of a PN, which describes the causal relationship between event substrings through paths including non-measurable places. Finally, the output changes provoked by events are described by marking changes in measurable places and then related to pertinent transitions in the PN. Notice that the number of nonobservable places is not predefined.

C. Sample processing 1) Event sequences

As stated before, the data obtained from the system to be identified is a set of sequences of output vectors , such that

, where refers to the jth observed vector in sequence . Sequences may have different length. From these sequences, strings of observed events are first computed. Definition 11. An observed event is the variation between two consecutive output vectors , ; it is computed as .

Then for every sequence , a sequence of observed events is obtained. During the process, if the difference has not been observed before, a new event is created ( .

The following example illustrates that the method copes also with DES exhibiting cyclic behaviour. The sequences i of the detected events e j associated to output changes are obtained:

Additionally consider that the input signals measured at the output changes are: a, , b, c, d, e, for e1, e2, e3, e4, e5, e6, respectively.

2) Sequences of -length event traces

In order to introduce the accuracy parameter for identification, from every sequence of events we compute sequences of event traces such that every is a substring of length of that finishes with . For the first elements of the trace sequence the event is used. Following with the previous example, the sequences of traces using are: for for

D. Building the basic structure 1) Representing event traces

Once the sequences of event traces have been obtained, every trace is related to a transition in the IPN through a function :T

; the firing of a transition implies that consecutive events related to such a transition have been observed. In order to preserve firing order between transitions, dependencies are created between them and associated with an observed marking through the function , which relates every non-measurable place with an observed marking, such that every transition has only one input place and one output place ( ); when an event trace is found again in a the associated dependency must be used if it leads to the same observed marking. Let e j be the last event in the trace ; the associated transition will be denoted as (more than one transition may have associated the same e j ). This strategy can be systematically performed following the next procedure. . //Create an initial empty set of transitions, an initial empty event trace set, and an initial set of places with only a marked place associated with the first observed vector.

Step 2. : 2.1 //Take as current. Since search operation has linear complexity and the algorithm implies the addition of a transition for every computed trace that has not been yet observed, then Algorithm1 is executed in polynomial time on the number and maximum length of observed output sequences.

. //

Property 1. The IPN G built through algorithm 1 represents all and only all the sequences (sub-sequences) in

Proof. By construction, every sequence is represented in G by a circuit starting from including the sequence of ; this circuit represents also the sequence i of events. Furthermore, the reuse of computed transitions having associated the same event traces, during the processing of subsequent , is done only when common paths of length are built, which does not introduce other sequences.

Using the previous algorithm the obtained IPN corresponding to the two sequences of event traces of the example is showed in figure 1. transitions with the same associated event, then merge and and their output places accordingly. This test is performed iteratively on the new obtained place, until there are no many output transitions of such a place sharing the same event. Rule 2:

//If has more than one input transition with the same associated event, then merge and and their respective input places accordingly.

Since merging has linear order and rules cannot be applied more times than the number of places in the net, Algorithm 2 is executed in polynomial time on the number of repeated transitions in the same place.

Property 2. G' preserves the sequences in G.

Proof. Let be the set of observable sequences from place . Consider . Before applying rule 1, , with , ; , . After applying rule 1, . Thus , and then language of G is not changed by the application of rule 1. A similar reasoning can be done for rule 2.

In the example the application of the rules lead to the model of figure 2. Since the initial place has two input transitions associated to e 6 , they can merge and their input places.

3) Concurrent transitions

Other transformations may be performed when there are transitions that appear in the sequences in different order In figure 2 notice that between the transition associated to e 3 and the new transition associated to e 6 there are paths with all possible permutations of e 4 and e 5 ; then, we can transform this into a concurrent component and we obtain the net showed in figure 3. Notice that this model preserves the same event sequences of the previous one. The simplification by analysis of concurrency is not strictly necessary for representing the event sequences; however the equivalent model with concurrent transitions may be simpler. Although the analysis could be inefficient when the number of paths in the subnet is large, usually this number is reduced.

E. Adding outputs and inputs 1) Representing outputs changes

Once the event sequences are represented in the basic model, it must be completed by adding the output changes represented by the events and their respective inputs. Recall that events are vectors computed from the difference of consecutive output vectors; thus e j relates measurable places representing the outputs yielding the incidence matrix corresponding to measurable places. The inputs are straightforward included in the model. This procedure is detailed below.

Algorithm 3. Representing outputs changes

Input: G', {e i (j)} Output; Q : the IPN including measurable places and inputs Step 1.

. Create q measurable places for every one of the components in the output vectors.

Step 2.

: if then and , if then and , if then and (add arcs to and from the measurable places according to its associated vector event ).

Step 3. If component of vector is 1 then , otherwise (put tokens in the measurable places to represent the first output vector).

Step 4. Associate to each the input symbol registered at the detection of e j .

The net with measurable places and inputs for the example is showed in figure 4.

2) Model simplifying

Implicit non-measurable can be removed: if there is a nonmeasurable place whose input and output transitions are exactly the same than any measurable place, then delete such a place and its input and output arcs.

The final model for the illustrative example is showed in figure 5; the incidence matrix, output function and initial marking of this IPN are given below. The associated inputs for transitions is given by the lambda function: (t 1 )=a, (t 2 )= , (t 3 )=b, (t 4 )=d, (t 5 )=c, (t 6 )=e.

, , Since every one of transitions in the net actually represents a sequence of events of length , the output language of length + 1 of the net is equal to the observed output language. Even, for the illustrative example, the output language of the IPN is equal to the observed output language, i.e. only the observed cyclic output sequences are represented by the evolution of the net.

F. Procedure

Now, we can summarise the stages of the method for IPN model identification. 2. Compute the event sequences from the observed vectors.

3. For every sequence of events, create traces of length .

4. Create the non-observable behaviour of the IPN (Algorithm 1) and simplify it (algorithm 2).

5. Complete the Petri net adding the observed behaviour and delete implicit places (algorithm 3).

Proposition 1. For a DES S that holds the hypothesis given in section III.A and an identification parameter , Algorithm 4 yields an IPN model (Q,M 0 ) which represents exactly .

Proof. Since the deletion of implicit places does not alter , we make the proof with the model obtained before this procedure. The firing of a transition in the system is not affected by the addition of arcs to, and from , since these arcs were computed from differences of vectors in . Then, also in this model, every event sequence of length less or equal than belongs to the language of the net iff it was observed.

The sequences of transitions of length less or equal than that can be fired lead to markings in the measurable places that also belong to (since the marking change provoked in the measurable places was obtained from the difference of observed vectors). Then, we have that sequences of observed output vectors of length less or equal than correspond to sequences of marking vectors in the net and .

IV. CONCLUSION

In this paper we provided a method to create an IPN model for a given DES described by a set of input-output sequences that may contain cycles representing infinite length behaviour. The proposed off-line procedure operates considering an accuracy parameter such that in the output language of the obtained IPN only and all observed output sequences of length + 1 are represented. This approach avoids including in the model more sequences than that observed.

The proposed algorithms have polynomial complexity on the input data size; they have been tested with numerous examples of diverse complexity including an industrial experimental system; in general must be greater than one and the obtained model allows representing complex behaviour such as concurrency and non determinism.

Current research deals with a depth study on the choice of the accuracy parameter and the definition of an efficient procedure for the analysis of concurrency in the basic nonmeasurable model.

  It has five output signals, = {A, B, C, D, E}, and five input signals = {a, b, c, d, e}. Two output sequences have been observed:

Algorithm 1 .

 1 Building the basic IPN structure Input: The set Output: An IPN model G composed by non-measurable places Step1.

Figure 1 .Algorithm 2 .

 12 Figure 1. Basic model describing the sequences of event traces 2) Simplifying the basic structure Additional node merging operations can be performed on the basic structure in order to obtain an equivalent trace model. Now we can take into account the event e j associated to transitions. Consider the following transformation rules. Algorithm 2. Simplifying basic structure. Input: G Output: G': an equivalent IPN Apply the following rule on the initial place Rule 1: //If has several output

  firing; this behaviour is exhibited by concurrent transitions. The analysis can be performed on a model component comprised between two transitions relied by several paths containing the concurrent transitions. If there are m! paths, we can explore if there exists m different transitions in the paths and every sequence is a permutation from each other. When it is verified, the subnet can be transformed into a concurrent component of G' preserving the same behaviour.

Figure 2 .

 2 Figure 2. Model after merging

Figure 3 .

 3 Figure 3. Simplified basic model

Figure 4 .

 4 Figure 4. IPN model including the measurable places

Figure 5 .

 5 Figure 5. Simplified IPN model

Algorithm 4 .

 4 Building an IPN model Inputs: A DES and an accuracy parameter Output: (Q,M 0 ): an IPN model 1. Obtain the input symbols and the cyclic sequences of observed output vectors.

Definition 1 :

 1 An ordinary Petri Net G is a bipartite digraph represented by the 4-tuple G=(P,T,I,O) where:

	P = {p 1 ,p 2 ,...,p n } and T ={t 1 ,t 2 ,...,t m } are finite sets of
	vertices named places and transitions respectively,
	I (O) : P × T	{0, 1}is a function representing the arcs
	going from places to transitions (transitions to places).
	Pictorially, places are represented by circles, transitions are
	represented by rectangles, and arcs are depicted as arrows. The
	symbol t j (t j ) denotes the set of all places p i such that I(p i ,t j ) 0
	(O(p i ,t j ) 0). Such places are called input (output) places of t j .
	Analogously, p i (p i ) denotes the set of all transitions t j such
	that O(p i ,t j ) 0 (I(p i ,t j ) 0). Such transitions are called input
	(output) transitions of p i .
	The incidence matrix of G is	, where
	;

ACKNOWLEDGEMENT

A.P Estrada-Vargas was sponsored by CONACYT Grant No. 50312, and by ENS de Cachan International Scholarship.