R Golovchak

A Kozdras

T Hodge

J Szlęzak

C Boussard-Pledel

Ya Shpotyuk
email: yashpotyuk@gmail.com

B Bureau

Optical and thermal properties of Sb/Bi-modified mixed Ge-Ga-Se-Te glasses

Keywords: Chalcogenide glass, Crystallization, Kinetics, DSC, Optical properties

 glasses are studied with differential scanning calorimetry method, and analysed using Frazer-Suzuki fitting function. The applicability of Johnson-Mehl-Avrami (JMA) model to describe the crystallization kinetics in these materials is verified, and JMA exponent is determined for the relevant crystallization processes. It is found, that substitution of Sb with Bi in glass composition leads to a decrease in optical gap and thermal stability of the glass, and increase in the density and crystallization activation energy for glass-ceramics formation.

Introduction

In a past decade, the compounds of chemical elements with large spin-orbit electron coupling attracted a substantial attention because of the possibility of quantum Hall effect [START_REF] Moore | The birth of topological insulators[END_REF][START_REF] Ando | Topological insulator materials[END_REF][START_REF] Wehling | Dirac materials[END_REF].

In the latter case, the spin-orbit coupling plays the same role as the external magnetic field: instead of being driven by a magnetic field, the electrons can form a quantum Hall state as a result of spindependent force when moving through the crystal lattice even in non-magnetic substances [START_REF] Ando | Topological insulator materials[END_REF][START_REF] Wehling | Dirac materials[END_REF][START_REF] Agarwala | Topological Insulators in Amorphous Systems[END_REF][START_REF] Mitchell | Amorphous topological insulators constructed from random point sets[END_REF].

Materials which possess this effect have been named topological insulators, because they are typically insulators in the 'bulk' but have exotic metallic states present at their surfaces due to the topological order (these surface states are one of the most important properties of the topological insulators) [START_REF] Moore | The birth of topological insulators[END_REF][START_REF] Ando | Topological insulator materials[END_REF][START_REF] Wehling | Dirac materials[END_REF]. For a topological insulator to form, the spin-orbit coupling must be strong enough to significantly modify the electronic structure -it should in fact be comparable to the bandgap of the material [START_REF] Moore | The birth of topological insulators[END_REF][START_REF] Ando | Topological insulator materials[END_REF][START_REF] Wehling | Dirac materials[END_REF]. Historically, one of the first discovered materials with topological insulating properties was the BixSb1-x compound [START_REF] Hsieh | Observation of Unconventional Quantum Spin Textures in Topological Insulators[END_REF]. Recently, some chalcogenide stoichiometric crystals, such as Bi2Te3, Bi2Se3, Sb2Te3 and Sb2Se3, as well as mixed/doped derived compounds have been proposed as perspective three-dimensional topological insulators [START_REF] Zhang | Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface[END_REF][START_REF] Chen | Experimental realization of a three-dimensional topological insulator[END_REF][START_REF] Seifert | Spin Hall photoconductance in a three-dimensional topological insulator at room temperature[END_REF][START_REF] Qin | Confined-path interference suppressed quantum correction on weak antilocalization effect in a BiSbTeSe2 topological insulator[END_REF][START_REF] Cho | Suppression of bulk conductivity and large phase relaxation length in topological insulator Bi2-dSndTe3 epitaxial thin films grown by Metal-Organic Chemical Vapor Deposition (MOCVD)[END_REF][START_REF] Usanmaz | Spinodal Superlattices of Topological Insulators[END_REF]. Some of them are characterized by extremely simple topological surface states described by a single gapless Dirac cone at the k = 0 Γ point in the surface Brillouin zone [START_REF] Zhang | Topological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surface[END_REF][START_REF] Chen | Experimental realization of a three-dimensional topological insulator[END_REF][START_REF] Usanmaz | Spinodal Superlattices of Topological Insulators[END_REF]. It is obvious, that Bi and Sb atoms play a key role in the formation of topological insulating properties.

Besides the prospect of forming topological insulators, incorporation of Bi or Sb into the structure of chalcogenides has a profound impact on their electrical properties [START_REF] Vautier | Role of metal impurity "Bi" in amorphous chalcogenide semiconductors[END_REF][START_REF] Feltz | Amorphous inorganic materials and glasses[END_REF]. These dopants are able to unpin the Fermi level and change the conduction from p to n type in chalcogenide glasses (which are known to be highly insensitive to metal doping [START_REF] Vautier | Role of metal impurity "Bi" in amorphous chalcogenide semiconductors[END_REF][START_REF] Feltz | Amorphous inorganic materials and glasses[END_REF]), and even cause a strong thermoelectric power effect [START_REF] Singh | Effect of bismuth on the electrical properties of a-Ge20Se80 glasses[END_REF][START_REF] Liu | A Study on Crystallization Kinetics of Thermoelectric Bi2Se3 Crystals in Ge-Se-Bi Chalcogenide Glasses by Differential Scanning Calorimeter[END_REF]. The combination of these properties with high optical transparency (extending up to ~18 µm) and non-linearity [START_REF] Zakery | Optical nonlinearities in chalcogenide glasses and their applications[END_REF][START_REF] Sanghera | Active and passive chalcogenide glass optical fibers for IR applications: a review[END_REF], makes Bi/Sb-containing chalcogenide alloys a very appealing multifunctional medium for numerous applications in (nano)optics and (nano)electronics [START_REF] Eggleton | Chalcogenide photonics[END_REF][START_REF]Chalcogenide Glasses: Preparation, properties and application[END_REF][START_REF] Lim | Thermoelectric properties of thin film topological insulators: A firstprinciples study[END_REF][START_REF] Liu | Crystallographically Textured Nanomaterials Produced from the Liquid Phase Sintering of BixSb2-xTe3 Nanocrystal Building Blocks[END_REF]. Thus, it was discovered recently that small variation in Bi content allows fabrication of ~mm thick fully transparent or fully opaque in ~3-16 µm wavelengths range semiconducting optical media [START_REF] Golovchak | Giant VIS-IR light attenuation effect in nanostructured narrow-bandgap glasses[END_REF], which can be used for a variety of applications in IR optics and sensors. The possibility to form Bi2Te3, Bi2Se3, Sb2Te3 and Sb2Se3 (nano)crystals in amorphous matrix by a controlled crystallization (heating the glass above its glass transition temperature) or laser exposure (inducing induced fluidity or local melting of the material), opens a whole range of new metamaterials with unique physical properties [START_REF] Agarwala | Topological Insulators in Amorphous Systems[END_REF][START_REF] Mitchell | Amorphous topological insulators constructed from random point sets[END_REF][START_REF] Lim | Thermoelectric properties of thin film topological insulators: A firstprinciples study[END_REF][START_REF] Liu | Crystallographically Textured Nanomaterials Produced from the Liquid Phase Sintering of BixSb2-xTe3 Nanocrystal Building Blocks[END_REF].

In this paper, we introduce new glassy materials within Bi-Sb-Ga-Ge-Se-Te family with varying Bi and Sb content, which can be used as universal host matrix having potential to combine superior optical properties (like IR transparency), thermo-electric effect, different types of conductivity and topological insulating behaviour in one medium. The optical and thermal properties for these complex metamaterials are analyzed and discussed.

Materials and methods

The BixGa5Ge20Sb10-xSe45Te20 (x=1,3,5,10) glasses were prepared by conventional meltquenching technique using high-purity elemental precursors (5N or more). Appropriate amounts of Ga, Bi, Sb, Ge, Se and Te (around 20 g in total) were vacuum-sealed in silica tube of 10 mm diameter. Then, the ampoules were heated up to 900 °C with 2 o C/min rate in a rocking furnace, homogenized for ~10-12 hours, quenched into room-temperature water from 700 °C and swiftly moved to the preheated furnace for annealing at ~30 °C below the glass transition temperature for 3 h. For optical measurements samples were cut into discs and polished to thickness of ~1.8 mm.

The overall glassy state of the obtained materials was confirmed by the absence of crystalline reflexes in their X-ray diffraction (XRD) patterns (Fig. 1a). Nevertheless, some minor reflexes have been noticed in the XRD pattern of Bi10Ga5Ge20Se45Te20 sample (Fig. 1a), which could be attributed to the formation of small-size nucleation sites right after the quench. The homogeneity and actual compositions of the prepared materials Ga5Ge20Sb10Se45Te20 (Bi0), Bi1Ga5Ge20Sb9Se45Te20 (Bi1), Bi3Ga5Ge20Sb7Se45Te20 (Bi3), Bi5Ga5Ge20Sb5Se45Te20 (Bi5) and Bi10Ga5Ge20Se45Te20 (Bi10) were verified by scanning electron microscopy (SEM) with energy-dispersive X-ray analysis (EDX), and X-ray photoelectron spectroscopy. No significant impurities or deviations (more than ±1 at.%) from nominal compositions have been found, as well as no significant inhomogeneities have been detected through SEM (Fig. 1b).

The high-energy XRD experiment was performed at the beamline 11-ID-B, Advanced Photon Source, Argonne National Laboratory, using X-rays with wavelength 0.2127 Å and position sensitive Perkin Elmer amorphous silicon image plate detector [START_REF] Chupas | Rapidacquisition pair distribution function (RA-PDF) analysis[END_REF]. These measurements were carried out at room temperature using coarse powdered (with mortar and pestle) samples sealed in Polyimide capillaries. Differential scanning calorimetry (DSC) measurements were performed using Mettler Toledo DSC1 instrument calibrated with Hg, In, Sn and Zn as standard elements. Bulk ~20 mg chunks of each glass were heated in 40 µl aluminum crucibles under nitrogen atmosphere at 2, 5, 10, 15, and 20 K/min constant heating rates. To assure the repeatability of the results, the DSC measurements were repeated minimum 3 times for each heating rate and for each composition (each time using a fresh sample). The DSC scan included 2 runs at the same heating rate: the first run revealed the glass-to-supercooled liquid transition and crystallization peaks; the second run of fully crystallized glass-ceramic sample provided us with a confident baseline, which was then subtracted from the first-run curve. The crystallization peaks in the repeatable DSC data sets were fitted with Fraser-Suzuki distributions [START_REF] Perejon | Kinetic Analysis of Complex Solid-State Reactions. A New Deconvolution Procedure[END_REF] in Origin. Each of these fitted curves was then integrated and normalized using Origin's intrinsic functions to obtain a fraction of crystallized volume (α).

