Laure Arbenz

Olivier Chadebec

Christophe Espanet

Youness Rtimi

Gilles Cauffet

Characterization of Permanent Magnet Magnetization

teaching and research institutions in France or abroad, or from public or private research centers.

I. INTRODUCTION

M

necessary to mesh the air [START_REF] Arbenz | Characterization of the local incremental permeability of a ferromagnetic plate based on a four needles technique[END_REF]. But if the magnetic circuit is open, both the air box and the mesh must be large. Moreover, this method is known to suffer from numerical noise when calculating magnetic fields around magnets. Another drawback is that the forward relation between magnetization and stray field in air is not explicit.

On the other hand, integral methods offer an explicit and direct link between the magnetization and the field in the surrounding air region. Moreover, this region does not have to be meshed, and integral approaches are well known to be very accurate if focused on the computation of stray fields. This is why it is suitable and efficient to use integral methods in the context of magnetostatic inverse problems.

A few papers have already been published on that topic. In [START_REF] Igarashi | Inverse inference of magnetization distribution in cylindrical permanent magnets[END_REF], an integral inverse formulation similar to the one used in this paper is proposed, but it is limited to 2-D cases. In [START_REF] Chadebec | Recent improvements for solving inverse magnetostatic problem applied to thin shells[END_REF] and [START_REF] Chadebec | A review of magnetostatic moment method[END_REF], the approach is improved, also covering ferromagnetic materials (with high permeability), but remains limited to magnetic sheets and thus 2-D magnetizations. In a more recent work, the determination of magnetizations of 3-D volumes is proposed. The forward relation is based on a coupling between FE and surface integral methods [START_REF] Bruckner | Solving large-scale inverse magnetostatic problems using the adjoint method[END_REF]. The uniqueness of the solution is ensured with an optimization process, coupled with a regularization function. In this paper, 3-D magnetization distribution is also identified, but the use of linear purely integral formulation is preferred. The uniqueness of the solution is then implemented by truncating the matrix spectrum.

The paper is structured as follows. First, the inverse problem formulation is defined; then, the experimental setup is described, and finally, the method is used for the characterization of different samples.

II. INVERSE PROBLEM FORMULATION

A. Forward Integral Formulation

Let us consider a ferromagnetic domain associated with an unknown distribution of magnetization. This magnetization creates an external flux density which can be expressed by the following integral equation [START_REF] Durand | Magnétostatique[END_REF]:

B(P) = - μ 0 4π ∇ M • r r 3 d (1)

Characterization of Permanent Magnet Magnetization

Laure Arbenz1 , Olivier Chadebec 2 , Christophe Espanet 1 , Youness Rtimi 1,2 , and Gilles Cauffet 2

The performances of electrical machines highly depend on the behavior of materials, such as steel, copper, and permanent magnet. This paper focuses on permanent magnets whose magnetization must be precisely known. Thus a non-destructive approach which includes a combination of experiments and inverse problem, using the magnetostatic moment method, is presented in order to determine the magnetization of the magnet. With a Hall sensor, a scan of the field around the sample is made. These measurements and a numerical approach allow us to deduce the sample magnetization. Studies of two samples illustrate the advantages of this method.

Index Terms-Magnetostatic moment method, measurement, permanent magnet.

ANY mechatronic systems include a permanent magnet whose magnetization must be precisely defined and known [START_REF] Loussert | An efficient and optimal moving magnet actuator for active vibration control[END_REF], [START_REF] Gilson | 2-D analytical subdomain model for high-speed permanent-magnet machines[END_REF]. For example, a fault on a magnet can drastically reduce the properties of an electric machine (through increased detent torque or iron losses…). The standard method to measure the magnetic moment of permanent magnets is the use of a Helmholtz coil in combination with a fluxmeter [START_REF]Magnetic Materials-Part 14: Methods of Measurement of the Magnetic Dipole Moment of a Ferromagnetic Material Specimen by the Withdrawal or Rotation Method[END_REF]. This approach only makes sense for a two-pole magnet, and the angle of the magnetization is quite complex to identify.

For a magnetization with more complex patterns than a mere bipolar orientation, an alternative approach must be applied. It is possible to measure the magnetic stray field of the specimen with a magnetic field sensor, typically a Hall effect probe. The magnetic moments and the magnetization angles can then be deduced from the measured signals by using an adapted algorithm. For example, this measurement process can be applied in-line during the manufacturing process to detect bipolar magnets with unknown orientation [START_REF] Meyer | Concept for magnet intra logistics and assembly supporting the improvement of running characteristics of permanent magnet synchronous motors[END_REF].

However, this task is not trivial and leads to the resolution of an inverse problem. Starting from the stray magnetic field measured around the magnet on several points (i.e., the effect), the algorithm has to reconstruct the local distribution of magnetization (i.e., its origin or sources). This kind of identification problem is known to be ill-posed, meaning that the solution is not stable versus the measurement noise: mathematically speaking, the solution is said to be "non-unique."

Before solving an inverse problem, the first step is to solve the forward one. In our case, Maxwell equations under magnetostatic assumption have to be considered. A key point is the choice of the numerical method. This choice will determine the numerical link between the magnetization and the generated stray magnetic field. Once this relation is determined, it has to be inverted.

With the finite element (FE), which is the most classical method used for static magnetic field calculation, it is where B is the flux density (T), P is a point located in the surrounding air region, M is the magnetization (A/m), and r is a vector linking the integration point in to point P. This expression has the advantage of explicitly linking the source which has to be determined (i.e., the magnetization) to its effect (i.e., the flux density in the surrounding air region).

B. Integral Computation

Let us mesh the domain into n volume elements. We assume that a uniform magnetization is associated with each of them. Integral (1) presents a (∼1/r 3 ) singularity. Computing it directly by numerical integration can lead to significant inaccuracies, especially if the flux density has to be expressed in close proximity to the magnetized material. An alternative is using the well-known equivalent magnetic charge distribution based on Coulombian approach [START_REF] Durand | Magnétostatique[END_REF]. The magnetization M is constant for each element, so its divergence vanishes such as no volume magnetic charge exists. Thus, the field is created only by surface charge distributions M•n, located on the surfaces delimiting volume elements, with vector n being the external normal of volume elements. The advantage of this approach is to present a reduced singularity (∼1/r 2 ) as compared to (1) (∼1/r 3 ). Thus, we get [START_REF] Chadebec | Recent improvements for solving inverse magnetostatic problem applied to thin shells[END_REF], [START_REF] Chadebec | A review of magnetostatic moment method[END_REF]

B(P) = - μ 0 4π n i n f i S j M i • n j r r 3 d S j (2)
where S j is the j th face of element i , n j is the external normal vector, and n f is the number of faces associated with one element of the mesh. The integral can be very accurately computed using well-known analytical expressions, even if r is very small (i.e., when sensors are located very close to the ferromagnetic material). Equation ( 2) can be interpreted as the linear relation linking the magnetization to the flux density measured by one sensor.

C. Linear Inverse Problem Formulation

Let us consider that m field probes are located in the "air" region surrounding the magnetic structure. In addition, a cartesian vector basis is associated with each volume element, leading to three degrees of freedom for each unknown magnetization. Thus, we get the following linear matrix system:

B mes = [A]M (3)
where B mes is the 3 ×m measurement vector, [A] a (3 × n, 3×m) matrix, and M represents 3×n vectors of the unknowns. This inverse problem is classically ill-posed [START_REF] Hansen | The use of the L-curve in the regularization of discretization of discrete ill-posed problem[END_REF], and its solution is not unique. In fact, even a very low noise on magnetic measurements can lead to huge variations of the solution. If so, the resolution process of the inverse problem needs to be stabilized. Many different techniques exist, but a good strategy to solve (3) is to use single value decomposition to compute the pseudoinverse of [A] and to truncate its spectrum in order to limit the noise amplification