Optical transmission spectra in the IR region were obtained with Bruker Vertex 70V spectrometer, while for the measurements in UV/VIS/NIR region the PerkinElmer LAMBDA 950 spectrophotometer was used. The VIS/NIR spectra and PARAV software [START_REF] Ganjoo | Computer program PARAV for calculating optical constants of thin films and bulk materials: Case study of amorphous semiconductors[END_REF] were used to estimate optical gaps of the glasses.

Results and discussion

Optical transmission spectra (Fig. 2a) and absorption coefficient, calculated in the fundamental optical absorption edge region (Fig. 2b), are typical for Se/Te-based glasses [START_REF]Chalcogenide Glasses: Preparation, properties and application[END_REF][START_REF] Cui | Selenide and telluride glasses for mid-infrared bio-sensing[END_REF].

The samples with low Bi content are characterized by a relatively wide transmitance window (~2-16 µm), which allows a variety of IR applications. Complete substitution of Sb with Bi in glass composition leads to a full opacity of Sb-free Bi10 sample in a wide VIS-IR spectral range (transmission less than 1% for ~2 mm thick sample, not shown in Fig. 2a). This opacity is not caused by scattering, though, since no crystallites of appreciable size or other inclusions have been found in this Bi10 sample neither with XRD nor with SEM. The mechanism is, the most probably, associated with giant attenuation coefficient as suggested for similar glasses earlier [START_REF] Golovchak | Giant VIS-IR light attenuation effect in nanostructured narrow-bandgap glasses[END_REF]. The optical gap estimated from Tauc plots using PARAV software [START_REF] Ganjoo | Computer program PARAV for calculating optical constants of thin films and bulk materials: Case study of amorphous semiconductors[END_REF] decreases with Bi content from the value of ~1 eV for Bi0 glass to ~0.7 eV for Bi5 sample (Table 1).

The DSC thermograms of the investigated glasses recorded at different heating rates q = 2, 5, 10, 20, 30 K/min show glass transition (Fig. 3) and crystallization (Fig. 4) features in all cases.

The onset temperature of the glass transition Tg on and the peak temperature of the first crystallization process Tc (both values determined for 10 K/min DSC runs) were used to calculate Dietzel criterion of glass stability ∆T = Tc -Tg on [START_REF] Dietzel | Glass Structure and Glass Properties[END_REF] shown in Table 1. It is found that Tg value as determined from DSC heating curves does not change significantly with Bi addition (Tg = 230±2 o C for 10 K/min runs), except a noticable decrease for Bi10 sample (Tg ~ 220 o C). The activation energy of viscous flow, ETg, detemined by Ozawa method [START_REF] Ozawa | Nonisothermal kinetics of crystal growth from pre-existing nuclei[END_REF][START_REF] Ozawa | A new method of analyzing thermogravinetric data[END_REF] from the data in Fig. 3 for the thermally annealed (structurally relaxed) glasses, demonstrates a maximum in compositional dependence at ~1-2 at.% of Bi (Table 1). Simultaneousely, the density of glass measured by Archimedes' method increases (Table 1) and the glass becomes less stable (smaller ∆T values) when Sb is substituted with Bi, the most unstable being the Bi10 composition with 10 at.% of Bi content (its still a glass, though, according to XRD data). Nevertheless, the Tg of all investigated samples (including Bi10) is higher than 200 o C, which makes these glasses attractive for applications in many IR instruments.

For the samples with low Bi content (x < 5 at.%) one family of crystallization peaks is observed, while for Bi5 and Bi10 samples the low-temperature crystallization emerges (Fig. 4).

Fraser-Suzuki fits of DSC heating curves (Fig. 5) suggests the crystallization of at least 2 phases for the Bi-free sample (Bi0 composition), 3 phases for intermediate Bi concentrations (Bi1, Bi3 and Bi5 compositions); and as many as 5 phases for Bi10 sample. It should be noted however, that shape of real crystallization peaks can be more complicated than Fraser-Suzuki form and usually depends on the crystallizing phase [START_REF] Perejon | Kinetic Analysis of Complex Solid-State Reactions. A New Deconvolution Procedure[END_REF]. So, the proposed fits should be considered as idealized models with precautions, especially if there are weak shoulders which can be obscured by assymmetric form of Frazer-Suzuki function. They are, however, our best estimates (due to the overlapped crystallization processes of unknown nature) to determine the crystallization temperatures and activation energies of crystallization processes.

Johnson-Mehl-Avrami (JMA) nucleation-growth model is usually used to study the nonisothermal crystallization kinetics with DSC [START_REF] Malek | Kinetic analysis of crystallization processes in amorphous materials[END_REF][START_REF] Sestak | Thermophysical Properties of Solids. Their Measurements and Theoretical Analysis[END_REF]. The fraction of crystallized volume α during non-isothermal crystal growth from preexisting nuclei can be determined from [START_REF] Malek | Kinetic analysis of crystallization processes in amorphous materials[END_REF][START_REF] Sestak | Thermophysical Properties of Solids. Their Measurements and Theoretical Analysis[END_REF]:

      - = RT E a e Af dt d ) (α α (1) ∫ ∆ = T c dT H 0 1 φ α . ( 2
)
Here φ is the specific heat flow measured with DSC (W/g) and ΔHc is the total enthalpy change associated with the crystallization process; the pre-exponential factor A and activation energy Ea are kinetic parameters that should not depend on the temperature T and α; and

m m f / 1 1 )] 1 ln( )[ 1 ( ) ( - - - - = α α α (3)
is an algebraic expression of the JMA model [START_REF] Malek | Kinetic analysis of crystallization processes in amorphous materials[END_REF][START_REF] Sestak | Thermophysical Properties of Solids. Their Measurements and Theoretical Analysis[END_REF].

The JMA exponent m is a characteristic parameter related to a crystal forming morphology: m ~ 1 means predominant surface crystallization, m ~ 3 corresponds to three-dimensional bulk crystallization, and m ≥ 4 usually describes the crystallization processes occurring with increasing nucleation rate [START_REF] Matusita | Kinetics of non-isothermal crystallization process and activation energy for crystal growth in amorphous materials[END_REF].

A simple test for the applicability of JMA model is based on the analysis of probe functions [START_REF] Malek | Kinetic analysis of crystallization processes in amorphous materials[END_REF]:

      - = RT E a e y φ α ) ( (4) 2 ) ( T z φ α = (5)
For JMA model to be valid, the maximum of z(α) function should occur around α = 0.63±0.02 value in the non-isothermal crystallization conditions [START_REF] Malek | Kinetic analysis of crystallization processes in amorphous materials[END_REF].

The fits with Fraser-Suzuki distribution function were used to build z(α) plots for the resolved crystallization processes. As could be seen from Fig. 6, only peaks III of Bi1 and Bi3 samples, as well as peak II of Bi5 sample have their z(α) maxima at roughly α ~ 0.63 value (taking into account a reasonable accuracy achievable with DSC instrument for rapidely crystallizing phases). So, only these crystallization processes can be potentially described with JMA model of crystallization. The JMA exponent m estimated from the double-logarithm plots (Fig. 7) is found to be close to 3 for the peaks III of Bi1 and Bi3 samples, and approximately 7 for the peak II of Bi5 sample. The latter means crystallization processes with increasing nucleation rate [START_REF] Matusita | Kinetics of non-isothermal crystallization process and activation energy for crystal growth in amorphous materials[END_REF].

The activation energies of crystallization (Ea) calculated using Ozawa's plots (ln q vs 1000/T) [START_REF] Ozawa | Nonisothermal kinetics of crystal growth from pre-existing nuclei[END_REF][START_REF] Ozawa | A new method of analyzing thermogravinetric data[END_REF] show overall increase with Bi content (Table 2). This can be explained by the fact that onset of crystallization shifts to lower temperatures in Bi-contaning glasses, where higher viscosity of supercooled liquid imposes more constraints on structural rearrangements needed for crystallization to occur. Simultaneously, crystallization peaks become sharper (Fig. 4), which means more rapid crystallization of the phases in glasses with greater Bi content. This result is of significant importance for the application of these materials in phase-change memory devices.

The nature of the crystallizing phases in each case requires further thorough investigations using combination of Raman/IR spectrocopy, electron microscopy, XRD analysis, and such complementary structural probes as X-ray photoelectron spectroscopy and extended X-ray absorption fine structure spectroscopy. Nevertheless, on the basis of previous studies [START_REF] Saturday | Devitrification of Bi-and Ga-containing germanium-based chalcogenide glasses[END_REF], we expect that in the low temperature region the predominant crystallization of Bi-based compounds (like Bi2Te3, Bi2Se3 or Bi2SeTe2) should occur. In the high-temperature range beyond ~320 o C, the crystallization processes should rely mostly on the crystallization of stable Ge-and Sb-based chalcogenide phases [START_REF] Saturday | Devitrification of Bi-and Ga-containing germanium-based chalcogenide glasses[END_REF].

Conclusions

The Ga-Ge-Se-Te matrix is found to be quite suitable for the incorporation of high concentration of Sb and Bi (up to 10 at.%), while preserving the overall glassy state. Substitution of Sb with Bi leads to an increase in density of the glass, a decrease in the optical gap of these materials and, finally, to the complete opacity of ~2 mm thick sample with 10 at.% of Bi in a wide visible and infrared range of spectrum. Analysis of the crystallization kinetics reveals that only few crystallization processes can be described with Johnson-Mehl-Avrami crystallization model, while majority of them do not pass JMA applicability test (at least per Frazer-Suzuki fits). The activation energies of crystallization show a tendency to increase with Bi content, which is associated with an overall shift of the crystallization processes to lower temperatures. Glassceramics formation by crystallization becomes more rapid in Bi-rich glasses, but simultaneously their Dietzel criterion of glass stability decreases. Owing to the obtained results, the developed glass matrix can be considered for multifunctional applications in IR optics, photonics, phasechange memory devices, and sensors. Table 1. Dietzel criterion of glass stability ∆T = Tc -Tg on [START_REF] Dietzel | Glass Structure and Glass Properties[END_REF], calculated for the investigated BixGa5Ge20Sb10-xSe45Te20 glasses using DSC curves recorded at 10 K/min heating rate (Tc is the peak value of the earliest crystallization peak, Tg on is the onset value of glass transition); density

  ρ; estimated optical gap Eg; and activation energies of viscous flow ETg, averaged for Ozawa's plots of onset and midpoint Tg values.

Fig. 1 .

 1 Fig. 1. (a) High-energy XRD patterns of the investigated BixGa5Ge20Sb10-xSe45Te20 glasses; (b)

Fig. 2 .

 2 Fig. 2. Transmission (a) and fundamental optical absorption (b) spectra of

Fig. 3 Fig. 3 .

 33 Fig.3

Fig. 4 .

 4 Fig. 4. DSC scans at different heating rates for BixGa5Ge20Sb10-xSe45Te20 glasses, showing

Table 2 .

 2 Average activation energy of crystallization Ea calculated for the investigated materials using Ozawa's plots[START_REF] Ozawa | Nonisothermal kinetics of crystal growth from pre-existing nuclei[END_REF][START_REF] Ozawa | A new method of analyzing thermogravinetric data[END_REF].

	(±0.005)

Acknowledgements

Work performed at Argonne National Lab and use of the Advanced Photon Source were supported by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357. JS acknowledge support of this work provided by CampusFrance via BGF joint-PhD Cotutelle scholarship programme. APSU team acknowledges help with high-energy XRD measurements Joerg C. Neuefeind, Katharine Page, Matthew Tucker and Jue Liu. The authors acknowledge support for this research from NSF Grant DMR-1409160.


Frédéric J P Richard

Anisotropy of Hölder Gaussian random elds: characterization, estimation, and application to image textures

Keywords: Hölder regularity, anisotropy, fractional Brownian eld, quadratic variations, texture analysis, photographic paper

The characterization and estimation of the Hölder regularity of random elds has long been an important topic of Probability theory and Statistics. This notion of regularity has also been widely used in Image Analysis to measure the roughness of textures. However, such a measure is rarely sucient to characterize textures as it does not account for their directional properties (e.g. isotropy and anisotropy). In this paper, we present an approach to further characterize directional properties associated to the Hölder regularity of random elds. Using the spectral density, we dene a notion of asymptotic topothesy which quanties directional contributions of eld high-frequencies to the Hölder regularity. This notion is related to the topothesy function of the so-called anisotropic fractional Brownian elds, but is dened in a more generic framework of intrinsic random elds. We then propose a method based on multi-oriented quadratic variations to estimate this asymptotic topothesy. Eventually, we evaluate this method on synthetic data and apply it for the characterization of historical photographic papers.

Introduction

In this paper, we focus on irregular Gaussian random elds (called Hölder random elds) whose realizations are continuous but non-dierentiable (see Section 2 for more details). The degree of Hölder regularity of these elds is quantied by a parameter H (called the Hölder index) in (0, 1). Hölderian elds include fractional Brownian elds (i.e. multidimensional versions of fractional Brownian motions [START_REF] Mandelbrot | Fractional Brownian motion, fractionnal noises and applications[END_REF], and their F. Richard Aix Marseille Univ, CNRS, Centrale Marseille, I2M, Marseille, France. E-mail: frederic.richard@univ-amu.fr anisotropic extensions [START_REF] Biermé | Hölder regularity for operator scaling stable random elds[END_REF][START_REF] Biermé | Operator scaling stable random elds[END_REF][START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF][START_REF] Clausel | Explicit construction of operator scaling Gaussian random elds[END_REF]Richard, 2016b). They have been widely used in Image Analysis to model rough image textures from engineering domains such as Medical Imaging (Biermé et al., 2009;[START_REF] Biermé | Analysis of texture anisotropy based on some Gaussian elds with spectral density[END_REF]Richard, 2016a;[START_REF] Richard | Statistical tests of anisotropy for fractional Brownian textures. Application to full-eld digital mammography[END_REF]. The Hölder index of these models has served for the quantication of texture roughness. [START_REF] Biermé | A turning-band method for the simulation of anisotropic fractional Brownian eld[END_REF]. Fields of rst and second rows have Hölder indices of 0.3 and 0.6, respectively. Fields of a same column have a topothesy function which is represented on the third row. Hurst functions of elds are constant and equal to the Hölder index. All simulations were obtained using a same pseudo-random number sequence so as to highlight texture dissimilarities due to variations of the simulated elds.

However, this index is not always sucient to characterize the texture aspect. In particular, since it is regardless of orientations, it does not account for directional properties of textures. This shortcoming is illustrated on Figure 1 with some simulated textures having both low and high Hölder indices. From a regularity viewpoint, textures of the two dierent rows can be distinguished while those of a same row cannot. Dierences between textures of a same row are only due to variations of their directional properties. In particular, the rst texture of each row is isotropic (i.e. it has same aspect in all directions) whereas the second and third ones are anisotropic (i.e. their aspect varies depending on the direction).

The main motivation for this work is to set a description of textures that would not only account for their Hölder regularity but also for relevant directional properties associated to this regularity. In that perspective, we rst propose a characterization of the anisotropy of Hölder elds. Then, we address the issue of estimating features derived from this characterization.

The anisotropy of Hölder elds is often characterized through parameters of a specic model [START_REF] Biermé | Operator scaling stable random elds[END_REF][START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF][START_REF] Clausel | Explicit construction of operator scaling Gaussian random elds[END_REF][START_REF] Roux | Self-similar anisotropic texture analysis: The hyperbolic wavelet transform contribution[END_REF]. As an example, anisotropic fractional Brownian elds (AFBF) introduced by [START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF] are ddimensional Gaussian elds with stationary increments whose second-order properties are determined by a spectral density of a form (see Section 2 for details) g τ,η (w) = τ (arg(w))|w| -2η(arg(w))-d .

(1) Such a density depends on two functions τ and η called the topothesy and Hurst functions, respectively. By assumption, these functions are both even, positive, π-periodic, and bounded functions. Depending on the spectral orientation arg(w), they can be used to characterize directional properties of AFBF; visual eects induced by the topothesy function are illustrated on Figure 1. However, such a characterization is specic to a model. Moreover, it concerns all eld frequencies (including its low frequencies) and is not exclusively associated to the eld regularity. More generic characterizations which are intrinsically linked to a notion of regularity were developed by Abry et al. (2015a); [START_REF] Hochmuth | Wavelet characterizations for anisotropic Besov spaces[END_REF]; Slimane and Braiek (2012). They rely upon the analysis of an anisotropic function space (typically, an anisotropic Besov space) by decomposition of elds into an appropriate basis (e.g. a basis of the hyperbolic wavelets).

In this paper, we investigate another characterization approach which is based on the eld spectral density. In a generic framework of Hölder random elds introduced by Richard (2016b), we dene a function, called the asymptotic topothesy, which quanties directional contributions of high-frequencies of a eld to its irregularity. This function is related to the topothesy of an AFBF. For such a eld, eld high-frequencies in directions where the Hurst function reaches a minimal value H are the largest. Due to these highfrequencies, the Hölder index of the eld is H (see Section 2 for details). The inuence of these high-frequencies on the eld regularity is further weighted by the topothesy function: the larger this function, the larger the high-frequencies and their contribution to the eld regularity.

In a second part, we focus on the estimation of the asymptotic topothesy. This issue is both functional and non-parametric. It diers from the estimation issues that are usually tackled in the literature. Indeed, the analysis of directional features associated to random elds often reduces to the estimation of parameters of a specic anisotropic model. For instance, the estimation procedure by [START_REF] Roux | Self-similar anisotropic texture analysis: The hyperbolic wavelet transform contribution[END_REF], which is an implementation of the anisotropy characterization by Abry et al. (2015a) using hyperbolic wavelets, targets the single parameter of a Gaussian operator scaling eld. The works of [START_REF] Biermé | Estimation of anisotropic Gaussian elds through Radon transform[END_REF]; [START_REF] Richard | Statistical tests of anisotropy for fractional Brownian textures. Application to full-eld digital mammography[END_REF] are devoted to the Hurst function of an AFBF. In these works, the estimation procedure is based on a Radon transform of images. Due to the discretization of this transform, the procedure can only be applied to the estimation of the Hurst function in a few directions. [START_REF] Richard | Analysis of anisotropic Brownian textures and application to lesion detection in mammograms[END_REF]Richard ( , 2016a,b) ,b) have already investigated an estimation issue within the random eld framework of this paper. In these works, an anisotropy analysis was developed using the so-called multi-oriented quadratic variations which are sums of squares of eld increments computed in dierent directions. This analysis relies upon the estimation of a directional function which is asymptotically and linearly linked to quadratic variations. In this paper, we show that this directional function is indirectly related to the asymptotic topothesy through a convolution with a specic kernel that we analytically compute. We thus propose to recover the asymptotic topothesy of a eld by solving an inverse problem associated to this convolution.

Eventually, we illustrate the interest of the asymptotic topothesy on an application to photographic papers. For the description of textures of these papers, we test a combination of two indices, an estimated Hölder index and an anisotropy index derived from the estimated asymptotic topothesy.

Theoretical Framework

Hölder Regularity

The denition of the asymptotic topothesy is associated to the notion of Hölder regularity which is dened as follows [START_REF] Adler | The geometry of random elds[END_REF].

Denition 1 A eld Z satises a uniform stochastic Hölder condition of order α ∈ (0, 1) if, for any compact set C ⊂ R d , there exists an almost surely nite, positive random variable A such that the Hölder condition

|Z(x) -Z(y)| ≤ A|x -y| α .
(2) holds for any x, y ∈ C, with probability one. If there exists H ∈ (0, 1) for which Condition (2) holds for any α < H but not for α > H, then we say that Z is Hölder of order H or H-Hölder. The critical parameter H is called the Hölder index.

Hölder elds are well-suited for the modelling of images with rough textures. In such modelling, we can interpret the texture as a visual eect of the eld irregularity. Then, the texture roughness depends on the degree of the eld regularity, and may be quantied from 0 to 1 by 1 -H.

For random elds with stationary increments such as the AFBF, the Hölder regularity can be characterized from the behavior of the spectral density at high-frequencies [START_REF] Bonami | Anisotropic analysis of some Gaussian models[END_REF].

Proposition 1 Let Z is a mean square continuous Gaussian with stationary increments. Assume that its semi-variogram is determined by a spectral density f as follows

1 2 E((Z(y + h) -Z(y) 2 ) = 1 (2π) d R d f (w)|e i h,w -1| 2 dw.
Let H ∈ (0, 1). (i) If for any 0 < α < H, there exist two positive constants A 1 and B 1 such that for almost all

w ∈ R d |w| ≥ A 1 ⇒ f (w) ≤ B 1 |w| -2α-d , (3)
then the eld Z is Hölder of order ≥ H.

(ii) If, for any H < β < 1, there exist two positive constants A 2 and B 2 and a positive measure subset

E of [0, 2π) d-1 such that |w| ≥ A 2 and arg(w) ∈ E ⇒ f (w) ≥ B 2 |w| -2β-d , (4)
then the eld Z is Hölder of order ≤ H.

(iii) If conditions (3) and (4) both hold, then the eld Z is H-Hölder.

Example 1 (Anisotropic fractional Brownian elds)

Let Z be an AFBF as dened by Equation (1). Assume that

H = ess inf{η(s) ∈ (0, 1), s ∈ [0, 2π) d-1 , τ (s) > 0} ∈ (0, 1), (5)
Then, the eld is Hölder of order H. More generally, assume that the eld density only satises

|w| > A ⇒ 0 ≤ f (w) -g τ,η (w) ≤ C|w| -2H-d-γ , (6)
for some positive constants A, C and a spectral density g τ,η of the form (1) [START_REF] Richard | Analysis of anisotropic Brownian textures and application to lesion detection in mammograms[END_REF](Richard, , 2016a,b),b), then it is Hölder of an order H given by (5).

According to Propositions 2.1.6 and 2.1.7 of [START_REF] Biermé | Champs aléatoires : autosimilarité, anisotropie et étude directionnelle[END_REF], Proposition 1 is still valid for intrinsic random elds (IRF) which extend elds with stationary increments (see [START_REF] Chilès | Geostatistics: modeling spatial uncertainty[END_REF] for an introduction). In what follows, we will work in the same framework of [START_REF] Richard | Analysis of anisotropic Brownian textures and application to lesion detection in mammograms[END_REF]Richard ( , 2016a,b) ,b) composed of IRF of an arbitrary order M whose density fullls Equation (6).

Asymptotic Topothesy

Denition 2 (Asymptotic topothesy) Let us assume that the density of a eld fullls Condition (6) for some H ∈ (0, 1). The asymptotic topothesy is a function dened, for almost all direction s of [0, 2π) d-1 as

τ * (s) = lim ρ→+∞ f (ρs)ρ 2H+d . ( 7
)
The asymptotic topothesy is a non-negative and bounded function τ * which basically gives a measure of the speed at which the spectral density converges to 0 at innity in each direction. Intuitively, such a measure quanties the magnitude of eld high-frequencies: the larger τ * in a direction, the slower the density convergence and the larger high-frequencies in that direction. Since τ * is bounded, the density convergence occurs at a speed which is not slower than a reference speed of order ρ -2H-d . According to Proposition 1 (i), this implies that the Hölder index of the eld cannot be below H. In directions s where τ * (s) = 0 the convergence speed is faster than any speed of order ρ -2H-d . Hence, in these directions, the condition of Proposition 1 (ii) is not satised. This means that high-frequencies of these directions are not large enough to make the Hölder index of the eld as low as H. By contrast, the convergence speed is of order ρ -2H-2 in directions of the set

E 0 = {s, τ * (s) > 0}. ( 8
)
On E 0 , the condition of Proposition 1 (ii) holds, which implies that the Hölder index is exactly H. In other words, high-frequencies in directions of E 0 contribute to the eld irregularity. Their contributions are further weighted by the asymptotic topothesy: the larger τ * in a direction of E 0 , the larger the contribution of high-frequencies to the eld irregularity.

The asymptotic topothesy linearly depends on the variance of the eld. Hence, in practice, we rather use a version of a topothesy normalized by its mean. This normalized topothesy measures relative contributions of highfrequencies to the eld irregularity, and provides us with a more intrinsic information about the eld anisotropy.

We will say that the eld texture is isotropic when these contributions are uniform (i.e. the asymptotic topothesy is constant), and anisotropic when they are not. Let us outline that these notions of isotropy and anisotropy describe high-frequencies of the eld, and are intrinsically related to its irregularity.

An IRF can always be decomposed into the sum of a smooth eld and a stationary eld which account for its low and high frequencies, respectively. Next, we state a proposition showing how the asymptotic topothesy characterizes the correlation structure of high-frequency part of this decomposition.

Proposition 2 Let Z be an IRF whose density f fullls Condition (6) for some Hölder index H ∈ (0, 1) and asymptotic topothesy τ * (see Equation ( 7)). Then, as A → +∞, for any

x ∈ R d , K A,f (x) -K A,f * (x) = o(A -2H ), (9)
where K A,f and K A,f * are two auto-covariances dened by

K A,g (x) = 1 (2π) d |w|>A e i x,w g(w)dw and f * (w) = τ * (arg(w))|w| -2H-d . ( 10
)
In this proposition, K A,f corresponds to the auto-covariance of the highfrequency stationary eld ZA of a decomposition of Z. Equation ( 9) means that, when A tends to ∞, the correlation structure of this eld gets approximately the same as the one of another stationary eld whose auto-covariance K * A only depends on the Hölder index and the asymptotic topothesy.

Estimation Method

In this section, we address the issue of the estimation of the asymptotic topothesy.

Multi-oriented Quadratic variations

Multi-oriented quadratic variations were introduced by Richard (2016b) to construct isotropy tests, and further used by [START_REF] Richard | Analysis of anisotropic Brownian textures and application to lesion detection in mammograms[END_REF]Richard ( , 2016a) ) to develop anisotropy indices. The denition of these variations is based on the computation of image increments. These increments give some information about image variations at highest observed scales which are not marked by trends. Moreover, as they are computed in dierent orientations, they provide us with relevant directional information.