M = [A] † B mes . ( 4
)
A discrete evaluation of the magnetization is then obtained for the defined mesh.

III. EXPERIMENTAL SETUP AND TEST PROCEDURE

To measure the magnetic stray field of a given specimen presented in Fig. 1(a), a three-axes Hall effect sensor is used (Senis, model F3A-03C02C-A.2T2K5J-E). The full scale of this sensor is ±0.2 T for a total measurement accuracy below ±0.1%. A sensitivity study shows that the relative position between this sensor and the magnet must be accurately known. For this purpose, a homemade precision positioning bench is used [Fig. 1(b)]. The Hall sensor is fixed on an aluminum fixture and a non-magnetic support made of Ertacetal (POM plastics) tailored to fit the magnet geometry is fixed to the positioning bench. These non-magnetic surroundings of the magnet minimize potential risks of disturbing the measurement. The control of the positioning bench and the acquisition of the data are completely automated, effectively mitigating operator influences on the measurement process and thus its results.

For samples with an annular (or circular segment) shape [Fig. 1(a)], the sensor is fixed in position, while the specimen is rotated. With this experimental protocol, it is possible to measure the field on equidistant curves (arcs) surrounding the sample (Fig. 2). This operation is repeated for 11 different Magnetization of the magnet deduced by the method (three measurements). curves around the sample, which allows us to properly map the surrounding field.

IV. RESULTS AND DISCUSSION

A. Deduced Magnetization

The magnetization of the sample can be deduced from the field measurement with a numerical method. The chosen method is the magnetostatic moments method (magnetostatic MoM) [START_REF] Chadebec | A review of magnetostatic moment method[END_REF], [START_REF] Harrington | Field Computation by Moment Methods[END_REF]. As explained in Section I, the MoM is more powerful than the classical FE because surrounding air does not require any meshing and especially since MoM can easily be inverted. This method is thus well adapted to identify magnetization patterns. In this paper, numerical problems are solved with a software package using the magnetic moment method and called LOCAPI. Results of the sample presented in Fig. 1(a) are shown in Fig. 3. We can conclude that the sample has a diametrical magnetization.

B. Uniqueness of Solution

It is important to determine whether the calculated magnetization is the unique solution of the problem. To check this, the non-observable space magnetization must be determined [START_REF] Penrose | A generalized inverse for matrices[END_REF].

The solution vector M is projected onto the null space of the matrix A using the adapted projection operator

M no = (I-[A] † A)M. ( 5
)
If M no is sufficiently low compared to M, the solution can be considered as unique. From a physical point of view, the component M no which satisfies A•M no = 0 is a magnetization which does not leave any traces on the sensors. Thus, M no is the non-returnable information using the inversion procedure.

It is interesting to note that A is only a function of the geometry of the problem (sample geometry and mesh, location of the measurements). Therefore, the only way to overcome this issue is changing the sensors' spatial positions or the sample mesh. Mathematically, this equals reducing the nonobservable space (i.e., null space of A). It is thus possible to determine whether the measurement definition and the mesh are properly adapted before even launching the test sequence.

C. Sensitivity Study

To illustrate the reproducibility of the method, the same sample has been characterized three times, using the same experimental and numerical methods. These three measurement runs are completely independent of each other. The magnetization results are given in Fig. 4.

The standard deviation is less than 0.04 T (0.06 T) for radial (tangential) remanent flux density on an angular aperture of ±60°(the angular aperture of the magnet is ±70°). So the reproducibility of the method is sufficient to clearly identify the magnetization of the magnet.

However, the standard deviation is five times larger at the tips of the magnet (-70°to -60°and +60°to +70°). This phenomenon is due to a difference between the physical magnet position during the measurements and the theoretical magnet position used in the numerical model. The magnetization at the end of the magnet is therefore either oversized or undersized.

V. APPLICATION OF MEASUREMENT METHOD

We then propose to study two samples A and B with the same physical shape [Fig. 1(a)] but different magnetization patterns. Both samples are what we call "three-pole magnets" (e.g., a N / S / N pattern), one of which has a radial magnetization and the other one a sinusoidal magnetization. The radial and tangential magnetizations are given in Fig. 5(a). Both samples are clearly different, but it is simpler to evaluate when comparing the flux density versus the angular position [Fig. 5(b)]. It is evident that sample A is radial (visible steps at the pole transitions) and sample B is sinusoidal (linear flux rotation).

VI. CONCLUSION

The proposed approach allows deducing 3-D magnetization distribution of a magnet from stray magnetic field measurements. It is based on the inversion of a pure volume integral formulation of the magnetic field regularized with a matrix spectrum truncation. Being a non-destructive approach, it involves a combination of experiments and inverse problem using the magnetostatic moment method. The procedure was applied to samples with identical geometry but different magnetization patterns. The results allow us to reliably distinguish samples without any ambiguity. A major advantage of the method is that it can be applied whatever the geometry of the sample. Further improvement may be possible by adding a probe positioning sensor on the test fixture. This should substantially increase the reproducibility of the method. Moreover, it may make sense to find a way to automate the localization of measurement points, depending on geometry and mesh of the sample to be investigated.

Fig. 1 .

 1 Fig. 1. (a) Shape sample. (b) Experimental setup of the proposed method.

Fig. 2 .

 2 Fig. 2. Measured field (magnetic flux density) versus angular position.

Fig. 3 .

 3 Fig. 3. Magnetization of the magnet deduced by the method. Overlay: visualization of magnetization vectors.

Fig. 4 .

 4 Fig. 4.Magnetization of the magnet deduced by the method (three measurements).

Fig. 5 .

 5 Fig. 5. (a) Radial and tangential magnetization of samples A and B. (b) Remanent flux density angle versus angular position.

Research and Development Department, Moving Magnet Technologies,

000 Besançon, France 2 G2Elab, CNRS, Grenoble INP, University of Grenoble Alpes, F-38000 Grenoble, France


domain: sdv spi chim phys

text:
Tijana Lainović

Jérémie Margueritat

Quentin Martinet

Xavier Dagany

Larisa Blažić

Dejan Pantelić

Aleksandar Mihailo D Rabasović

J Krmpot

Thomas Dehoux
email: thomas.dehoux@univ-lyon1.fr

Aleksandar J Krmpot

Micromechanical imaging of dentin with Brillouin microscopy

come

INTRODUCTION

Human dentin is an organized, hard, mineralized tissue of the tooth, composed of 70 wt% calcified tissue (hydroxyapatite), 20 wt% organic phase (mostly composed of collagen type I as well as other fibrils), and 10 wt% water [START_REF] Goldberg | Dentin: structure, composition and mineralization[END_REF]. It is perfused with microtubules that allow sensory communication to the underlying pulp and nerves. [START_REF] Stock | Bovine and equine peritubular and intertubular dentin[END_REF] This fiber-composite-like structure yields an anisotropic distribution of mechanical properties that support masticatory stress concentrated in the surrounding enamel. [START_REF] Kinney | Hardness and Young's modulus of human peritubular and intertubular dentine[END_REF] Tooth mechanical strength can be altered by diet, age or diseases such as caries or sclerosis. [START_REF] Zaytsev | Deformation behavior of root dentin under Sjogren's syndrome[END_REF] It is of utmost importance to characterize mechanical properties during tooth life in order to predict and understand tooth decay, to design restorative dental procedure, and to investigate the tribological behavior of teeth.

Most of the characterization means are inherited from material science, and involve destructive testing such as indentation or fatigue tests. [START_REF] Kinney | Hardness and Young's modulus of human peritubular and intertubular dentine[END_REF][START_REF] Harley | Phonons and the elastic moduli of collagen and muscle[END_REF] Due to the variety of techniques and the complexity of the dentin microstructure, there is a lack of quantified consensus on the mechanical properties of dentin. [START_REF] Kinney | The mechanical properties of human dentin: a critical review and re-evaluation of the dental litterature[END_REF] Moreover, such approaches are clearly difficult to put into practice in vivo. In the clinic, dentinal pathological changes can be detected by visuo-tactile and radiographic methods to guide clinicians in establishing a relevant diagnostic, and in proposing an adequate therapy. [START_REF] Makhijaa | Dentists' decision strategies for suspicious occlusal caries lesions in a National Dental PBRN study[END_REF] However, such approaches cannot give precise information on the microstructure, cannot delineate precisely caries margins, [START_REF] Schwendicke | Effects of using different criteria for caries removal: A systematic review and network meta-analysis[END_REF] and are highly depend on the clinician [START_REF] Joves | Nanoindentation hardness of intertubular dentin in sound , demineralized and natural caries-affected dentin[END_REF].

Acoustic techniques have long been used to probe the mechanics of dentin, initially analyzing the reflection or transmission of ultrasonic pulses emitted from a piezoelectric transducer in the MHz range. [START_REF] Kossoff | Examination of the contents of the pulp cavity in teeth[END_REF] Later implementations in tooth tissues in various healthy and pathological states used resonant ultrasonic spectroscopy, [START_REF] Kinney | The mechanical properties of human dentin: a critical review and re-evaluation of the dental litterature[END_REF] acoustic microscopy, [START_REF] Maev | New data on histology and physico-mechanical properties of human tooth tissue obtained with acoustic microscopy[END_REF][START_REF] Raum | Preservation of microelastic properties of dentin and tooth enamel in vitro-A scanning acoustic microscopy study[END_REF][START_REF] Shelke | Visualisation of localized elastic properties in human toth and jawbone as revealed by scanning acoustic microscopy[END_REF] surface waves [START_REF] Toda | Nondestructive testing in human teeth using a leaky Lamb wave device[END_REF] or laser ultrasonics. [START_REF] Sun | Experimental and numerical studies for nondestructive evaluation of human enamel using laser ultrasonic technique[END_REF], [START_REF] Wang | Laser ultrasonic evaluation of human dental enamel during remineralization treatment[END_REF] However, these techniques are limited to the MHz range, i.e. a ∼ 100 µm resolution, unable to reveal the micromechanics of the tooth, in particular at tissue interfaces.

State-of-the-art microscopies, such as multiphotonic or non-linear technologies, [START_REF] Cloitre | Multiphoton imaging of the dentine-enamel junction[END_REF] can offer optical resolution but cannot provide any information on the mechanics of the sample. Laser spectroscopies have recently rose to prominence, offering quantitative imaging of the chemical structure of tooth tissues and restorative materials. [START_REF] Salehi | Functional mapping of human sound and carious enamel and dentin with Raman spectroscopy[END_REF][START_REF] Ramakrishnaiah | Applications of Raman Spectroscopy in Dentistry: Analysis of Tooth Structure[END_REF][START_REF] Lopes | Fourier transform infrared spectroscopy (FTIR) application chemical characterization of enamel, dentin and bone[END_REF][START_REF] Steelman | Light-scattering methods for tissue diagnosis[END_REF] In this paper we propose a new approach to probe acoustic waves in the GHz range by laser spectroscopy, thus offering quantitative microscopy with a contrast based on the mechanical properties of tooth tissues.

We developed a microscope based on Brillouin light scattering (BLS) to probe mechanical changes in tooth tissues. BLS is an inelastic process similar to Raman scattering. Raman-based techniques rely on the scattering of light by optical phonons, and operate at THz frequencies. They provide information on the molecular structure of the sample. By contrast, Brillouin scattering uses the scattering of light by acoustic phonons. [START_REF] Brillouin | Diffusion de la lumière et des rayonnés X par un corps transparent homogène; influence de l'agitation thermique[END_REF] Since acoustic phonons have lower energies than optical phonons, Brillouin scattering produces frequency shifts in the GHz range that reveal the mechanical properties of the sample. Developed in the 60's for solid-state physics, BLS has only recently emerged as a key technology in life sciences due to its ability to produce non-contact, label-free microscopic images of the mechanical properties of cells and tissues. [START_REF] Palombo | Brillouin Light Scattering: Applications in Biomedical Sciences[END_REF] It has notably been applied for the analysis of biological fibers, connective tissues and muscles [START_REF] Harley | Phonons and the elastic moduli of collagen and muscle[END_REF][START_REF] Randall | Brillouin scattering in systems of biological significance[END_REF][START_REF] Palombo | Biomechanics of fibrous proteins of the extracellular matrix studied by Brillouin scattering[END_REF], eye tissues [START_REF] Vaughan | Brillouin scattering, density and elastic properties of the lens and cornea of the eye[END_REF][START_REF] Scarcelli | Confocal Brillouin microscopy for three-dimensional mechanical imaging[END_REF], bone [START_REF] Lees | Studies of compact hard tissues and collagen by means of brillouin light scattering[END_REF][START_REF] Mathieu | Micro-Brillouin Scattering Measurements in Mature and Newly Formed Bone Tissue Surrounding an Implant[END_REF] or tumors. [START_REF] Margueritat | High-frequency mechanical properties of tumors measured by Brillouin light scattering[END_REF] We prepared thick slices of teeth extracted in accordance with the ethical requirements for ex-vivo studies. We obtained maps of Brillouin frequency shift and linewidth that can be interpreted as maps of sound velocity and viscosity. Our results show significant changes between healthy tissues and pathological lesions. Such results allow delineating precisely destructed dentin, paving the way to minimally invasive strategies. We also show maps of the sagittal and transversal planes of healthy tubular dentin that reveal its anisotropic microstructure with a 1 µm resolution, several orders of magnitude below previous reports. [START_REF] Kinney | The mechanical properties of human dentin: a critical review and re-evaluation of the dental litterature[END_REF] Our observations indicate that the collagen-based matrix of dentine is the main load bearing structure, which can thought of as fiber-reinforced composite. In the vicinity of polymeric tooth-filling materials, we observed fingering of the adhesive complex into the opened tubules of healthy dentine. We confirmed our observations with two-photon excitation fluorescence (2PEF) and second harmonic generation (SHG) microscopies that reveal teeth microstructure due to two-photon excited autofluorescence and collagen-specific second order nonlinear optical susceptibility [START_REF] Oheim | Principles of two-photon excitation fluorescence microscopy and other nonlinear imaging approaches[END_REF]. The ability to probe the quality of this interfacial layer could lead to innovative designs of biomaterials used for dental restorations in contemporary adhesive dentistry, and would have direct repercussions on the decision-making during clinical work.

MATERIALS AND METHODS.

Sample preparation and staining.

Teeth were collected at Faculty of Medicine, School of Dentistry, at the Department of Oral Surgery of Dental Clinic of Vojvodina, Serbia. They were extracted for medical reasons and, after the patient's consent, approved to be used for research purposes according to the protocol of Ethical board of Dental Clinic of Vojvodina, Serbia. All types of permanent teeth were analyzed, including incisors, canines, premolars and molars. The teeth were classified in the following groups: healthy teeth with intact hard tissues (two subgroupsteeth cut perpendicularly to their longitudinal axis, teeth cut longitudinally), teeth with caries lesions in dentin [score 5 on the International Caries Detection and Assessment System (ICDAS) scale; distinct cavity with visible dentin], restored teeth with composite resin.

Teeth restored with composite used for the imaging of the dentin-adhesive interface were prepared according to the standard dental adhesive placement procedure, after the etching procedure by using 37% orto-phosporic acid for 15 seconds, rinsing it, and placing a single bond universal adhesive (3M Deutschland GmbH, LOT 663414), before restoration with the Filtek Ultimate Body composite. For 2PEF imaging of dentin-adhesive interface, a water solution containing 0.5% eosin-Y was added in small quantities into the adhesive in order to distinguish the adhesive from the surrounding dentinal tissue.