Let us assume that an image is a realization of a random eld Z on a grid

[[1, N ]] d . Let us denote by Z N [m] = Z(m/N ) the image intensity at position m ∈ Z d . Given a vector u in Z d , increments in direction arg(u) at scale |u| are obtained by a discrete convolution ∀m ∈ Z d , V N u [m] = k∈Z 2 v[k]Z N [m -T u k], (11)
with an appropriate convolution kernel v and a transform T u which is a combination of a rotation of angle arg(u) and a rescaling of factor |u|. The kernel is chosen so as to ensure that the convolution annihilates any polynomial of a predened order K (kernel of order K) (Richard, 2016b).

Example 2 Some two-dimensional kernels selected by Richard (2016b) for their optimality are given for L ∈ N\{0} by

v[l 1 , l 2 ] = (-1) l1 L l 1 , (12) if (l 1 , l 2 ) ∈ [[0, L]] × {0} and 0 otherwise, n k standing for the binomial coe- cient. Such a kernel is of order K = L -1.
The information provided by increments are summarized into a single random variable called quadratic variations

W N u = 1 N e m∈E N (V N u [m]) 2 , ( 13
)
where E N is a set of cardinal N e containing positions m where increments can be computed on grid points. To get information at dierent scales and Richard orientations, we compute quadratic variations for dierent vectors u indexed in a set I of size n I . We gather all these variations into a single random vector Y N = log(W N u k ) k∈I of log-variations. The following theorem species the asymptotic probability distribution of Y N .

Theorem 1 For some integer M ≥ -1, let Z be a mean-square continuous Gaussian M -IRF. Assume that its spectral density f fullls Condition (6) for some H ∈ (0, 1). Let τ * asymptotic topothesy dened by Equation ( 7). Consider a log-variation vector Y N constructed using a kernel v of order K > M and K ≥ M/2 + d/4 if d > 4. For all i ∈ I, dene random variables N i such that

Y N i = H x N i + log(β H,τ * (arg(u i ))) + N i , (14)
with

x N i = log(|u i | 2 /N ) and β H,τ * (θ) = 1 (2π) d [0,2π) d-1 τ * (ϕ) Γ H,v (θ -ϕ) dϕ = τ * Γ H,v (θ), (15)
where stands for a circular convolution product over [0, 2π) d-1 , and Γ H,v is dened by

Γ H,v (θ) = R + |v (ρθ)| 2 ρ -2H-1 dρ, (16)
with v the discrete Fourier transform of v. Then, as N tends to +∞, the random vector (N

d 2 N i ) i∈I tends in distribution to a centred Gaussian vector.
This theorem is proved by Richard (2016b) (Theorem 3.4). The expression of β H,τ * in Equation ( 15) was slightly changed to highlight a convolution product.

Inverse Problem

In Theorem 1, the topothesy appears as a solution of Equation ( 15). In other words, it can be theoretically recovered by solving this equation. In practice, functions β H,τ * and Γ H,v involved in this equation are unobserved. However, they can both be estimated from log-variations. Indeed, according to Equation ( 14), log-variations are linearly related to the Hölder index H, which determines the closed form ( 16) of Γ H,v , and an intercept function which is equal to log(β H,τ * ). Hence, following Richard (2016b), we can estimate the Hölder index and the intercept function by linear regression on log-variations. Then, we can deduce estimates of β H,τ * and Γ H,v .

According to these remarks, the problem can be stated as follows. Let {θ j } j∈J be an indexed set of orientations in [0, 2π) d-1 formed by arguments of vectors {u i } i∈I . Let βj be some unbiased estimates of β H,τ * (θ j ) for j ∈ J , and H an unbiaised estimate of the Hölder index H. Then, for all j ∈ J , Theorem 1 implies that

βj = Γ H,v τ (θ j ) + ν j (17)
for some correlated Gaussian random variables ν j . Let Σ be the covariance matrix of the random vector formed by the ν j . One way to recover the topothesy τ would be to minimize a generalized least square criterion

C 0 (τ ) = i∈J j∈J Σ -1 ij βi -Γ H,v τ (θ i ) βj -Γ H,v τ (θ j ) (18)
over a space of real π-periodic functions of L 2 ([0, 2π)). Unfortunately, this minimization problem may be ill-posed; as we shall see in practice, it may lack stability, especially when H is close to 1. So as to x this issue, we choose a subspace W of L 2 ([0, 2π) d-1 ) equipped with an inner product and its associated norm | • | W . Over the space W , we then propose to minimize the penalized least square criterion

C λ (τ ) = C 0 (τ ) + λ|τ -τ 0 | 2 W , (19)
where λ > 0 and τ 0 is the mean value of τ over [0, 2π) d-1 . The second term of this criterion penalizes variations of τ by constraining the solution to remain as close as possible to a constant. The parameter λ sets a trade-o between this prior constraint and the delity of the recovered τ * to observations (as measured by C 0 ). In Section 3.4, we will present a method to set this parameter optimally.