Teeth were sectioned into ~0.5 mm thick slices along the horizontal and longitudinal axes, by a watercooled, low speed diamond saw using a hard tissues microtome. Samples were wet-polished with abrasive sandpaper sheets of silicone-carbide of increasing grit numbers (400-2000). This was followed by polishing with a diamond paste on a soft fabric, in order to eliminate the surface corrugations. The smear layer produced by polishing was removed by 10% polyacrylic acid used for 5 seconds, after which the samples were rinsed, slightly dry-aired and disinfected in 0,5% chloramine solution and preserved in plastic bags until examined.

Nonlinear laser scanning microscopy

Teeth specimens were observed by a nonlinear optical microscopy setup at the Institute of Physics in Belgrade. The homemade nonlinear laser-scanning microscope used in this study has been described in detail elsewhere. [START_REF] Rabasović | Nonlinear microscopy of chitin and chitinous structures: a case study of two cave-dwelling insects[END_REF][START_REF] Bukara | Mapping of hemoglobin in erythrocytes and erythrocyte ghosts using two photon excitation fluorescence microscopy[END_REF] Briefly, a Ti:Sapphire femtosecond laser (Coherent, Mira 900-F) tunable in 700-1000nm range, and a Yb:KGW (Time-Bandwidth Products AG, Yb GLX) femtosecond laser at 1040nm were used as a laser light sources. The samples were imaged using Carl Zeiss, EC Plan-NEOFLUAR, 40×1.3 oil immersion objective or Plan-APOCHROMAT 20x/0.8 air objective for laser focusing and collection of the signal. A visible interference filter (415 nm-685 nm) positioned in front of the detector is used to remove scattered laser light in fluorescence images. Thus, the whole visible range has been detected in two-photon fluorescence images. The tooth auto-fluorescence was excited by the Ti:Sapphire laser at 730nm. The Eosin-Y fluorescence was excited by the Ti:Sapphire laser at 820nm. Narrowband interference filters (10 nm FWHM) were used for second harmonic generation imaging. The SHG filter central wavelength was 520nm (FB520-10, Thorlabs) for the Yb:KGW laser.

The SHG filter central wavelength was 420nm (FB420-10, Thorlabs) for Ti:Sapphire laser when tuned at 840nm. The SHG excitation wavelength and filter were selected so that fluorescence leakage to SHG channel is minimal.

Brillouin spectroscopy

BLS uses the scattering of light by phonons. Thermal phonons (incoherent vibrations) produce a broad peak at low frequencies due to Rayleigh scattering, while acoustic phonons (coherent vibration) with wavenumbers that fulfill the Bragg condition produce well-defined peaks (Figure 1). Using a monochromatic light beam of wavelength 𝜆 in a normal incidence backscattering geometry, forward propagating phonons produce a positive shift at a frequency 𝑓 = 2𝑛𝑉/𝜆 (anti-Stokes peak), where 𝑉 is the sound velocity and 𝑛 is the refractive index, typically in the GHz range (backward propagating phonons creates a negative Stokes peak at -𝑓).

Assuming that the sample is homogeneous over the scattering volume, the BLS spectra can be interpreted as the response of a damped harmonic oscillator with frequency position 𝑓 and full width at half maximum Γ. In this approximation, the storage modulus is ′ = 𝜌(𝑓𝜆/2𝑛) 2 , where 𝜌 is the mass density. The linewidth Γ is ascribed to the loss modulus 𝑀 ′′ = 𝜌Γ𝑓(𝜆/2𝑛) 2 . In other words, the frequency shift and linewidth can be interpreted as the sound velocity (or stiffness) and viscosity of the sample. Note that at the high magnifications we use, Γ can artificially increase due to scattering.

Brillouin microscope

The Brillouin spectrometer is based on a tandem Fabry-Pérot interferometer (JRS Scientific Instruments). [START_REF] Margueritat | High-frequency mechanical properties of tumors measured by Brillouin light scattering[END_REF] In order to reduce phototoxicity, we used a 𝜆 = 647.1 nm continuous-wave laser as a light source (instead of classical designs based on a frequency-doubled YAG laser at 532~nm) which gives a 𝜆 2𝑛 ≈ 200 nm acoustic resolution, where 𝑛 = 1.6 is the typical refractive index of dentin. The spectrometer offers a 40 GHz free spectral range sampled with a 30~MHz sampling rate. Appropriate signal analysis allows reaching a 6 MHz accuracy at a -10 dB noise level. It is equipped with two sets of mirrors with a 95% reflectivity, and an avalanche photodiode (Count ® Blue, Laser Components) with a dark count lower than 10 counts/sec (typically 2counts/sec) and a quantum efficiency of 65% at this wavelength. The spectrometer is coupled with an inverted life science microscope (Eclipse Ti-U, Nikon) equipped with a micro-positioning stage (Mad City Labs). The collimated linearly-polarized laser light is focused onto the sample using a 100 × objective lens (NA 0.9), offering a ~1 µm resolution. The typical power used at the focus was 5~mW, resulting in a power density over the scattering volume a hundred times smaller than that previously reported in tissues and cells, again to ensure low phototoxicity. The backscattered light is collected by the same objective lens and is collected outside the microscope using a telescope. A particular attention has been paid to respect the f-number (f/18) of the spectrometer in order to maximize photon collection on the photodiode. Spectra were acquired by averaging typically over ~100 acquisitions.

Data analysis

Raw spectra sampled at 30 MHz are recorded by the proprietary software of the spectrometer (Ghost software, Table Stable Ltd., Switzerland). When only one peak is observed on each side of the Rayleigh line (see Fig. 1b for instance), we fit separately the anti-Stokes and Stokes peaks to a Lorentzian function using Matlab and extract amplitude 𝐴, frequency shift 𝑓 and linewidth Γ. Values for each peak are then averaged. When two closely-spaced peaks are observed, we window sequentially around each peak and fit them individually with a single Lorentzian function. This procedure might underestimate slightly the true spacing between the two peaks, but it reducing the error of the fit and hence allows automating the fitting procedure to analyze large images.

The data collected on the different groups of samples were statistically analyzed using Matlab software. The level of significance (𝑝-value) was evaluated according to an unpaired two-tailed 𝑡-test.

RESULTS.

Mechanical signature of caries lesions

Brillouin spectra were recorded for healthy tissues (𝑛 = 34 in 5 different teeth) and caries lesions (𝑛 = 39 in 5 different teeth). Care was taken to focus the laser beam in the regions that didn't show tubules.

For illustration, we plot in Fig. 1 a typical spectrum. Figure 1a shows that the raw spectrum is dominated by Rayleigh scattering (grey box). Figure 1b shows a zoom-in on the Stokes (negative frequencies) and anti-Stokes (positive frequencies) peaks for the healthy tissue and caries lesion. We fit these peaks with Lorentzian functions (dashed lines) to extract the frequency shift 𝑓 and linewidth Γ (Methods). The distributions of frequency and linewidth values are plotted in Figs. 2a and2b, respectively. Note that the spreading of the data is larger in the case of caries because BLS features depend on the position within the highly heterogeneous lesion, as discussed below. Frequency shift is clearly different between healthy and diseased conditions. We don't however observe any significant variation in the linewidth. This observation suggest that BLS features could be used to delineate lesions.

Figure 1: typical spectrum. (a) Raw spectrum where the contribution of Rayleigh scattering is indicated by the grey box. (b) Zoom-in on the Stokes (negative frequencies) and anti-Stokes (positive frequencies) peaks for a healthy tissue (black, intensity multiplied by 2 for clarity) and a caries lesion (blue line). Fits by Lorentzian functions are superimposed (red dashed lines).

To investigate this, we scanned the laser along a line across the frontier between healthy tissue and a caries lesion. Figure 3 shows the frequency shift 𝑓 and linewidth Γ measured from the healthy region (left) to the core of the lesion (right). We see that, as already observed in Fig. 2b

Mechanical anisotropy of the dentin

Healthy dentin is composed of a mineralized matrix perfused by microtubes called tubules. This structure leads to an anisotropic distribution of mechanical properties, largely documented in the literature on a macroscale. [START_REF] Kinney | The mechanical properties of human dentin: a critical review and re-evaluation of the dental litterature[END_REF] To investigate the anisotropy of the healthy tubular dentin on a microscopic scale, healthy teeth were cut longitudinally or transversally (Methods). The microstructure of the healthy dentin in the sagittal and transversal planes is well revealed by the combination of 2PEF image obtained with the Ti:Sapphire tuned at 730 nm and SHG image obtained with the Ti:Sapphire tuned at 840 nm (Figures 4 and5). The typical tubular structure of dentin appears clearly on the 2PEF image. Collagen type I, as a triple-helical molecule, the organic component on which the dentin is mostly built, is aligned in a non-centrosymmetric manner, and is responsible for the contrast in the SHG images. The SHG images of dentin thus reveal the intertubular dentin, which is less mineralized than the peritubular dentin. The exposed surfaces in the sagittal and transversal planes were mapped with BLS. The laser was scanned in a raster pattern with a 1 µm resolution and spectra recorded at each pixel. Resulting images obtained in the sagittal (along the tubules) and transversal planes (perpendicular to the tubules) of the tooth are shown in Figs. 6 and7, respectively. In the sagittal plane, the tubules are clearly visible on the optical image (Fig. 6a), and they are characterized by a lower amplitude (Fig. 6b) and a lower 𝑓 (Fig. 6c). We define an amplitude threshold 𝐴 𝑡ℎ = 0.5𝐴 𝑚 , where 𝐴 𝑚 is the maximum amplitude, that separates the tubules from the surrounding matrix. Pixels in the 𝑓 and Γ images that correspond to an amplitude 𝐴 < 𝐴 𝑡ℎ are ascribed to the tubules, and the remaining pixels to the matrix. The distribution of pixels for 𝑓 and Γ are shown in Figs. 6e and6f. While the linewidth doesn't allow distinguishing the two components (as observed on Fig. 6d), 𝑓-values are slightly higher for the tubules.

In the transverse plane, the tubules appear as disks on the optical image (Fig. 7a). They are also clearly visible on the BLS images, characterized here too by a low amplitude (Fig. 7b). In this case however,

Characterization of the dentin-resin interface

In this section we examine the adhesive interface between healthy dentin and tooth-filling resin. We probed the interface between healthy dentin and resin with BLS. Figure 9a shows a typical spectrum where we recognize the signature of healthy dentin at around 20 GHz (blue area). In addition, we also observe a second peak at around 13 GHz (pink area), revealing the presence of the adhesive.

This result demonstrates that BLS is able to detect the presence of the two components mixed in the hybrid layer. To investigate this, we recorded BLS images in a typical hybrid layer. The frequency shift (Fig. 9b) and linewidth (Fig. 9c) both allow to observe the adhesive, hybrid layer and dentin. This clearly highlights the ability of BLS to quantify adhesive interfaces, and to reveal the fingering of the adhesive in the hybrid layer. Such results should pave the way for novel strategies for tooth repair. 𝜌 ∼ 2100 kg/m 3 is the mass density (Methods). [START_REF] Gradl | Mass Density Measurement of Mineralized Tissue with Grating-Based X-Ray Phase Tomography[END_REF] This matches very well with previous results, ~37-39 GPa, obtained by ultrasonic interferometry. [START_REF] Lees | Anisotropy in hard dental tissues[END_REF] Our data shows an anisotropy of 6% and stiffer values in the transverse plane. Since the tubules are aligned in the softer sagittal plane, this observation demonstrates that the tubules do not contribute to mechanical strength of dentine. Instead, the mineralized collagen fibrils, aligned in the stiffer transverse plane, act as a load bearing structure, as suggested previously. [START_REF] Kinney | The mechanical properties of human dentin: a critical review and re-evaluation of the dental litterature[END_REF] This means that the anisotropy of dentine originates from the anisotropy of the matrix itself, and that its structure could be thought as a fiber-reinforced composite material in the transverse plane.

The difference between frequencies of healthy dentin and those of caries lesions is very clear (𝑝 < 0.0001). The transition from healthy dentin to caries allows delineation with a 20 µm precision for the samples we probed. The linewidth Γ = 2𝜋𝜂(2𝑛/𝜆) 2 is constant across this region, meaning that both 𝑛 and kinematic viscosity 𝜂 are constant. The change in 𝑓 can thus be ascribed to a reduction in the stiffness of the tooth in the lesion. This observation demonstrates the high specificity of BLS, and the ability to produce quantitative indicators of tooth decay. In contemporary clinical practice, the visuotactile method is the dominant approach to caries delineation, and is obviously less resolved than BLS and operator-dependent. In order to avoid inter-examiner variations, visual examination can be aided by dyes, light-fluorescence based methods, or fiber-optic transillumination. In this context, BLS has a unique potential to provide quantitative information at a microscale, in a label-free, non-invasive, nonionizing manner. Further comparison with histological sections of lesions with different scores should reveal the sensitivity of BLS, and grant new minimally invasive strategies in dentistry. As has been done for cancer therapy [START_REF] Margueritat | High-frequency mechanical properties of tumors measured by Brillouin light scattering[END_REF], ophthalmology [START_REF] Yun | Brillouin microscopy: assessing ocular tissue biomechanics[END_REF] or bone, [START_REF] Tsubota | Local ultrasonic wave velocities in trabeculae measured by micro-Brillouin scattering[END_REF] BLS could be used as a platform for early diagnosis of tooth decay, micrometric precise caries detection during clinical work, but also for the monitoring of efficacy of the therapeutic techniques. [START_REF] Sauro | Therapeutic bonding approaches to remineralize dentin-bonded interfaces[END_REF][START_REF] Sauro | Influence of phosphoproteins ' biomimetic analogs on remineralization of mineral-depleted resindentin interfaces created with ion-releasing resin-based systems[END_REF] BLS images of the hybrid layer bring new information on the spatial distribution of mechanical properties in the adhesive junction. In future, BLS could be used to test the adhesion strength of different types of adhesives under preclinical investigation, or the weakening of the hybrid layer by enzyme-caused degradation of collagen or adhesive, mineral-depleted dentin [START_REF] Sauro | Therapeutic bonding approaches to remineralize dentin-bonded interfaces[END_REF][START_REF] Hashimoto | In vitro degradation of resindentin bonds analyzed by microtensile bond test, scanning and transmission electron microscopy[END_REF] or caries lesion [START_REF] Oliveir | The use of micro-Raman spectroscopy for the analysis of caries-aff ected dentin adhesive interfaces[END_REF]. Moreover, the hybrid layer is subjected to multiaxial stresses during functional use, and it has been predicted that stress concentration in the underlying hard dental tissues depends on the stiffness of this interface [START_REF] Misra | Material-tissue interfacial phenomena: Challenges in mathematical modeling[END_REF]. In future, BLS could also help studying the impact of wear on this interface at a microscale. Beyond dentistry, BLS could also provide new tools to investigate the link between microwear, stiffness and diet, offering important information in paleobiology.

  , Γ does not vary. On the other hand, 𝑓 decreases slowly, making the precise delineation difficult. As a reference, we indicate the median (red dashed line) and the 25 th and 75 th percentiles (grey box) measured in Fig. 2a. The upward arrow indicates where 𝑓 falls below the 25 th percentile. This criterion might be used to define the limit of the healthy tissue. In this example, 𝑓 decrease almost linearly by ∼ 0.01 GHz/µm within the lesion. Given the spreading of the data for healthy dentin measured in Fig. 2a, ∼ 0.2 GHz between the mean and the 25 th percentile, such decrease would lead to a 20 µm precision in the identification of the lesion margins.