Numerical Resolution

In order to minimize the criterion in Equation ( 19), we rst expand τ in an orthogonal basis (ϕ m ) m∈N of W with ϕ 0 ≡ 1 :

∀θ ∈ [0, 2π) d-1 , τ (θ) = m∈N τ m ϕ m (θ), (20)
for some scalar coecients τ m . From now on, we will use the same notation τ for the function and its expansion coecients. We then approximate the solution in a nite-dimensional subspace W a of W spanned by the rst a basis functions. On W a , the penalized criterion reduces to

Ca λ (τ ) = (L H τ -β) T Σ -1 (L H τ -β) + λ τ T Rτ, (21)
where β is a column vector containing estimates of β values, τ a vector gathering a + 1 coecients of decomposition in the basis of W a , L H a matrix of size |J | × a having terms ϕ m Γ H,v (θ j ) on the jth row and mth column, and R a diagonal matrix having (0,

|φ 1 | W , • • • , |φ a | W ) on the diagonal.
The minimum of the approximated criterion Ca λ is reached at

τ * λ = (L T H Σ -1 L H + λR) -1 Σ -1 L T H β. (22)
3.4 Choice of λ.

The solution τ * λ given by Equation ( 22) is intended to approach the solution τ * of the linear system

L H τ = β (23)
using some estimates β and H of β and H, respectively. Next, we give some bounds for the relative bias and variance of this estimation. To get these bounds, we neglect the eect of estimating H, and set L H L H ; in practice, such an approximation is compensated for by an accurate estimation of H (see Section 4).

Theorem 2 Let τ * be the solution of the linear system (23). Dene the matrix

B = L T H Σ -1 L H . Let ν -,
and ν + , κ be the lowest and highest eigenvalues, and the 2-norm condition number of B, respectively. Then, the relative bias and standard deviation of τ * λ taken as an estimator of τ * are bounded as follows.

BIAS = |E(τ * λ ) -τ * | |τ * | ≤ λ κ |R| 2 λ + ν + , (24)
STD = trace(V(τ * λ )) |τ * | ≤ κ ν + √ ν -trace(B -1 ) Σ -1 β, β 2 (λ + ν + ) . ( 25
)
Using this theorem, we can nd an optimal value of λ minimizing the bound of the relative mean square error.

Corollary 1 In the estimation of τ * by τ * λ , the relative mean square error

RMSE(λ) = E(|τ * λ -τ * | 2 ) |τ * | 2
is bounded by

g(λ) = κ 2 (λ + ν + ) 2 λ 2 |R| 2 2 + ν 2 + ν -trace(B -1 ) Σ -1 β, β 2 ,
for all λ > 0. This function reaches a global minimum at

λ * = ν + ν -trace(B -1 ) Σ -1 β, β 2 |R| 2 2 . ( 26
)
In practice, we set the penalization weight λ to λ * . According to Equation ( 26), this implies that the penalization hardens as the condition number of B increases. Such an increase occurs when the Hölder index H gets close to 1 (see Section 4).

Implementation

In this section, we describe an implementation of the estimation procedure in two dimensions (d = 2). The space W is dened as follows. Let p > 0. We set λ 0 = 1 and ϕ 0 ≡ 1, and for any m ∈ N * , ν 2m-1 = ν 2m = (1 + m p ), ϕ 2m-1 (θ) = cos(2mθ) and ϕ 2m (θ) = sin(2mθ). We dene a Sobolev space

W =    τ = m≥0 τ m ϕ m , m≥0 ν m τ 2 m < +∞    , (27)
and equip it with the inner product τ, τ W := m≥0 ν m τ m τm and its associ- ated norm |τ | W := τ, τ W . To compute the solution ( 22), we need to evaluate terms of L H . For that, we use a closed-form of Γ H .

Proposition 3 Let v be a mono-directional increment lter, i.e. v is of the

form v[l 1 , l 2 ] = v 1 [l 1 ], for (l 1 , l 2 ) ∈ [[0, L]] × {0} and 0 otherwise. Then, for θ ∈ [0, 2π), Γ H,v (θ) = c γ H (θ), (28)
where

γ H (θ) = | cos(ϕ)| 2H and c = 2 R+ |v 1 (ρ)| 2 ρ -2H-1 dρ.
To apply this formula, we compute the discrete Fourier transform γ H of γ H , and set the term L H of the jth row and mth column as γH [ m 2 ]ϕ m (θ j ). Let us quote that we do not evaluate the constant c. Hence, the solution is obtained up to a constant. However, this is not a matter for our application since we only use the normalized topothesy.

To solve the inverse problem, estimates H and β of H and β are also required. In our implementation, we took ordinary least square estimates of parameters of the linear model ( 14) [START_REF] Richard | Analysis of anisotropic Brownian textures and application to lesion detection in mammograms[END_REF](Richard, , 2016a,b),b). Using the ∆method and Theorem 1, it could be shown that these estimates are both unbiased and asymptotically Gaussian. We further replace Σ -1 by the inverse of a covariance matrix estimate of β.

Numerical Study

We evaluated our estimation procedure using 10000 realizations on a grid of size 800 × 800 of anisotropic fractional Brownian elds (see denition in Example 1). These realizations were simulated using the turning-band method developed by [START_REF] Biermé | A turning-band method for the simulation of anisotropic fractional Brownian eld[END_REF]. The Hurst function η of each simulated eld was set to a constant, which was sampled from a uniform distribution on (0.05, 0.95). Its topothesy was dened in the approximation space W ā for ā = 47. Its expansion coecients τ m were sampled from independent centered Gaussian distributions of decreasing variances

1 (1+ m 2 2
) . We set the coecient τ 0 = a m=1 |τ m | so as to ensure that the topothesy was non-negative. On each simulated eld, we computed increments and quadratic variations (see Equations ( 11) and ( 13)) at scales |u| and in directions arg(u) prescribed by some vectors u selected in the set {v ∈ N × Z, |v| ∈ [1, 20]}. For each represented direction, we took the three vectors with smallest scales if they existed. The set of transforms is represented on Figure 2. To compute increments, we used a kernel v of the form (12) with L = 2. Next, using these quadratic variations, we could estimate the Hölder index H and the intercept function β H,τ * of Equation ( 15) at 96 dierent angles. The root mean square error for the estimation of H was about 0.5%. Eventually, we computed several estimates of the asymptotic topothesy by solving Equation ( 22) in approximation spaces W a of dierent dimensions a ∈ [[1, ā]]; we used Sobolev spaces dened by Equation ( 27) with dierent values of p. In Equation ( 22), the parameter λ was set to the optimal value λ * given by Equation ( 26). For comparison, we also computed solutions obtained without penalization for λ = 0.

For each type of estimation, we evaluated the mean square error (MSE) by averaging squares of the quadratic distances between the estimated and true topothesy. We also evaluated the part of errors due to the approximation by averaging square distances between the original topothesy in the space W ā and its projection into an approximation space W a of lower dimension. These errors are plotted on Figure 3; they are expressed in percent of the mean square norm of true topothesy function. We rst comment Figure 3 (a) and (b). As the dimension of the approximation space was increased, the estimation error increased for both methods (with and without penalization). At lowest dimensions, the increase of these errors was compensated for by a decrease of approximation errors, leading to a decrease of MSE. Above a critical dimension, estimation errors became predominant and MSE started to increase. This critical dimension was only a = 6 without penalization and a = 44 with penalization. At this dimension, the minimum reached by the MSE was higher for the method without penalization (1.27%) than for the penalization one (0.82%). After this dimension, the MSE quickly went above 20% without penalization while it remained below 5% until the highest dimension with penalization. As a conclusion, without penalization, it was not possible to estimate correctly coecients of high-frequency components of the topothesy. Such an estimation could be accurately achieved with a penalization, showing the benet of the proposed method. However, as shown in Figure 3 (c), changing penalizations had an eect on the estimation. Taking a Sobolev penalization with p = 1 slightly reduced the performance, especially at lowest dimensions, probably due to a lack of penalization. Taking a penalization with p = 3 made the method failed at highest dimensions. For a same value of λ, the estimation bias induced by the penalization is higher for p = 3 than for p = 2. Hence, the optimal parameter λ has to be set at lower values for p = 3 than for p = 2. As a result, the penalization can be decreased when increasing p. To further analyse results of the penalization method (with p = 2), we computed MSE with an approximation dimension a = 44 by range of Hölder index values H; see Table 1. The Hölder index had an eect on the estimation performance. When it was close to 1, the MSE was large. In this case, the condition number κ of the matrix L T H Σ -1 L H was large, leading to a strong penalization and a high estimation bias (see Theorem 2 and Corollary 1). However, when H was not too large (< 0.8) the MSE was moderate (< 1%). Hence, in such a situation, we could obtain good estimates of high frequency components of the topothesy.

An application photographic papers

In this section, we present an application to historical photographic prints. The texture of these prints is a feature which is critical for the works of artists, manufacturers, and conservators [START_REF] Johnson | Pursuing automated classication of historic photographic papers from raking light images[END_REF][START_REF] Messier | Automated surface texture classication of inkjet and photographic media[END_REF]. In particular, conservators rely upon the texture to investigate the origin of an unknown print [START_REF] Johnson | Pursuing automated classication of historic photographic papers from raking light images[END_REF]. At present, such investigations are manually done by comparing the texture of the unknow print to those of identied references. They could be eased by an automated classication tool that would select relevant references and measure texture similarities between prints. To test the feasibility of automated classications of historic photographic paper, [START_REF] Johnson | Pursuing automated classication of historic photographic papers from raking light images[END_REF]; [START_REF] Messier | Automated surface texture classication of inkjet and photographic media[END_REF] assembled two datasets. The rst one, named the inkjet dataset, gathers 120 photomicrographs of non-printed inkjet papers collected from the Whilhelm Analog and Digital Color Print Material Reference Collection [START_REF] Messier | Automated surface texture classication of inkjet and photographic media[END_REF]. The second one, named the b&w dataset, is composed of 120 photomicrographs of non-printed silver gelatin photographic papers [START_REF] Johnson | Pursuing automated classication of historic photographic papers from raking light images[END_REF]. These datasets are publicly available at www.PaperTextureID.org. Some classication attempts were reported in conference papers by Abry et al. (2015b); [START_REF] Roux | Multiscale anisotropic texture unsupervised clustering for photographic paper[END_REF].