Figure 2 :

 2 Figure 2: Comparison between healthy tissues (n=34) and caries lesions (n=39). Distributions of (a) frequency shift and (b) linewidth (****p<0.0001, *p<0.05, unpaired two-tailed t-test). The bottom and top edges of each grey box indicate the 25th and 75th percentiles, respectively. The red central line indicates the median. The whiskers extend to the most extreme data points not considering outliers ('+' symbol).

Figure 3 :

 3 Figure 3: scan across a caries lesion. (a) frequency shift and (b) linewidth. The grey box and red dashed line indicate the 25th and 75th percentiles and median identified in Fig. 2a. Upward arrow indicates the position where the frequency shift falls below the 25 th percentile.

[ 28 ]Figure 4 :

 284 Figure 4: 2PEF and SHG images of healthy dentin in the sagittal plane. (a) 2PEF image reveal the tubular structure of dentin, (b) SHG reveals the presence of collagen type I, (c) merged. Scale bars: 20 µm.

Figure 5 :

 5 Figure 5: 2PEF and SHG images of healthy dentin in the transversal plane. (a) 2PEF image reveal the tubular structure of dentin, (b) SHG reveals the presence of collagen type I, (c) merged. Scale bars: 20 µm.

Figure 6 :

 6 Figure 6: images of healthy dentin in the sagittal plane. (a) optical image, (b) amplitude of the Brillouin peaks, (c) frequency shift and (d) linewidth. Panels (e) and (f) show the distribution of 𝑓 and Γ-values measured in panels (c) and (d), respectively. The tubules are identified in panel (b) by amplitude values lower than 50% of the maximum amplitude, and the matrix correspond to remaining pixels.

  the low amplitude of the Brillouin peaks does not allow to identify 𝑓 or Γ-values in the tubular regions, and only the matrix appears in the 𝑓 and Γ images, Figs.7d and 7f. This is because the tubules cause a depression in the transverse plane that defocuses the laser beam. Considering pixels where 𝐴 < 𝐴 𝑡ℎ thus only reveals the distributions of values for the matrix, plotted in Figs.7c and 7e, respectively. The 𝑓 and Γ values differ slightly in the sagittal and transverse planes, suggesting that the matrix itself is anisotropic.

Figure 7 :

 7 Figure 7: images of healthy dentin in the transversal plane. (a) optical image, (b) amplitude of the Brillouin peaks, (d) frequency shift and (f) linewidth. Panels (c) and (e) show the distribution of 𝑓 and Γ-values measured in panels (d) and (f), respectively.

Figure 8a shows a

  typical sample with adhesive interface indicated by the arrow. Merged 2PEF image obtained with the Ti:Sapphire tuned at 820 nm and SHG image obtained with the Yb:KGW (Figure8b) reveals the typical structure of the interface between dentin (SHG, blue pseudocolor) and adhesive (labelled by eosin in 2PEF, red pseudocolor). This image reveals that the low-viscous monomer adhesive fingers in-between the exposed dentinal collagen fibers of intertubular dentin, or in the tubules themselves, forming the so-called hybrid layer. Extensions of monomer entering a few individual dentinal tubules, i.e. resin tags, are also observable (a typical one is pointed at by a white arrow).

Figure 8 :

 8 Figure 8: structure of a typical dentin adhesive interface. (a) photograph of dentin-resin adhesive interface shown by dashed line and the arrow. (b) Merged 2PEF and SHG images of a dentin-adhesive interface. 2PEF micrograph reveals the adhesive colored by eosin Y (red pseudocolor) and SHG micrograph reveals dentin (blue pseudocolor). Arrow points to a typical resin tag. Scale bar: 20 µm.

Figure 9 : 4

 94 Figure 9: imaging the hybrid layer. (a) typical spectrum at the dentin-adhesive interface. The peaks arising from dentin and adhesive are shown in blue and pink, respectively. (a) Frequency shift and (b) linewidth maps of the hybrid layer.

ACKNOWLEDGEMENTS.

The work was partly supported by the COST Action CA16124 "Brillouin Light Scattering Microspectroscopy for Biological and Biomedical Research and Applications", the Agence Nationale de la Recherche (Grant No. ANR-17-CE11-0010-01), and the Ministry of Education, Science and Technological Development of Serbia (projects TR35020, III 45016, ON 171038). The authors thank Tatjana Maravić and the team of professor Bojan Petrović (Faculty of Medicine, School of Dentistry, University of Novi Sad) for their support in samples preparation. REFERENCES.


M Goldman

An ε-regularity result for optimal transport maps between continuous densities

The aim of this short note is to extend the recent variational proof of partial regularity for optimal transport maps to the case of continuous densities.

Introduction

The aim of this short note is to extend the partial regularity result for optimal transport maps obtained in [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF] to the case of continuous densities (rather than Hölder continuous). The interest lies in the proof rather than in the result in itself since it is known to hold under the weaker assumption that the densities are bounded from above and below [START_REF] Figalli | Partial regularity of Brenier solutions of the Monge-Ampère equation[END_REF]. Indeed, we show that for the squared Euclidean cost, both the variational approach to regularity theory for the Monge-Ampère equation recently developed in [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF][START_REF] Goldman | Quantitative linearization results for the Monge-Ampère equation[END_REF][START_REF] Miura | Sharp boundary ε-regularity of optimal transport maps[END_REF] and the one of [START_REF] Philippis | Partial regularity for optimal transport maps[END_REF] lead to the same result. We must however emphasize that the major achievement of [START_REF] Philippis | Partial regularity for optimal transport maps[END_REF] is the treatment of arbitrary cost functions (see [START_REF] Otto | Variational approach to regularity of optimal transport maps: general cost functions[END_REF] for the extension of the variational approach to that setting). Our main ε-regularity theorem is the following (compare with [START_REF] Philippis | Partial regularity for optimal transport maps[END_REF]Th. 4.3]): Theorem 1.1. Let ρ 0 and ρ 1 be densities with compact support i and T be the optimal transport map from ρ 0 to ρ 1 for the squared Euclidean cost on R d . For every α ∈ (0, 1), there exists ε(α, d) > 0 such that if for some R > 0,

1 (2R) d+2 B 2R |T -x| 2 ρ 0 + 1 -ρ 0 2 L ∞ (B 2R ) + 1 -ρ 1 2 L ∞ (B 2R ) ≤ ε, then T is of class C 0,α in B R . * Université
de Paris, CNRS, Sorbonne-Université, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France, michael.goldman@u-paris.fr i we assume compactness of the supports for simplicity. The statement is valid as soon as an optimal transport map exists (see [START_REF] Villani | Topics in optimal transportation[END_REF]).

1

With this ε-regularity result at hand and arguing as for [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF]Th. 1.1], it is not hard to prove that T is a C 0,α homeomorphism outside of a set of measure zero if ρ 0 and ρ 1 are continuous.

Theorem 1.2. For E and F two bounded open sets, let ρ 0 : E → R + and ρ 1 : F → R + be two continuous densities with equal masses, both bounded and bounded away from zero and let T be the optimal transport map between ρ 0 and ρ 1 . Then, there exist open sets E ⊆ E and F ⊆ F with |E\E | = |F \F | = 0 and such that for every α ∈ (0, 1), T is a C 0,α homeomorphism between E and F .

The proof of Theorem 1.1 follows very closely the proof of [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF]Th. 1.2]. It is based on a Campanato iteration scheme which uses at his heart an harmonic approximation result (see Proposition 3.1). The main difference with [5, Th. 1.2] lies in the iteration argument (see Theorem 3.5 below). Indeed, for continuous densities the linear part of the affine transformations introduced in the excess improvement by tilting estimate do not necessarily converge to the identity. This causes the possible blow-up of the C 1,α -norms.

Notation

In the paper we will use the following notation. The symbols ∼, , indicate estimates that hold up to a global constant C, which typically only depends on the dimension d and the Hölder exponent α (if applicable). For instance, f g means that there exists such a constant with f ≤ Cg, f ∼ g means f g and g f . An assumption of the form f 1 means that there exists ε > 0, typically only depending on dimension and the Hölder exponent, such that if f ≤ ε, then the conclusion holds. We write |E| for the Lebesgue measure of a set E. When no confusion is possible, we will drop the integration measures in the integrals. For R > 0 and x 0 ∈ R d , B R (x 0 ) denotes the ball of radius R centered in x 0 . When x 0 = 0, we will simply write B R for B R (0). We will also use the notation -

B R f := 1 |B R | B R f.
Let ρ 0 and ρ 1 be two densities with compact support in R d and equal mass. We say that T is an optimal transport map between ρ 0 and ρ 1 if it minimizes

W 2 (ρ 0 , ρ 1 ) := min T ρ 0 =ρ 1 R d |T -x| 2 ρ 0 , (2.1)
where by a slight abuse of notation T ρ 0 denotes the push-forward by T of the measure ρ 0 dx. We refer the reader to [START_REF] Villani | Topics in optimal transportation[END_REF] for the existence, uniqueness and characterization of such maps.

3 Proof of Theorem 1.1

Let T be the minimizer of (2.1). As in [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF], the proof of Theorem 1.1 is based on the decay properties of the excess energy

E(ρ 0 , ρ 1 , T, R) := R -2 - B R |T -x| 2 ρ 0 . (3.1)
As already alluded to, the main ingredient for the proof of Theorem 1.1 is the following harmonic approximation result (which by scaling we state for R = 1).

Proposition 3.1. For every 0 < τ 1, there exists ε(τ, d) > 0 and C(τ, d) > 0 such that if E(ρ 0 , ρ 1 , T, 1) + 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 ) ≤ ε, (3.2)
then there exists a function ϕ harmonic in B 1/2 , such that

B 1/2 |T -(x + ∇ϕ)| 2 ρ 0 ≤ τ E(ρ 0 , ρ 1 , T, 1) + C 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 ) (3.3)
and sup

B 1/2 |∇ϕ| 2 E(ρ 0 , ρ 1 , T, 1) + 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 ) . (3.4)
Proof. To simplify notation let E := E(ρ 0 , ρ 1 , T, 1) and

D := 1-ρ 0 2 L ∞ (B 1 ) + 1-ρ 1 2 L ∞ (B 1
) . The claim is an almost direct application of [START_REF] Goldman | Quantitative linearization results for the Monge-Ampère equation[END_REF]Th. 1.5]. Notice first that for i = 0, 1,

W 2 ρ i B 1 , ρ i (B 1 ) |B 1 | χ B 1 + ρ i (B 1 ) |B 1 | -1 2 1 -ρ i 2 L ∞ (B 1 ) . (3.5)
Therefore, [START_REF] Goldman | Quantitative linearization results for the Monge-Ampère equation[END_REF]Th. 1.5] gives the existence of a radius R ∈ (3/4, 4/5), a constant c ∈ R and a couple (ρ, j) solving in the distributional sense the continuity equation ii

∂ t ρ + ∇ • j = 0 on R d × (0, 1) and ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 (3.6) such that the following holds. If Φ solves ∆Φ = c in B R and ν • ∇Φ = ν • 1 0 jdt on ∂B R ,
where ν denotes the external normal to ∂B R , then

B 1/2 |T -(x + ∇Φ)| 2 ρ 0 ≤ τ E + CD
ii note that (ρ, j) is actually the solution of the Eulerian version of (2.1), see [START_REF] Goldman | Quantitative linearization results for the Monge-Ampère equation[END_REF].

B 1/2 |T -(x + ∇ϕ)| 2 ρ 0 E(ρ 0 , ρ 1 , T, 1) d+2 d+1 + 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 ) . Lemma 3.3. For g ∈ L ∞ (B 1 ), every solution solution ϕ of ∆ϕ = g in B 1 and ν • ∇ϕ = 1 H d-1 (∂B 1 ) B 1 g on ∂B 1 , satisfies sup B 1 |∇ϕ| 2 g 2 L ∞ (B 1 ) .
Proof. This follows from global Schauder estimates [START_REF] Troianiello | Elliptic differential equations and obstacle problems[END_REF]Th. 3.16 (iii)] and the fact that if g ∈ L ∞ (B 1 ), then g is in the Morrey space L 2,d-2(1-α) (B 1 ) for every 0 < α < 1 (see [START_REF] Troianiello | Elliptic differential equations and obstacle problems[END_REF]).

We now prove that as in [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF]Prop. 3.6], this estimate implies an "excess improvement by tilting"-estimate. Even though the proof is similar to the one in [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF], we include it for the reader's convenience. Proposition 3.4. For every β ∈ (0, 1) there exist ε(d, β) > 0, θ = θ(d, β) > 0 and C θ (d, β) > 0 with the property that for every R > 0 such that

E(ρ 0 , ρ 1 , T, R) + 1 -ρ 0 2 L ∞ (B R ) + 1 -ρ 1 2 L ∞ (B R ) ≤ ε, (3.7)
there exist a symmetric matrix M with det M = 1 and a vector b with

|M -Id| 2 + 1 R 2 |b| 2 E(ρ 0 , ρ 1 , T, R) + 1 -ρ 0 2 L ∞ (B R ) + 1 -ρ 1 2 L ∞ (B R ) , (3.8)
such that, letting x := M -1 x, ŷ := M (y -b) and then T (x) := M (T (x) -b), ρ0 (x) := ρ 0 (x) and ρ1 (ŷ) := ρ 1 (y), (3.9)

we have

E(ρ 0 , ρ1 , T , θR) ≤ θ 2β E(ρ 0 , ρ 1 , T, R) + C θ 1 -ρ 0 2 L ∞ (B R ) + 1 -ρ 1 2 L ∞ (B R ) . (3.10)
Proof. By a rescaling x = R -1 x, which amounts to the re-definition T ( x) := R -1 T (R x) (which preserves optimality) and b := R -1 b, we may assume that R = 1.

As above, we introduce the notation

E := E(ρ 0 , ρ 1 , T, 1) and D := 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 )
. Let τ ∈ (0, 1) to be fixed later and then ϕ be the harmonic function given by Proposition 3.1. Define b := ∇ϕ(0), A := ∇ 2 ϕ(0) and set M := e -A/2 , so that det M = 1. Using (3.4) from Proposition 3.1 and the mean value property for harmonic functions, we see that (3.8) is satisfied.

Defining ρi and T as in (3.9) we have by (3.8) and (3.7)

- B θ | T -x| 2 ρ0 = - M B θ |M (T -b) -M -1 x| 2 ρ 0 - B 2θ |T -(M -2 x + b)| 2 ρ 0 - B 2θ |T -(x + ∇ϕ)| 2 ρ 0 + - B 2θ |(M -2 -Id -A)x| 2 ρ 0 + - B 2θ |∇ϕ -b -Ax| 2 ρ 0 - B 2θ |T -(x + ∇ϕ)| 2 ρ 0 + θ 2 |M -2 -Id -A| 2 + sup B 2θ |∇ϕ -b -Ax| 2 .
Recalling M = e -A/2 , A = ∇ 2 ϕ(0), and b = ∇ϕ(0), we obtain

θ -2 - B θ | T -x| 2 ρ0 (3.3) θ -(d+2) (τ E + C τ D) + |∇ 2 ϕ(0)| 4 + θ 2 sup B 2θ |∇ 3 ϕ| 2 (3.4) θ -(d+2) (τ E + C τ D) + (E + D) 2 + θ 2 (E + D) τ θ -(d+2) + θ 2 E + C τ θ -(d+2) D,
where we used the harmonicity of ∇ϕ and the fact that E + D θ 2 (recall that θ has not been fixed yet). We may thus find a constant C(d) > 0 such that