They are both organized into 4 groups containing sample sets of an increasing heterogeneity. In group 1, there are 3 sets of 10 samples obtained from a same sheet and expected by experts to have a high degree of similarity. In group 2, there are 3 sets of 10 samples from dierent sheets of a same manufacturer package which should also show a strong similarity albeit to a lesser extent. In group 3, there are 3 sets of 10 samples made from same manufacturer specications over time and expected to be more dissimilar. The group 4 is composed of 30 samples selected to show the diversity of papers. Datasets were documented by [START_REF] Johnson | Pursuing automated classication of historic photographic papers from raking light images[END_REF]; [START_REF] Messier | Automated surface texture classication of inkjet and photographic media[END_REF] with metadata including manufacturer, brand, date, type of texture and reectance. Photomicrographs of these datasets were acquired using a microscope system under the illumination of a single light placed at a 25 degree raking angle to the surface of the paper [START_REF] Johnson | Pursuing automated classication of historic photographic papers from raking light images[END_REF][START_REF] Messier | Automated surface texture classication of inkjet and photographic media[END_REF]. This specic illumination produces highlights and shadows reecting reliefs of the paper surfaces (see examples in Figures 4 and5). Depending on the paper properties, image textures look more or less rough and anisotropic. So as to characterize these paper properties, we used two texture features derived from our estimation procedure: the Hölder index H and an anisotropy index I dened as

I = [0,π) τ * (s) - [0,π) τ * (u)du 2 ds. (29)
These indices were computed by replacing H and τ * by their estimate (see Section 3.2). Let us notice that the estimated topothesy was normalized. Consequently, the anisotropy index was invariant to both the eld variance and the increment lter used for the estimation.

For the estimation of the Hölder index and the topothesy, we computed quadratic variations using a mono-directional increment lter of the form (12). The length L of the lter was adapted to each image so that increments satisfy conditions of Theorem 1. It was set to L = M + 2 (lter of order M + 1) Richard using an estimate M of the order of the IRF underlying the image. This order M was taken as the lowest one for which quadratic variations of image increments became almost constant at large scales. For most of the images, the maximal scale at which increments were computed was set to 20 pixels. But, for some papers (e.g. glossy inkjet papers), the relationship between logarithms of quadratic variations and scales was poorly linear at scales above 7 pixels. So, for these papers, the maximal scale was automatically set to 7 pixels. Between minimal and maximal scales, we used all possible increments in directions where at least two increments could be computed. Using quadratic variations of these increments, we estimated parameters of the linear model ( 14), including the Hölder index and the intercepts β. Using scales below 20 pixels (resp. 7 pixels), we could estimate intercepts in 96 (resp. 17) directions. Using these estimates, we set the procedure for the estimation of the topothesy function. In this procedure, the parameter a was set to 47 (resp. 8) so to ensure that it was about the half of the number of intercepts. We took a Sobolev penalization with p = 2. On Figures 6 and7, we ploted couples ( H, Ĩ) of textures of groups 1 to 3 for inkjet and b&w datasets. For textures of group 1 ("same sheet") and 2 ("same package") of the inkjet dataset, indices were both homogeneous within each set and separated across sets, showing a stability of the manufacturing process. Sets 61-70 and 71-80 of group 3 ("same manufacturer") were also quite homogeneous, but not the set 81-90. The variability of this set could be due to the variety of paper brands. Besides, we observed three point clusters corresponding to papers with a same reectance: A rst cluster composed of glossy papers with low Hölder and anisotropy indices, a second one containing semi-glossy papers with larger Hölder and anisotropy indices, and a third one formed by matte papers with larger Hölder indices. Some samples of these clusters are shown on Figure 4 with their associated indices. Let us note that the anisotropy index of glossy papers could be underestimated due to the use of a reduced number of topothesy coecients in the estimation procedure. Textures of the b&w dataset showed a larger intra-set variability than those of the inkjet dataset, even for papers of groups 1 and 2. Textures of these photographic papers might be less homogeneous than the inkjet ones. The variability was particularly large on the set 81-90 of group 3. This was probably due to the variety of paper reectance within this set. On the b&w semi-glossy papers (sets 31-60), which is considered as good by experts, was underestimated at level fair or poor. The anity between glossy papers (sets 1-10 and 61-90) and matte papers (set 21-30) was correctly estimated at level very poor whereas the one between semi-glossy papers (sets 61-90) and matte papers was slight overestimated at level poor. The anity between glossy papers of dierent brands (Canon: 1-10, HP: 61:70, Epson: 71-80, Kodak: 81-90) was mostly estimated at levels excellent or perfect whereas it is only qualied as good by experts. Globally, the anity matrices obtained by experts and the automated classication matched at a level of 40%.

The anity between b&w papers of a same group was well-estimated at level good or excellent, except for the group 81-90 which mixes glossy and lustre papers. The anity evaluated as poor by experts were globally wellestimated. This is for instance the case between glossy papers of Ilford group 41-50 and lustre papers from Kodak groups 11-20 and 31-40. Some anities qualied as good by experts were underestimated at levels poor or very poor between the lustre and glossy papers of the group 81-90 and groups 11-20, 41-50. For b&w papers, the overall match between anity matrices obtained by experts and the automated classication was of 60%.

Richard

We also computed the anity matrices using only the Hölder index. The degrees of match with the expert ones were 39% and 52% for inkjet and b&w papers, respectively. The benet of the anisotropic index was thus larger for the classication of b&w papers than inkjet papers.

Discussion

We presented an approach for the characterization of the anisotropy of Hölder random elds. Using the eld spectral density, we rst dened an asymptotic notion of topothesy which quanties the directional contributions of highfrequencies to the eld irregularity. We then designed and evaluated a procedure based on multi-oriented quadratic variations for the estimation of this asymptotic topothesy. Eventually, we used this procedure for the classication of textures of photographic papers. This classication was done by combining two features, a usual estimate of the Hölder index of the eld and a new anisotropy index derived from the estimated asymptotic topothesy. It led to some clusters gathering papers with similar reectances.

The anisotropy index we used are related to the ones proposed by [START_REF] Richard | Analysis of anisotropic Brownian textures and application to lesion detection in mammograms[END_REF]Richard ( , 2016a)). These indices are also obtained from the directional intercept function β H,τ * that appears in the linear asymptotic relationship between multi-oriented quadratic variations and the Hölder index of the eld (see Equation (15) of Theorem 1). They measure a dispersion of the intercept function and, indirectly, of the asymptotic topothesy. However, they depend on the Hölder index of the eld and the order of increments required for its analysis (this order has to be set with respect to the order of the eld). Hence, their variations do not exclusively reect dierences between eld directional properties. They also account for changes of eld regularity or order. By contrast, the anisotropy index of this paper is a direct measure of dispersion of the asymptotic topothesy. This measure is intrinsically related to directional properties of the eld and invariant to its order and Hölder index.

For the classication, information carried by the asymptotic topothesy was reduced to a single anisotropy index, but it is much richer. Depending on the interest for the application, other indices such as kurtosis or skewness of the asymptotic topothesy could be computed. Directions where the asymptotic topothesy reaches optima could also be of interest for nding main orientations of an image texture. Eventually, some dimensionality reduction algorithms could be directely applied to topothesy coecients to extract the most relevant information to be used for classication.

Besides, according to Proposition 2, the asymptotic topothesy provides us with an information about the correlation structure of eld high-frequencies. Such an information could be used to set an IRF model where the image would be decomposed into a trend eld and a texture eld with a specied correlation structure. Using such a model, it would become possible to achieve other image processing tasks such as separation of trend and texture of an image, examplar-based texture simulations, or inpainting of missing or occluded parts of an image.

Our characterization and estimation approach is only devoted to the global analysis of homogeneous random elds whose Hurst and topothesy functions remain spatially constant. It would be challenging to extend it to the local analysis of heterogeneous anisotropic elds. Such an extension is analogous to the one that passed from a global regularity analysis of fractional Brownian motions to a local one of multifractional Brownian motions (see for instance [START_REF] Istas | Quadratic variations and estimation of the local Holder index of a Gaussian process[END_REF] To achieve it, it would be necessary to set a suitable framework of heterogeneous anisotropic elds, which can not be intrinsic anymore. It would also imply to reduce the analysis area to small neighborhoods of points, which can be critical for the estimation.

The anisotropy characterization we proposed is valid for any eld of the framework, no matter how it is observed. However, our estimation procedure is specically designed for elds observed on a lattice. Due to this particular structure, observations can be analyzed without any interpolation using lattice-preserving transforms. An extension of this procedure to partial observations could be envisaged. However, it would probably require an interpolation of data, which could have an eect on the estimation.

A Proofs

Proof (Proposition 2) For some positive constant c, we have

∆ A (x) = |K A,f (x) -K A,f * (x)| ≤ c |w|≥A |f (w) -f * (w)|dw.
Since the eld density satises Condition (6), we further obtain Then, let us decompose the integral over E c 0 into the sum of two integrals, one over a set F δ = {s, η(s) -H > δ/2} dened for δ > 0, and the other over a set E δ = {s, 0 < η(s) -H < δ/2}. It follows that where µ(F δ ) is the measure of F δ on the unit sphere of R d . But, as shown by Richard (2016b), lim δ→0 µ(F δ ) = 0. Hence, letting δ = log(A) 1-α for some 0 < α < 1, we obtain Consequently, ∆ A (x) = o(A -2H ) as A tends to +∞.

Proof (Proposition 3) When the increment eld is mono-directional, the expression of Γ H,v of Equation ( 16) reduces to

Γ H,v (θ) = R + |v 1 (ρ cos(θ))| 2 ρ -2H-1 dρ.
The expression of β H,τ * in Equation ( 28) follows from the simple coordinate change u = ρ cos(θ).

Proof (Theorem 2) We aim at estimating the solution τ * of a linear system L H τ = β with a random vector

τ * λ = (B + λR) -1 L T H Σ -1 β
where B = L T H Σ -1 L H , R is a diagonal matrix, λ > 0, and β is an unbiased estimate of β of variance V (β).

Since β is unbiased, the expectation of τ * λ satises

(B + λR)E(τ * λ ) = L T H Σ -1 β.
Moreover, Bτ = L T H Σ -1 β. Hence,

(B + λR)(E(τ * λ ) -τ * ) = -λRτ * .
Thus,

|E(τ * λ ) -τ * | 2 = |(B + λR) -1 λRτ * | 2 ,
where | • | 2 denotes the 2-norm. Hence, using norm properties, we get

|E(τ * λ ) -τ * | 2 ≤ λ|(B + λR) -1 | 2 |R| 2 |τ * | 2 .
Therefore,

|E(τ * λ ) -τ * | 2 |τ * | 2 ≤ λ|B -1 | 2 |(I + λB -1 R) -1 | 2 |R| 2 ≤ λ|R| 2 ν - |(I + λB -1 R) -1 | 2 , ( 30
)
where ν -is the lowest eigenvalue of B.

Next, we establish a bound for |(I + λB -1 R) -1 | 2 . For any vector u, we have

|B -1 Ru| 2 ≥ 1 ν 2 + |Ru| 2 ≥ 1 ν + 2 |u| 2 .
Therefore, the lowest eigenvalue of B -1 R is above ν -1 + . Thus, the one of I + λB -1 R is above 1 + λν -1 + . Consequently,

|(I + λA -1 R) -1 | 2 ≤ ν + (ν + + λ) . (31)
Using Equation ( 30), we eventually obtain the inequality (24). Now, we turn to the variance of the estimator. We have

V (τ * λ ) = E (τ * λ -E(τ * λ ))(τ * λ -E(τ * λ )) T ) = (B + λR) -1 B(B + λR) -1 .
Hence, trace(V (τ * λ )) = trace (I + λB -1 R) -2 B -1 .

Moreover, since B -1 is a covariance matrix, any term B -1 ij of B -1 is bounded by B -1

ii B -1 jj . Hence, it follows that

trace(V (τ * λ )) ≤ |(I + λB -1 R) -1 ∆| 2 2 ≤ |(I + λB -1 R) -1 L T H | 2 2 |∆| 2 2 ,
where ∆ is a vector formed by terms (B) -1 ii . Noticing that |∆| 2 2 = trace(B -1 ) and using Equation (31), we get trace(V (τ * λ )) = ν + ν + + λ trace(B -1 ).

(32)

Besides, since Bτ * = L T H Σ -1 β, we have

|B| 2 |τ * | 2 ≥ |B -1/2 Σ -1/2 β| 2 ≥ √ ν - Σ -1 β, β 2 .
Hence,

1 |τ * | 2 ≤ ν + √ ν - Σ -1 β, β 2 .
Combined with Equation (32), Equation (25) follows.

Proof (Corollary 1) Using expressions of bias and variance in Theorem 2, we clearly see that the function g bounds the relative mean square error. Then, a simple variation analysis of this function suces to show that it reaches a global minimum at λ * .

Fig. 1 :

 1 Fig.1: Realizations of anisotropic fractional Brownian elds simulated using the turningband method of[START_REF] Biermé | A turning-band method for the simulation of anisotropic fractional Brownian eld[END_REF]. Fields of rst and second rows have Hölder indices of 0.3 and 0.6, respectively. Fields of a same column have a topothesy function which is represented on the third row. Hurst functions of elds are constant and equal to the Hölder index. All simulations were obtained using a same pseudo-random number sequence so as to highlight texture dissimilarities due to variations of the simulated elds.

Fig. 2 :

 2 Fig. 2: Visualization of the transforms used for computing quadratic variations. Each point corresponds to a transform. The norm and argument of a point give the rescaling factor and the rotation angle of the transform, respectively. The rotation angles are further represented by segment orientations.

Fig. 3 :

 3 Fig. 3: Estimation errors: (a) without penalization, (b) with a penalization from Sobolev space with p = 2, (c) with dierent penalizations.

Fig. 4 :

 4 Fig. 4: Patches of size 500 × 500 extracted from the inkjet dataset, associated to their metadata (sample number, manufacturer, reectance) and computed indices H and Ĩ.

Fig. 5 :

 5 Fig. 5: Patches of size 500 × 500 extracted from papers of the b&w dataset, associated to their metadata (sample number, manufacturer, reectance) and computed indices H and Ĩ.

Fig. 6 :

 6 Fig. 6: Hölder and anisotropy indices of texture samples of dierent groups from the inkjet dataset.

Fig. 8 :

 8 Fig. 8: Anity matrices from experts (left) and computations (right) for the inkjet dataset (top) and b&w dataset (bottom).

∆

  A (x) ≤ c |w|≥A |f (w) -gτ,η(w)|dw + |w|≥A |f * (w) -gτ,η(w)|dw, ≤ c |w|≥A |w| -2H-d-γ dw + |w|≥A |f * (w) -gτ,η(w)|dw, ≤ o(A -2H ) + |w|≥A |f * (w) -gτ,η(w)|dw, as A tends to +∞. Now, in directions s of the set E 0 = {s, τ * (s) > 0}, we notice that η(s) = H and τ * (s) = τ (s). Hence, f * (w) = gτ,η(w) whenever arg(w) is in E 0 , f * (w) = gτ,η(w), and 0 otherwise. Consequently, in polar coordinate, we have |w|≥A |f * (w) -gτ,η(w)|dw = E c 0 +∞ A gτ,η(ρs)ρ d-1 dρ ds, (θ)-1 dρ ds.

  |w|≥A |f * (w) -gτ,η(w)|dw = O(A -2H )(A -δ + µ(F δ )).

  |w|≥A |f * (w) -gτ,η(w)|dw = o(A -2H ).

Table 1 :

 1 Errors as a function of the Hölder index.

	Richard

dataset, we also observed three main clusters corresponding to groups of papers with common reectance. A rst cluster included the semi-matte and glossy papers (low Hölder and anisotropy indices), a second one the matte and glossy papers (low Hölder index and large anisotropy index), and a third one the lustre, champois and half-matte papers (large Hölder and anisotropy indices).

Papers from these dierent clusters are compared on Figure 5.

For each sample pair, we further determined a degree of anity ranging in ve levels from "very poor" to "perfect". This level of anity was computed by thresholding the euclidean distance between index values ( H, Ĩ) of samples. In this procedure, thresholds were automatically set to optimize over a whole dataset the matching between the computed anities and the ones established by an expert from meta-data. In Figure 8, we display in the form of a matrix the sample anities computed for each dataset and compare them to the expert ones.

The anity between inkjet papers of a same group was globally well estimated at level excellent, except for the group 81-90 (papers of the manufacturer Kodak). The level of anity between glossy papers (sets 1-10 and 61-90) and


Pascale Royer

Pascal Swider

Pauline Assemat

Homogenization of advection-diffusion and solid diffusion in poroelastic media for modelling transport of soluble factors in biological tissues

Homogenization of advection-diffusion and solid diffusion in poroelastic media for modelling transport of soluble factors in biological tissues

Pascale Royer, Pascal Swider and Pauline Assemat

This work has been motivated by the study of solid tumor growth and more specifically by the examination of the most frequent pediatric primary bone tumors (osteosarcoma), which are characterised in their early stage by the formation of non-mineralised bone tissue, called osteoid [1]. As the tumor evolves in time, mineralisation of this growing tissue can take place. The purpose of this work thus is to investigate the advection-diffusion of soluble factor such as calcium and phosphate [2] in the context of tumor mineralisation, how this mineralisation of bone tissue may become a barrier to treatment in the context of drug transport [3] but also how proteins and growth factors can play a role in tumor growth [4]. For this purpose, the method of asymptotic homogenization is used so as to derive the macroscopic models for describing solute transport in poroelastic media, with an advective-diffusive regime in the fluid-saturated pores and diffusion in the solid phase, by starting from the description on the pore scale. The fluid/solid equations lead to Biot's model of poroelasticity. Then, homogenization of the transport equations leads to three macroscopic models that relate to three orders of magnitude of the diffusivity ratio: a model in which the solid diffusion only influences the accumulation term; ii) a model with memory effects; iii) a model without solid diffusion. Initially expressed by means of orders of magnitude of the diffusivity ratio, the domains of validity of each of these three models can be expressed in terms of relative orders of magnitude of two characteristic times. The three models contain a solute-solid interaction term, due to the advection regime and they are coupled to the poroelaticity model via the advection term.

[1] Fletcher, C. D.M., Bridge, J.A., Hogendoorn, P., Mertens, F. ( 2013

  ): WHO Classification of Tumours of Soft Tissue and Bone, Volume 5. [2] Giacaman, R., Perez, V., Carrera, C. (2016): Mineralization processes in hard tissues, Biomineralization and Biomaterials, 147-185. [3] Sciumè, G., Shelton, S.E., Gray, W.G.T, Millers, C.,Hussain, F., Ferrari, M., Decuzzi, P., Schrefler, B. (2012): Tumor growth modelling from the perspective of porous media mechanics, Molecular and cellular biomechanics : MCB. 9. 193-212. [4] Ambard, D., Swider, P. (2006): A predictive mechano-biological model of the bone-implant healing, European Journal of Mechanics -A/Solids. 25. 927-937.


Jiyeon Choi

Matthieu Bellec

Arnaud Royon

Kevin Bourhis

Gautier Papon

Thierry Cardinal

Lionel Canioni

Martin Richardson

Three-dimensional direct femtosecond laser writing of second-order nonlinearities in glass

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Three-dimensional direct femtosecond laser writing of second-order nonlinearities in glass Jiyeon Choi, 1,2 Matthieu Bellec, 2, * Arnaud Royon, 2 Kevin Bourhis, 3 Gautier Papon, 2 Thierry Cardinal, 3 Lionel Canioni, 2 and Martin Richardson 1 Second-order nonlinear optical phenomena such as second-harmonic generation (SHG) and the Pockels effect are widely used in many photonic devices, enabling new laser source wavelengths and the high speed external modulation of light for multiple purposes [START_REF] Saleh | Fundamentals of Photonics[END_REF]. These devices are, however, inherently restricted to using noncentrosymmetric materials such as nonlinear optical crystals or poled materials that exhibit second-order nonlinearities [START_REF] Boyd | Nonlinear Optics[END_REF]. Despite the existence of a number of nonlinear optical crystals and the improving capabilities of poled materials, there are still essential gaps in the range and operability of these devices. Commercial applications relying on the second-order nonlinearity are still limited due to this narrow choice of suitable materials.

There is thus a lingering need for a new class of nonlinear optical materials possessing second-order nonlinearities sufficiently large over a wide spectral range.