θ -2 - B θ | T -x| 2 ρ0 ≤ C τ θ -(d+2) + θ 2 E + C τ θ -(d+2) D.
We now fix θ(d, β) such that Cθ 2 ≤ 1 2 θ 2β , which is possible because β < 1. We finally choose τ 1 such that also Cτ θ -(d+2) ≤ 1 2 θ 2β , which concludes the proof of (3.10).

We may finally prove our main ε-regularity result. As already pointed out in the introduction, it is in this iteration argument that the proof departs from the one in [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF]. Indeed, under the assumption that the densities are merely continuous, the distance to the identity of the linear transformations M in (3.8) are not decaying and we need to compensate the possible blow-up of the cumulated linear transformations by downgrading the C 1,α estimates to C 0,α estimates. A similar argument is used in [START_REF] Philippis | Partial regularity for optimal transport maps[END_REF]. Notice that the Campanato iteration itself is somewhat simpler here compared to [START_REF] Goldman | A variational proof of partial regularity for optimal transportation maps[END_REF]Prop. 3.7] since we do not need to introduce an extra dilation factor at every step to propagate the smallness assumption on the data. Theorem 3.5. For every α ∈ (0, 1), if

E(ρ 0 , ρ 1 , T, 2R) + 1 -ρ 0 2 L ∞ (B 2R ) + 1 -ρ 1 2 L ∞ (B 2R ) 1, (3.11) then T is of class C 0,α in B R .
Proof. By scale invariance, we may assume that R = 1. Let us fix α ∈ (0, 1). By Campanato's theory, see [3, Th. 5.5], we have to prove that (3.11) implies sup

x 0 ∈B 1 sup r≤ 1 2 min b 1 r 2α - Br(x 0 ) |T -b| 2 1. (3.12)
Let us first notice that (3.11) implies that for every x 0 ∈ B 1

E := - B 1 (x 0 ) |T -x| 2 ρ 0 1 and 1 -ρ 0 2 L ∞ (B 1 (x 0 )) + 1 -ρ 1 2 L ∞ (B 1 (x 0 )) 1. (3.13)
Therefore, in order to prove (3.12), it is enough to show that (3.13) implies that for r ≤ Without loss of generality we may now assume that x 0 = 0. To simplify notation, we let

ε := E + 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 )
.

Fix from now on β ∈ (0, 1) and let θ(d, β) be given by Proposition 3.4. Thanks to (3.13), Proposition 3.4 applies and there exist a (symmetric) matrix M 1 of unit determinant and a vector b

1 such that T 1 (x) := B 1 (T (M 1 x)-b 1 ), ρ 1 0 (x) := ρ 0 (M 1 x) and ρ 1 1 (x) := ρ 1 (M -1 1 x+b 1 ) satisfy E 1 := E(ρ 1 0 , ρ 1 1 , T 1 , θ) ≤ θ 2β E + C θ 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 ) ≤ (θ 2β + C θ )ε. (3.15)
If T is a minimizer of (2.1), then so is T 1 with (ρ 0 , ρ 1 ) replaced by (ρ 1 0 , ρ 1 1 ). Indeed, because det M 1 = 1, T 1 sends ρ 1 0 on ρ 1 1 and if T is the gradient of a convex function ψ then T 1 = ∇ψ 1 where ψ 1 (x) := ψ(M 1 x) -b 1 • M 1 x is also a convex function, which characterizes optimality [9, Th. 2.12]. Moreover, since by (3.8), we have

|M 1 -Id| 2 ε and |b 1 | 2 θ 2 ε, if ε is small enough then M 1 B θ ⊆ B 1 and M -1 1 B θ + b 1 ⊆ B 1 , 1 -ρ 1 0 2 L ∞ (B θ ) + 1 -ρ 1 1 2 L ∞ (B θ ) ≤ 1 -ρ 0 2 L ∞ (B 1 ) + 1 -ρ 1 2 L ∞ (B 1 ) ≤ ε. (3.16)
Therefore, we may iterate Proposition 3.4, K > 1 times to find a sequence of (symmetric) matrices M k with det M k = 1, a sequence of vectors b k and a sequence of maps T k such that setting for 1

≤ k ≤ K, T k (x) := M k (T k-1 (M k x) -b k ), ρ k 0 (x) := ρ k-1 0 (M k x) and ρ k 1 (x) := ρ k-1 1 (M -1 k x + b k ), we have 1 -ρ k 0 2 L ∞ (B θ k ) + 1 -ρ k 1 2 L ∞ (B θ k ) ≤ ε (3.17) E k := E(ρ k 0 , ρ k 1 , T k , θ k ) ≤ θ 2β E k-1 + C θ ε, (3.18) |M k -Id| 2 E k-1 + ε, (3.19) 1 θ 2(k-1) |b k | 2 E k-1 + ε.
(3.20)

A simple induction argument shows that from (3.18) we get

E k ≤ θ 2kβ E + C θ k-1 j=0 θ 2βj ε ε (3.21)
and so (3.19) and (3.20) lead to

max(|M k | 2 , |M -1 k | 2 ) ≤ (1 + C √ ε) and |b k | 2 θ 2k ε. (3.22)
In particular, this implies that M K B θ K ⊆ B θ K-1 and M -1 K B θ K + b K ⊆ B θ K-1 so that we may keep iterating Proposition 3.4.

Letting,

A k := M k M k-1 • • • M 1 and d k := k i=1 M k M k-1 • • • M i b i , we see that T k (x) = A k T (A * k x) -d k . By (3.22), max(|A k | 2 , |A -1 k | 2 ) ≤ (1 + C √ ε) k . (3.23)
We first estimate by definition of T k , the fact that det A k = 1 and 1

-ρ k 0 L ∞ (B θ k ) 1, - A * k (B θ k ) |T + A -1 k d k | 2 - A * k (B θ k ) |T -A -1 k A - * k x + A -1 k d k | 2 + - A * k (B θ k ) |A -1 k A - * k x| 2 - B θ k |A -1 k (T k -x)| 2 + θ 2k |A -1 k | 2 |A -1 k | 2 (E k + 1) θ 2k
A * k (B θ k ) |T + A -1 k d k | 2 (1 + C √ ε) k(d+2) θ 2k .
Since α and θ are fixed, if ε is small enough, then 1 + C √ ε ≤ θ -2(1-α) d+2(1+α) so that

θ 2 (1 + C √ ε) d+2 ≤ θ 1 + C √ ε 2α
From this (3.14) follows, which concludes the proof of (3.12).

2 |∇Φ| 2 ERemark 3 . 2 .

 2232 + D. Using (3.6) and integration by parts, it is readily seen that we must have c = -B R (ρ 0 -ρ 1 ) and thus |c| 2 D. Taking ϕ := Φ -c 2d |x| 2 and using triangle inequality we get (3.3) and (3.4). Instead of appealing to [4, Th. 1.5] whose proof is quite long and intricate, one could alternatively give a direct proof of Proposition 3.1 following almost verbatim the proof of [5, Prop. 3.5]. One would only need to replace the use of [5, Lem. 2.2], which required the densities to be Hölder continuous, by Lemma 3.3 below. With respect to (3.3), this would lead to the slightly more quantitative statement

  0 )|T -b| 2 r 2α .(3.14)

( 3 .

 3 23)&(3.21)(1 + C √ ε) k θ 2k .Now if ε is small enough, (3.23) yields B

Acknowledgment

I warmly thank G. De Philippis for suggesting that our proof from [5] could cover the case of merely continuous densities. I also thank A. Zilio for pointing out the reference [8]. Part of this research was supported by the project ANR-18-CE40-0013 SHAPO financed by the French Agence Nationale de la Recherche (ANR) and by the LYSM LIA AMU CNRS ECM INdAM. I also thank the Centro de Giorgi in Pisa for its hospitality.