Glasses have been widely used as optical materials for many applications, in part because of the ease in tailoring their properties by changing their composition. Their remarkable properties, such as good mechanical and chemical durability, allow ease of preparation in various forms, such as fibers and curved substrates [START_REF] Zarzycki | Glasses and the Vitreous State[END_REF]. The possibility to incorporate in these vitreous materials noble metal ions allows one to increase their photosensitivity and to access to new or improved optical properties. If one can induce second-order nonlinearity by breaking their centrosymmetry, then perhaps glasses would become promising candidates for second-order nonlinear optical applications.

There are various techniques to break the centrosymmetry. One way is to polarize the glass, the so-called poling process, either optically in optical fibers [START_REF] Österberg | [END_REF] or thermally [5]. There have been many attempts to produce thermally poled glass systems in silicate [6], tellurite [7], and phosphate [8] matrices. The basic idea of thermal poling is to create a frozen electric field E dc inside a non-linear layer as an anodic depleted-ion region caused by the localized migration of mobile ions. The interaction between this permanent electric field and the third-order nonlinear susceptibility χ 3 of the glass gives rise to an effective second-order nonlinear response χ 2 eff 3χ 3 E dc , the well-known electric field-induced secondharmonic generation (EFISHG) effect.

Femtosecond laser direct writing is a technique that can be utilized in transparent materials to fabricate microstructures in three dimensions with modified physical properties thanks to multiphoton absorption [9]. This technique can be applied to induce localized secondorder nonlinearities in glass, either by changing the crystalline phase from amorphous into crystallite [10,11], or by modifying the morphology of metal nanoparticles embedded in glass from spherical to ellipsoidal [12].

In this Letter, we propose another approach to create inside a silver-containing zinc phosphate glass a laserinduced frozen charge gradient, similar to the thermal poling process, but localized with three-dimensionality. We show that the inhomogeneous charge distribution produced in the laser-exposed glass can give rise to a second-order nonlinearity that can be exploited for nonlinear optical applications. We suggest a model for this process and demonstrate SHG from a silver microstructure to prove the concept. The polarization of the SHG is analyzed to confirm the theoretical EFISHG model. Finally, the SHG efficiency of a silver microstructure is compared to quartz to estimate the induced χ 2 eff . The glass used for this study is a photosensitive zinc and gallium phosphate glass (40P 2 O 5 -55ZnO -4Ga 2 O 3 -1Ag 2 O) containing 1% (mol.%) of silver oxyde. Details of the glass fabrication process and the material properties are described elsewhere [13]. A Yb-doped fiber 450 fs ultrafast laser (μJewel, IMRA America Inc.) operating at 1045 nm with a repetition rate of 1.5 MHz and a maximum pulse energy of 200 nJ was used for the irradiation.

A microscope objective (NA 0.25) focused the laser beam into the glass sample mounted on a computer controlled three-dimensional (3D) stage. Silver microstructures consisting of fluorescent silver clusters Ag x m m < 10 were created by irradiating the samples for 1 second at an irradiance of 5.6 TW:cm -2 in a 4 μm diameter spot. The key mechanisms contributing to this clustering is a combination of photo-reduction (2Ag → Ag 0 Ag 2 ), photo-dissociation, and thermal and chemical diffusions, due to heat accumulation resulting from the high repetition rate laser pulses [14,15]. As a result, these fluorescent clusters are arranged into a ring shape in 2D (see Fig. 1(a)), a pipe shape in 3D. Figure 1(b) depicts schematically and hypothetically the laser-induced distribution of the silver species. Permanent static electric fields with circular symmetry are formed between the hole-trap color centers, such as Ag 2 located inside the rim and the electron-trap silver clusters Ag x m situated on the outside of the rim of the microstructure [15]. Although the spatial distribution of the silver species is clearly radial, their exact repartition inside the ring in not yet known and needs to be measured by near-field experiments.

Treating this polled structure as a nonlinear medium irradiated now by a probe beam E probe of angular frequency ω, the induced second-order polarization in the structure depicted in Fig. 1(a) can be represented as decomposed forms in x and y [START_REF] Boyd | Nonlinear Optics[END_REF]:

P x 2ω 3ε 0 χ 3 xxxx E 2 probe x E dc cos θ; P y 2ω 3ε 0 χ 3 yxxy E 2 probe x E dc sin θ; ( 1
)
where θ is the angle between an electric field vector and the polarization of the probe beam. P x and P y represent the polarizations along the x (horizontal, H) and y (vertical, V) axes. The probe beam E probe x is polarized along the x-direction. The second-harmonic (SH) irradiances associated to the above second-order polarizations are given by I HH 2ω ∝ jχ 3 xxxx j 2 jE 2 probe x j 2 jE dc j 2 j cos θj 2 ; I HV 2ω ∝ jχ 3 yxxy j 2 jE 2 probe x j 2 jE dc j 2 j sin θj 2 .

(

) 2
The SHG measurements were performed with the same laser (see Fig. Figure 3 shows the experimental and theoretical derived images of the SH signal for parallel (HH in Fig. 3 versus the angle of the electric field vector for parallel (HH) and crossed (HV) polarizers and analyzers. As expected, the ratio of the intensity maxima for each polarization is I max HH ∶I max HV 9∶1 for horizontal and vertical directions, as given by the Kleinmann symmetry relation χ 3

xxxx ∶3χ 3 yxxy 3∶1 in glass [START_REF] Boyd | Nonlinear Optics[END_REF]. The theory is in good agreement with the experimental images, except for the central region where a signal is observed. This signal seems to be insensitive to the analyzer orientation and is quite directional. The physical origin of this signal is not clear and is still under investigation. It could be due to a pump leakage, a contribution of the longitudinal component of the focused beam, a guiding effect that depolarizes the SH signal, or the 3D organization of the structure (pipe), which is not taken into account in our 2D model (ring).

Finally, in order to estimate the SHG efficiency of this microstructure, it was compared with that of z-cut quartz and normalized to the SH intensity of the quartz crystal exposed to the same probe irradiance (Fig. 4). The estimated χ 2 eff is 2.5 that of quartz, corresponding to ∼1.2 pm:V -1 . Thus the laser-induced value of χ 2 is comparable to that achieved with thermal poling of glass of similar composition, i.e., ∼1.8 pm:V -1 [16]. The thirdorder susceptibility of the silver-containing zinc phosphate glass is of the same order of magnitude as that of fused silica, i.e., χ 3 2.5 × 10 -22 m 2 :V -2 . The internal electric field can therefore be estimated to be ∼1.6 × 10 9 V:m -1 . This value is extremely large and is close to the electric field of an electron-ion dipole in an atom. Note that one possibility to get a higher χ 2 eff value could be to modify the glass composition and the introduction of hyperpolarizable ions in order to increase the third-order nonlinearity of the material.

In conclusion, we have proposed and demonstrated a new method based on femtosecond laser-induced depletion to create high second-order nonlinear optical properties in silver-containing photosensitive glass. This nonlinearity was produced in μm-scale structures, consisting of charged silver clusters, inducing an inhomogeneous space charged distribution with an associated permanent static electric field. We have estimated the induced second-order susceptibility of these structures to be ∼2.5 times that of quartz. This new laser-based method of creating new nonlinear optical materials opens a new pathway in the development of active and nonlinear optical components. Moreover, since this laser writing method allows for the micron-scale structuring of this nonlinearity in precision-designed 3D configurations, it opens the door to the development of many new classes of optical systems and components including multifunctional lab-on-chip devices, active 3D optical storage devices [13], and high damage threshold optical converters.

2

 2

  (a)), but with a power reduced by ∼57% of the power needed for structuring of the material. The SH signal from the structured sample was collected through a microscope objective (NA 0.4) and a set of filters. The SH signal was coupled into an intensified CCD (ICCD) camera attached to a spectrometer via a fiber-optic delivery. A CCD beam profiler was used to visualize the SHG distribution. A polarizer was located in front of the CCD beam profiler to analyze the polarization dependence of SH signal. The presence of SHG was confirmed by the spectral line centered at 522 nm for different probe irradiances (Fig.2(b)) and by the quadratic dependence of the SH intensity relative to the fundamental irradiance (Fig.2(c)).

  (a) and (c)) and crossed (HV in Fig. 3(b) and (d)) polarizers and analyzers. The experimental images show a central spot encircled by two meniscal lobes with a scale corresponding to what we believe to be the nonlinear region formed by the space charge distribution.

Figure 3 (Fig. 1 .Fig. 2 .

 312 Fig. 1. (Color online) (a) Confocal fluorescence image of a photo induced ring structure with an excitation at 405 nm. The scale bar is 5 μm. (b) Schematic representation of a possi ble radial distribution of the silver species inside the ring. The permanent electric field is E dc .

Fig. 3 .

 3 Fig. 3. (Color online) Experimental and associated theoretical images of the SH signal for parallel HH (a),(c) and crossed HV (b),(d) polarizers and analyzers. (e) Evolution of the SH inten sity versus the angle of the electric field vector for parallel (HH) and crossed (HV) polarizers and analyzers.

Fig. 4 .

 4 Fig. 4. (Color online) Experimental images of the SH signal for a parallel HH polarizer and analyzer of a (a) single ring struc ture and of (b) bulk quartz attenuated 5000 times. (c) Evolution of the normalized SH intensity of the ring structure (green full line) and of bulk quartz (green dotted line) versus the radial coordinate. The gray zones correspond to the integration areas where the χ 2 of the ring is measured.

1

  Townes Laser Institute, CREOL, University of Central Florida, Orlando, Florida 32816 2700, USA 2 LOMA CNRS/University of Bordeaux, 351 Cours de la Libération, 33405 Talence Cedex, France 3 ICMCB CNRS/University of Bordeaux, 87 Avenue Dr. Schweitzer, 33608 Pessac Cedex, France *Corresponding author: bellec@unice.fr We demonstrate that direct femtosecond laser writing in silver containing zinc and gallium phosphate glass enables generation of three dimensional (3D) optical second order nonlinear microstructures having an χ 2 value about 2.5 times that of quartz. The proposed physical model involves photo reduction, photo dissociation, and migration of silver species within the glass matrix. 3D laser written second order nonlinear structures could become a new class of nonlinear optical components. OCIS codes: 190.4400, 160.5335, 160.1245, 140.3390.

The authors acknowledge the Conseil Régional d'Aquitaine, the GIS AMA, and the ANR (grant 2010 BLAN 946 03) for their financial support, as well as P. Legros and the BIC for the confocal fluorescence microscopy images.