James Badro

Julien Siebert

Francis Nimmo

An early geodynamo driven by exsolution of mantle components from Earth's core

come

An early geodynamo driven by exsolution of mantle

components from Earth's core

At the present day, the geodynamo is powered primarily by compositional buoyancy9-11 due to the crystallization of the inner core from the outer core, which started around ~1 billion years ago8,12. This creates a conundrum as to the origin of the early field; the inner core is certainly much younger than 3.45 Ga, so a different process for driving an early field is required.

Whether an early dynamo could have been driven by thermal buoyancy alone depends on the power extracted from the core by the mantle, which is uncertain14. It has been suggested7,15-17 that light elements dissolved in the core during core formation could exsolve early in Earth's history as the core cools; the resulting compositional buoyancy would generate enough energy to fuel an early geodynamo. Magnesium exsolution prior to inner core growth has been proposed7,16 as a mechanism paralleling oxygen and/or silicon exsolution after inner core crystallization. The perquisite however is that magnesium must dissolve in iron during core formation.

In order to assess the plausibility of that mechanism, we experimentally investigated the solubility of magnesium in molten iron in equilibrium with basaltic and pyrolitic silicate melts at extremely high temperature. The experiments were performed in a laser-heated diamond anvil cell, where thin pure iron disks were sandwiched between two pyrolite or tholeiite glass disks of identical composition, thickness, and diameter. The assembly was compressed to 35-74 GPa, and laser-heated between 3300 and 4400 K for 30 to 60 seconds. After quench and decompression, thin sections were removed from the center of the laser heated spot using a crossbeam focused ion beam microscope. The thin sections were imaged by high-resolution field-emission scanning electron microscopy, and all showed a coalesced spherical iron ball surrounded by molten silicate (Extended Data Fig. 1), confirming the sample (metal and silicate) was fully molten during equilibration. The composition of the metal and silicate was analyzed using high-resolution electron probe microanalysis (see Methods).

Magnesium solubility in iron takes place according to [START_REF] Ringwood | Chemical Evolution of Terrestrial Planets[END_REF] and that reaction's equilibrium constant is , where T is temperature in K, and P pressure in GPa. The parameters (see Methods) were determined from a leastsquares fit to our data to obtain [START_REF] Rubie | Mechanisms of metal-silicate equilibration in the terrestrial magma ocean[END_REF] where numbers in parentheses are the standard errors on the parameters. Parameter c was found statistically irrelevant, showing that MgO solubility is independent of pressure. The regression is plotted along with the experimental data in Fig. 1 and shows an exceptionally good fit with an R 2 of 0.96. This confirms that reaction (Eq. 1) accurately describes the process of MgO dissolution in iron, and that pressure has no observable effect. Aluminum solubility also takes place and can similarly be quantified (Extended Data Fig. 2), as discussed in Methods. At extreme temperatures, however, the two-component system vanishes to a single homogeneous miscible (solvus) metal-silicate phase4. In that case, reaction (Eq. 1) ceases to describe the system because neither of the phases (metal and silicate) is present. The MgO content of the homogeneous melt is then solely a function of the original bulk composition of the two-phase system.

In order to estimate the amount of MgO that can be dissolved in the core during formation, we ran a series of multistage core formation models5 where the planet was grown to its present mass by iterative accretion and core-mantle differentiation of material (see Methods). The magnesium concentrations in the growing core and mantle were calculated iteratively along with other lithophile (O, Si, Al) and siderophile (Ni, Co, Cr, and V) elements. Over 8000 simulations were performed sampling parameter space, and only geochemically consistent models (where the final concentrations of Ni, Co, Cr, and V in the silicate matches present-day mantle abundances) were retained5.

For core formation without a giant impact, we found a maximum of 0.8 wt.% MgO in the core, in the most favorable (hottest geotherm and deepest magma ocean) case. For a presentday temperature18 at the core-mantle boundary (CMB) of 4100 K, the MgO equilibrium value (saturation threshold) in the core is 1.1 wt.% (Fig. 1b). The core is therefore undersaturated in MgO so that any primordial magnesium dissolved during formation would not exsolve to the mantle.

Then we ran a series of core formation models that involve a final giant impact. In the Moon-forming giant impact scenario6, the impactor is thought of typically as a Mars-sized planetary embryo, but models range from 2.5% to 20% Earth masses19,20. With such a size, the impactor is a differentiated object with a core and mantle, and the temperatures during the impact are sufficiently high that the impactor core and the surrounding silicate mantle turn into a single miscible metal-silicate phase (see Methods). As this dense silicatesaturated metallic object (hereafter called the "hybridized impactor core", HIC) merges with the Earth's core, it strongly increases the latter's lithophile element content. We calculated the composition of the HIC as a function of impactor size (Fig. 2a) by assessing its dilution ratio21 (see Methods) in the magma ocean, i.e. the relative mass of magma ocean that the impactor core interacts with. The amounts of Mg, Si, and O brought by the HIC to Earth's core are plotted in Fig. 2b.

The total MgO dissolved in the core (Extended Data Fig. 3) ranges between 1.6 and 3.6 wt. %. Those values are higher than the 1.1 wt.% saturation value at the present-day CMB, implying that the core became over-saturated in MgO as it cooled. The excess MgO must have exsolved to the mantle and provided a large source of potential energy7,16 to drive an early dynamo. Since MgO solubility depends on temperature but not on pressure, MgO exsolution in the core takes place at the CMB, where the temperature is lowest. As MgO exsolves from the metal, the residue becomes denser and sinks, and is replaced by lighter MgO-bearing metal. This process ensures the entire core is processed at the CMB, so that the equilibrium concentration at the CMB (Fig. 1b) sets the concentration in the whole core.

We estimated the energy released by MgO exsolution by calculating the difference in gravitational energies (ΔE grav ) of the core before and after exsolution, where the gravitational energy in each state is given by [START_REF] Wood | Accretion of the Earth and segregation of its core[END_REF] where G is the gravitational constant, M(r) is the mass of the core comprised below radius r, ρ(r) is the density of the core at radius r, and R the radius of the core.

The energy release depends on how the HIC mixes with Earth's core, as shown by the dependence on ρ(r) in Eq. (3). We investigated two extreme models of mixing: (i) full mixing of the HIC with Earth's core producing a homogeneous core, and (ii) full layering where the HIC sits atop Earth's core (see Methods). The energy release as a function of impactor size is plotted in Fig. 3. In the mixed case, it yields between 1 and 5.5x1029 J. For comparison, the total energy release from inner core growth (latent heat and buoyancy) ranges8,12 between 0.9 and 1.7x1029 J. The layered model provides less energy for small impacts (Fig. 3), but again reaches and exceeds the energy released by inner core growth for Mars-sized impactors or larger.

Since MgO solubility depends only on temperature, the power release and onset time of MgO exsolution depend on the temperature evolution at the CMB, which itself depends on the initial MgO concentration in the core (see Methods). Although the early evolution of CMB temperature is uncertain, as an example we adopt an a priori CMB temperature18 model. Prior to inner core growth, the exsolution rate is high as shown in Extended Data Fig. 4, and generates power in excess of ~3 TW (a conservative estimate of how much power is required to run a geodynamo by compositional buoyancy22) over the course of exsolution (Extended Data Fig. 5). With the onset of inner core growth, the cooling and exsolution rates decrease, and the power drops to ~1 TW (see Methods). In terms of timing, the onset of exsolution occurs once the (decreasing) MgO saturation value at the CMB reaches the concentration in the core (Extended Data Fig. 5). For our nominal model, this occurs ~1 Gyr after Earth's formation with a Mars-sized impact and shifts up to ~2.3 Gyr in the case of a small "fast spinning" impact (see Methods).

Rapid initial cooling following a giant impact may have driven an early thermal dynamo. However, our experimental results show that MgO exsolution likely dominated the core's energy budget in the intermediate period between early, rapid cooling and the onset of inner core growth. This provides a tangible basis for an exsolution-driven dynamo7, as well as a plausible mechanism explaining the uninterrupted geological record of magnetism13 in Earth's rocks and minerals dating to 3.5 Ga or earlier. This mechanism should be relatively ineffective in smaller planets such as Mars or on Earth-sized planets that haven't experienced a giant impact, but for super-Earths, where pressures and temperatures could remain super-solvus for extended periods, it represents a novel method of driving potentially detectable present-day dynamos.

Methods

Magnesium and Aluminum Solubility

The thermodynamic process of mantle component solubility involves the solubility of mantle components in the metal phase (Eq. 4), rather than redox exchange as in the case of siderophile element partitioning. The magnesium concentration in the metal ranges between 0.2 mol% and 1 mol% in our experiments. The equilibrium constant of the dissolution reaction ( 4) is [START_REF] Badro | Core formation and core composition from coupled geochemical and geophysical constraints[END_REF] and its logarithm is proportional to the Gibbs free energy change of reaction (Eq. 4) and can be written [START_REF] Hartmann | Satellite-Sized Planetesimals and Lunar Origin[END_REF] where parameters a, b, and c correspond to the entropy, enthalpy, and volume changes of reaction (Eq. 4), respectively. Those parameters were fit to the data by linear regression, and c was found statistically irrelevant (no pressure dependence), to yield

Similarly, the aluminum concentration of in the metal ranges from 0 (below detection limit, explaining 2 fewer points for the Al plot in Extended Data Fig. 2) to 1.1 mol%. The equilibrium constant of the dissolution reaction ( 8) is ( 9)

and its logarithm is proportional to the Gibbs free energy change of reaction (Eq. 8) and can be written again in the same form as equation (Eq. 6); fitting to the data by linear regression shows once more that c as statistically irrelevant, and we find [START_REF] Labrosse | On cooling of the Earth's core[END_REF] Saturation Conditions at the Core-Mantle Boundary Equations ( 5), ( 7), [START_REF] Lister | The Strength and Efficiency of Thermal and Compositional Convection in the Geodynamo[END_REF], and (10) allow us to calculate the Mg and Al concentration in molten iron as a function of temperature and silicate composition. An important case is that of the equilibrium value in the core at the core-mantle boundary (CMB). As shown above, MgO dissolution in iron has no pressure dependence. This means that MgO exsolves in the coldest part of the core, which is the CMB. The equilibrium value at the CMB is therefore is the MgO saturation value; if the MgO concentration in the core is above saturation, MgO will be exsolved until it reaches that value. Fig. 1b shows the equilibrium value of MgO concentration in the core as a function of CMB temperature, for a core buffered by (i.e. in local equilibrium with) a pyrolitic magma ocean (50 mol.% MgO in the mantle).

Experimental and Analytical

The silicate glasses were produced in an aerodynamic levitation laser furnace. The starting mixes were made by grinding and mixing from pure oxide (SiO 2 , MgO, Fe 2 O 3 , Al 2 O 3 ) and carbonate (CaCO 3 ) components, pressing them into pellets, before fusing them at constant fO 2 at 1900-2100 °C for 5 minutes in a laser furnace using a 120 W CO 2 laser. The fused samples were quenched to glasses, and analyzed for recrystallization, homogeneity, and composition on a Zeiss Auriga field-emission scanning electron microscope (IPGP, Paris). The glass beads were thinned down to 20 μm thick double-parallel thin section, and were processed using a femtosecond laser machining platform to cut disks of identical size for loading in the diamond anvil cell. Spherical iron balls 1-3 μm in size were flattened between two such silicate disks, and constituted the layered starting sample. Pressure was measured from the frequency shift of the first order Raman mode in diamond, measured on the anvil tips. Temperature was measured every second, simultaneously from both sides by spectroradiometry. Electronic laser shutdown operates in ~2-4 μs, and temperature quench occurs in ~10 μs (owing to thermal diffusion in the sample) ensuring an ultrafast quench of the sample.

After decompression, a thin section (20-by-10 μm wide, 1-3 μm thick) was extracted from the center of the laser-heated spot using a Zeiss Auriga crossbeam focused ion beam microscope (IPGP, Paris). The sample was imaged and then transferred to a TEM copper grid, and the metal and silicate phases were analyzed using a Cameca SX-Five electron microprobe (CAMPARIS, Paris) with 5 large-area WDS analyzers. Metal and silicate phases of the run products are large enough (> 5 μm) to perform reliable analyses with an electron probe micro-analyzer (EPMA) on FIB sections.

Metal and silicate phases were analyzed using a Cameca SX100 and Cameca SX FIVE (CamParis, UPMC-IPGP) electron probe micro-analyzers. X-ray intensities were reduced using the CITZAF correction routine. Operating conditions were 15 kV accelerating voltage, and 10-20 nA beam current and counting times of 10-20 s on peak and background for major elements and 20-40 s for trace elements (including Mg and Al in the metallic phases). Pure Fe metal was used as standard for metal. Fe 2 O 3 , SiO 2 , MgO, and Al 2 O 3 were used as standards to measure solubility of oxygen, silicon, magnesium, and aluminum in metal. Diopside glass (Si), wollastonite (Ca), orthoclase (K), anorthite (Al), albite (Na), rutile (Ti) and pure oxides (Fe 2 O 3 , MgO, SiO 2 , CaO and Al 2 O 3 ) were used as standards for the silicate. We verified that the geometry of metal and silicate phases was identical from both sides of the FIB sections, so that the excitation volume of EPMA analyses only samples a single phase. EPMA analyses with 1-2 μm beam size are large enough to integrate the small quench features of metal and silicate phases (<200 nm) and determine their bulk compositions. When a few small metallic blobs were present in the silicate (500 nm-2 μm diameter) special care was taken to avoid them during analysis of the silicates.

Core Formation Modeling

The core of the Earth formed in the first ~50 million years24,25 of the Solar System, by an iterative addition of material to the proto-Earth. The accreting material, consisting of mixtures of iron-rich metals and silicates similar to those found in extra-terrestrial bodies (chondrite parent bodies, HEDs, angrites, etc.), impacted the growing planet. The heat generated by those impacts maintained the outermost portion of the planet in a molten state known as a magma ocean26. At temperatures below the solvus of iron and silicate, the two phases unmix and the metal (twice denser) segregates towards the center and forms the core.

Along with the segregating metal, the siderophile elements are stripped to the core, among which are light elements such as Si and O. The depletion of siderophile elements from the mantle has been widely used to constrain the P-T-composition path of core formation, and has shown that the core formed in a deep magma ocean27,28. As the planet accretes, the magma ocean grows deeper; recent models5 show that the concentrations of Ni, Co, Cr, and V in the mantle satisfy terrestrial observables for final magma ocean depth between 1000 and 1700 km, corresponding to final pressures between 40 and 75 GPa and final temperatures between 3000 K and 4180 K, respectively.

We ran a series of traditional multistage core formation models5 where the planet was accreted to its present mass by 0.1% Earth mass increments, without giant impacts. At each stage the planet grows, and the pressure and temperature of equilibration increase accordingly. The concentrations of Ni, Co, V, Cr, O, Si, and Mg in the core were calculated iteratively during the 1000 steps of the accretion process. The simulations were run for a variety of redox paths (ranging from very reduced to very oxidized), several geotherms (between the solidus and the liquidus of peridotite), and for all possible magma ocean depths ranging from 0% (magma lake) to 100% (fully molten Earth) of the mantle. We forwardpropagated all uncertainties on the thermodynamic parameters governing the partitioning equations, using Monte Carlo simulation. Most models (very deep or very shallow) do not satisfy, within uncertainties, the observed geochemical abundances of Ni, Co, V, and Cr in the mantle and therefore aren't relevant. We only selected the models that do reproduce the geochemical abundances of Ni, Co, V, and Cr in the present-day mantle, and found that the maximum MgO concentration in the core at the end of accretion is 0.8 wt.%.

Giant Impact Modeling

In the Moon-forming giant impact scenario6, the impactor is thought of typically as a Marssized planetary embryo, but models range from 2.5% to 20% Earth masses19,20. With such a size, the impactor is a differentiated object with a core and mantle, and therefore (as opposed to small accretionary building blocks) it does not fully equilibrate with the entire magma ocean, but rather partially equilibrates29 with a small portion2,21 of that magma ocean. The impactor and the magma ocean (in the impact zone) reach tremendous temperatures during the impact, as shown by smoothed-particle hydrodynamic simulations19,20. Even though the temperatures from those simulations can be inaccurate because of intrinsic inaccuracies in the equations of state that they are based on, the minimum temperature19 for the impactor core is 8000 K and that of the magma ocean in the impacted area is 7000 K. Therefore, the system consisting of the impactor core and the surrounding silicate mantle is necessarily always hotter than 7000 K, and turns into a single miscible metal-silicate phase.

We calculated the composition of Earth's core after the giant impact in two steps. First, we modeled the pre-giant impact accretionary phase. The Earth was partially accreted as described in the previous paragraph, until it reached 80 to 99% Earth mass, leaving the planet in the state it was in prior to the giant impact. We only considered the models that reproduce the present-day geochemical abundances of Ni, Co, V, and Cr in the mantle. Then the final accretion event took place, consisting of the giant impact bringing in the remaining 1 to 20% of Earth's mass. We calculated the composition of the "hybridized impactor core" (HIC) as a function of its size (Fig. 2a) considering the fact that as opposed to small accretionary building blocks, the core of the giant impactor does not fully equilibrate with the entire magma ocean; it rather partially equilibrates29 with a small portion2,21 of that magma ocean (see Partial Core Equilibration and Turbulent Fragmentation and Mixing section below). It is clear from Fig. 2a the bigger the impactor, the smaller the relative mass of magma ocean it interacts and equilibrates with, and consequently the less mantle components (Mg, O, Si) the HIC contains. The net effect on Earth's core, once the HIC is added, is mitigated as shown in Fig. 2b; it is the result of the balance between larger HICs being less enriched in mantle component, but contributing more mass to the whole core.

Partial Core Equilibration and Turbulent Fragmentation and Mixing

The composition of the HIC was calculated by taking into account two main parameters that are usually neglected in traditional core formation models2,3,5,28,30.

First, the degree of partial equilibration, i.e. the fraction of the core that equilibrates with the mantle, has been constrained by geochemical modeling, from the combined analysis of the Hf-W and U-Pb isotopic systems, and shown to be at least25,29,31 40%. We used that conservative lower bound, meaning that 60% of the impactor core merges with Earth's core without equilibration (and therefore with no compositional effect) whereas the other half equilibrates in the magma ocean before merging with the core.

Second, the impactor core only "sees" a portion2 of the magma ocean and that fraction involved in the equilibration was estimated from fragmentation and turbulent mixing scaling laws21; those show that the ratio of equilibrated silicate to equilibrated metal (dilution ratio Δ) in the magma ocean is given by [START_REF] Gubbins | Can the Earth's dynamo run on heat alone?[END_REF] where ρ silicate and ρ metal are the densities of silicate and metal, and δ is the ratio of impactor to Earth mass.

MgO Exsolution Energy

The energy release depends on how the HIC mixes with Earth's core, as shown by the dependence on ρ(r) in Eq. ( 3). Even though simulations20 and energetic arguments32 both suggest that the HIC should thoroughly mix with Earth's core, we investigated two extreme models of mixing: (i) full mixing of the HIC with Earth's core producing a homogeneous core, and (ii) full layering where the HIC sits atop Earth's core.

In the mixed case, the HIC is diluted in the bulk of Earth's core and therefore the Si and O content delivered by the impactor are below the saturation limit of those elements5,30,33 (Fig. 2b); those concentrations are under-saturated with respect to the overlying conditions imposed by the magma ocean at the core-mantle boundary, and there is no chemical drive to force those components out of the system. In that case, we considered that MgO is the only phase to exsolve so that the associated energy release is a conservative lower bound.

In the layered case, the HIC is concentrated atop the proto-core, and all three mantle components (MgO, SiO 2 , and FeO) are highly concentrated in the layer, and over-saturated with respect to core-mantle boundary conditions prevailing atop that layer. In that case, all of those components would exsolves and remix with the overlying magma ocean.

In our energy calculations, we fixed the present-day CMB temperature to 4100 K. Lower temperatures imply a lower saturation level in the core, and mean that more MgO exsolves and more energy is produced, and conversely. The final density and radius of the core are the present-day values (10.6 g/cm 3 and 3485 km, respectively).

Impactor Core Mixing

We considered a uniform core of density ρ and radius R; it subsequently undergoes unmixing into an inner (dense) region with density ρ c and radius R c (the present-day values above), and an outer buoyant layer with density ρ layer . The volume fraction of the outer layer is f, which we take to be <<1. We may write [START_REF] Labrosse | Thermal evolution of the core with a high thermal conductivity[END_REF] and [START_REF] Tarduno | Solar Wind, and Magnetopause 3.4 to 3.45 Billion Years Ago[END_REF] where equation ( 13) is correct to first order in f. In practice, we specify ρ c and ρ layer (4.8 g/cm 3 for MgO) and calculate ρ and R for a given value of f, with the current core boundary taken to be R c . The gravitational energy E of the core in either state may be derived using (Eq. 3), and the change in energy ΔE in going from the uniform to the unmixed state can be available to do work (e.g. to drive a dynamo). Making use of equations ( 3), ( 12) and ( 13), it may be shown that, to first order in f: [START_REF] Labrosse | A crystallizing dense magma ocean at the base of the Earth's mantle[END_REF] For f = 20%, equation ( 14) overestimates the full calculation (plotted in the figures) by about 5%; the discrepancy is smaller with smaller f, and equation ( 14) can be used to a good approximation to estimate the amount of energy released by lithophile element exsolution from the core. This equation shows the correct limiting behavior in the cases of f=0 and ρ c =ρ layer .

Impactor Core Layering

In this case we take the mass fraction of the Earth's core added by the HIC to be f m .

Assuming the HIC density to have a density equal to ρ i and the present-day total core mass to be M c , the radius of the base of the impactor layer R 1 prior to unmixing of this layer is

given by ( 15)

The HIC layer then undergoes unmixing into two components: "mantle components" (ρ 2 , 5.6 g/cm 3 ) and "core material" (ρ 1 , 10.6 g/cm 3 ). The HIC density ρ i may then be derived using [START_REF] Stevenson | Core Exsolution: A Likely Consequence of Giant Impacts and a Likely Energy Source for the Geodynamo[END_REF] where R 2 is the radius of the base of the light element layer after unmixing. To make the total core mass correct, the density of the pre-impact core, ρ c , is also calculated. Once ρ i , R 1 and R 2 have been calculated, the energy change due to unmixing within the layer can be calculated using successive applications of equation ( 14) as before.

Thermal Evolution and Exsolution Power

Using a CMB temperature evolution model, we can estimate the MgO exsolution rate, and hence an exsolution power, as a function of time. A typical CMB temperature evolution is shown in Extended Data Fig. 4a, along with the associated MgO content of the core (Extended Data Fig. 4b) obtained by rewriting the MgO equilibrium curve (Fig. 1b) as a function of time. The time derivatives are the core's cooling rate and its MgO exsolution rate as a function of time, and are plotted in Extended Data Fig. 4c and4d, respectively.

Very early in Earth's history, the equilibrium MgO concentration at the CMB (Extended Data Fig. 4b) is higher than the core's MgO content, and no exsolution occurs. The reverse reaction, i.e. the potential for MgO to be dissolved from the mantle to the core is limited; it is prone to only affect a thin layer below the CMB that is enriched in MgO, becomes light and stably stratified, and therefore incapable to recycle and affect the entire core. As the core cools, exsolution starts once the temperature at the CMB reaches a critical value corresponding to an MgO equilibrium concentration equal to that in the core. This is shown in Extended Data Fig. 5, and is highlighted for two models: the Mars-size impact19 leaving behind a core containing 2.9 wt% MgO and a small "fast-spinning" impact20 producing a core containing 2.1 wt% MgO (see Fig. 2b and Extended Data Fig. 3). The power produced by MgO exsolution is linked to exsolution rate, and it can be estimated from the energy release (Fig. 3 and Extended Data Fig. 8) to be between 5.5 and 7 TW/wt.%/Ga. This allows us to translate an exsolution rate (Extended Data Fig. 4d) into exsolution power, as shown in Extended Data Fig. 5b and8.

What is noteworthy is that initial MgO core content doesn't directly affect exsolution power. The latter is only a function of exsolution rate, itself a function of core cooling rate. Initial MgO content only sets the onset of exsolution as shown in Extended Data Fig. 5. Of course, higher MgO contents in the core entail an earlier onset of exsolution, a longer duration for buoyancy-driven exsolution power, and hence much higher total exsolution energies as clearly shown in Fig. 3. Note that this dichotomy could be mitigated had we self-consistently included MgO exsolution in the thermal evolution model of the core. MgO exsolution power dramatically drops with the onset of inner core growth, as a consequence of the drop in core cooling rate. At the present day, MgO exsolution should still produce ~1 TW of power, much lower than the ~3 TW produced by inner core growth and driving the geodynamo. However, prior to inner core growth, exsolution power is always higher than ~3 TW, demonstrating that MgO exsolution can conceivably drive a geodynamo as early as ~1 Ga after core formation, and until the onset of inner core growth.

The Geodynamo

Assuming an entirely bottom-driven present-day dynamo, corresponding to a core-mantle boundary heat flow exactly at the adiabatic value Q ad of 15 TW34,35, the convective power sustaining the geomagnetic field P = ε Q ad is 3 TW, where ε=0.2 is the thermodynamic efficiency of latent heat and light element release at the inner core boundary22. Power-based scaling laws of the magnetic intensity36 then predict an internal magnetic field of about 1-4 mT, the higher estimate being in agreement with the observation of magnetic Alfvén waves in the core37 coupled to length-of-day variations at periods close to 6 years38.

Dynamo strength increases as buoyancy flux increases39,40, so the MgO exsolution mechanism represents a potent driver of an early geodynamo7. Although a giant impact might cause thermal stratification in the core12,41, the stabilizing thermal buoyancy will be completely overwhelmed by the compositional buoyancy associated with MgO exsolution.

Extended Data

Extended Data Figure 1. A fully molten metal-silicate sample recovered from the laser-heated diamond anvil cell.

A backscattered electron scanning electron microscopy image of a thin section recovered from a laser heated diamond anvil cell experiment. The section is excavated and lifted out from the center of the heated region then thinned down to 3 microns using a focused ion beam instrument. The metal and the silicate are both fully molten, as indicated by the coalesced metallic ball in the center and the circular rim of silicate around it. This sample was compressed to 55 GPa and heated to 3,600 K for 60 seconds.

Extended Data Figure 2. Magnesium and aluminum solubility in metallic iron melt at high pressure and temperature.

(top) Equilibrium constant for MgO dissolution in molten iron as a function of reciprocal temperature. The circles correspond to the experimental data (Extended Data Table 1) and error bars to standard error, while the squares correspond to the low-temperature extrapolation of DFT calculations4. The thick line corresponds to the least-squares linear fit to the experimental data alone (Fig. 1); it shows the remarkable agreement between the theoretical and experimental datasets, especially at high temperature where the theoretical dataset (which is extrapolated from higher temperatures) is the least influenced by extrapolation. (bottom) Equilibrium constant for Al 2 O 3 dissolution (see Methods) in molten iron as a function of reciprocal temperature. The circles correspond to the experimental data (Extended Data Table 1) and error bars to standard error. The thick line corresponds to the least-squares linear fit to the data (R 2 =0.92), and we find .

Extended Data Figure 3. Total MgO dissolved in the core after the giant impact.

A companion to Fig. 2, it is the sum of the MgO component dissolved in the core prior to the impact (0.8 wt.%) and that brought by the HIC. The 10% Mars-size impact19 and 2.5%

"fast-spinning" impact20 are highlighted by circles, and provide 2.9 wt.% and 2 wt.% MgO to the core, respectively. (a) The MgO equilibrium concentration in the core (same figure as Extended Data Fig. 4b), corresponding to our nominal CMB temperature evolution. The onset of MgO exsolution from the core occurs when the MgO equilibrium concentration drops below the MgO content in the core, which is reported here in two cases: 2.9 wt.% for the Mars-sized impactor, and 2.1 wt.% for the "fast-spinning" impactor. For the thermal evolution model in Extended Data Fig. 4a, this onset is at 1.1 Ga and 2.3 Ga, respectively. (b) and (c) Exsolution power for these two cases which is proportional to the MgO exsolution rate plotted in Extended Data Fig. 4d. It is noteworthy that power at a given time is independent of initial MgO content (as long as MgO is being exsolved). The latter only affects the onset of exsolution and therefore the duration of energy release. We also note that the power produced is in excess of 3 TW, and therefore sufficient to drive a dynamo by compositional buoyancy. Finally, the power drops dramatically with the onset of inner core growth, because of the associated drop in core cooling rate and MgO exsolution rate. The core at the present day is still exsolving MgO, and should produce ~1 TW of power, significantly less than the power produced by inner core growth.

Extended Data Figure 6. Equilibrium Mg and MgO concentration in the core as a function of CMB temperature.

This is obtained by rewriting into with X MgO =0.5 (pyrolitic mantle). This curve (red for MgO, blue for Mg) allows determining magnesium saturation in the core at a given temperature. This threshold is important to (i) estimate the present-day MgO content of the core, and hence how much MgO was lost by exsolution over geologic time (Extended Data Fig. 4), and to (ii) estimate the temperature at which MgO exsolution started after core formation (Extended Data Fig. 5). For instance, for a core containing 2.9 wt.% MgO (for a Mars-sized impact, see Extended Data Fig. 3), exsolution is not bound to occur until the temperature at the CMB cools below 5030 K. Moreover, if the present-day CMB temperature is 4100 K, the MgO saturation in the present-day core is 1.1 wt.%, so that the total amount of MgO that can be exsolved from the core isn't the total initial MgO content, but that amount minus the present-day saturation value.

Extended Data Figure 7. Chemical effect of equilibration of the impactor's core in Earth's magma ocean.

Another companion to Fig. 2, a plot showing the "swelling" of the impactor core to form the hybridized impactor core (HIC). The HIC is larger than the impactor core because of the dissolved mantle components therein, which can represent up to two times its initial mass. This y-axis shows the "swelling" factor, e.g. the ratio of HIC to impactor core . This is equivalent to an effective dilution ratio. Small impactors interact with larger relative fractions of the magma ocean; therefore they incorporate more mantle components per unit mass than large impactors, and "swell" more. The HIC of a "fast-spinning" impactor20 (2.5% Earth mass) is 2.2 times larger that the original impactor core, with 45% of its mass made of initial impactor core material (iron) and the remaining 55% consisting of magma ocean components, as shown in Fig. 2a. On the other hand, the core of a Mars-sized impactor19 (10% Earth mass) is 60% larger after equilibration with the magma ocean. A companion to Fig. 3, the gravitational energy released by exsolution is converted into average power, assuming a characteristic time of exsolution of 1 Gyr. Again, the red curve corresponds to the energy released if the HIC fully mixes with Earth's core, and the blue curve corresponds to the energy released if the HIC forms a layer on top of Earth's core. The gray horizontal band corresponds to 3 TW, the power driving the dynamo today, and thus provides a conservative estimate as to how much power is required to run a geodynamo by compositional buoyancy22 [ref 37]. The 10% Mars-size impact19 and 2.5% "fast-spinning" impact20 are highlighted by circles. Note that the blue curve represents a lower-bound to the energy released in case of layering of the HIC, because the layer contains so much lithophile elements that it would exsolve much faster, producing more power, albeit during a shorter period. By proportionality, this plot can be used to infer the power release for any characteristic exsolution time.

Extended Data Table 1 Analyses of the Mg and Al concentrations in the metal and silicate phases of the experimental runs.

Experimental conditions (pressure in GPa, temperature in K, uncertainties in parentheses) and phase composition; all compositions are in molar fractions, and standard errors are 1sigma. The values for log K D are plotted in Fig. 1 and Extended Data Fig. 2. Full chemical analyses of the samples can be found in an excel spreadsheet available online. to obtain , where T is temperature in K, and then converting Mg molar fractions to MgO weight fractions. This is the saturation MgO concentration in the core at a given temperature, and shows that for a present-day CMB temperature of 4100 K, the core cannot contain more than 1.1 wt.% MgO: any MgO dissolved in the core (during core formation) in excess of that value must have exsolved. For an extended version of this graph, see Extended Data Fig. 6. The compositional imprint of the giant impact on the core; between 2 and 8% of the core's total mass consist of mantle material transported by the HIC. The Si and O concentrations added to the core are lower than the amounts present in the core prior to the impact23. This shows that the giant impact's significant contribution to core chemistry is the magnesium influx. The 10% Mars-size impact19 and 2.5% "fast-spinning" impact20 are highlighted by circles. Calculated following equation (3), the red curve corresponds to the energy released if the HIC fully mixes with Earth's core. In that case, MgO exsolution occurs up to the current saturation limit (Fig. 1b). The blue curve corresponds to the energy released if the HIC forms a layer on top of Earth's core. In that case, the layer is so rich in lithophile elements (Fig. 2a) that the exsolution of all dissolved mantle components (MgO and SiO 2 ) takes place. The 10% Mars-size impact19 and 2.5% "fast-spinning" impact20 are highlighted by circles. The gray horizontal band corresponds to the energy release by inner core growth (gravitational + latent heat) since its inception, and is the main driver for the geodynamo today. The energies released by MgO exsolution are on the order of, if not significantly higher than, those released by inner core growth and show the effectiveness of lithophile element exsolution to drive an early dynamo. The average power of exsolution can be estimated assuming an exsolution time (Extended Data Fig. 9) or a temperature evolution model of the core (Extended Data Fig. 5).

Extended Data Figure 4 .Extended Data Figure 5 .

 45 Thermal evolution of the core and MgO exsolution rate.(a) Example CMB temperature evolution as a function of time (after core-formation) from Fig 4a of ref. 18, and its derivative (c), which is the cooling rate. (b) The associated MgO equilibrium concentration in the core, obtained by turning the temperature dependence in Fig. 1b into time dependence; and its derivative (d) which is the exsolution rate. MgO will only start exsolving from the core when the MgO equilibrium concentration (panel b) drops below the MgO content in the core inherited from core formation and the giant impact. Note that core cooling rate and therefore MgO exsolution rate drop dramatically with the onset of inner core growth. The core at the present day is still exsolving MgO, albeit at a much slower rate than that prior to inner core growth. Onset of MgO exsolution and associated exsolution power for two typical models.

Extended Data Figure 8 .

 8 Power released by exsolution if it occurs over 1 Gyr.

Figure 1 .

 1 Figure 1. Magnesium solubility in metallic iron melt at high pressure and temperature.(a) Equilibrium constant for MgO dissolution in iron as a function of reciprocal temperature (Eq. 2). The experimental data is from Extended Data Table1. The line corresponds to the least-squares linear fit to the data. A comparison with extrapolation from DFT calculations4 is shown in Extended Data Fig.2. (b) The resulting MgO concentration in iron in equilibrium with pyrolite as a function of temperature. This is obtained by rewriting (Eq. 2)

Figure 2 .

 2 Figure 2. Composition of the core of the giant impact after equilibration in the magma ocean, and its effect of Earth's core composition.(a) The composition of the hybridized impactor core (HIC), plotted as a function of impactor mass. Smaller impactors interact and equilibrate with larger relative amounts of magma ocean material; they "swell" (see Extended Data Fig.7) and become very enriched in Mg, Si, and O. (b) The compositional imprint of the giant impact on the core; between 2 and 8% of the core's total mass consist of mantle material transported by the HIC. The Si and O concentrations added to the core are lower than the amounts present in the core prior to the

Figure 3 .

 3 Figure 3. Gravitational energy released by the exsolution of buoyant lithophile components from the core after the giant impact.

Nature. Author manuscript; available in PMC 2017 February 18.

Nature. Author manuscript; available in PMC 2017 February 18.Europe PMC Funders Author Manuscripts

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n° 207467. Parts of this work were supported by the UnivEarthS Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02), IPGP multidisciplinary program PARI, and by Paris-IdF region SESAME Grant no.

12015908. JS acknowledges support from the French National Research Agency (ANR project VolTerre, grant no. ANR-14-CE33-0017-01). We thank Julien Aubert, Sarah Stewart, and Paul Asimow for fruitful discussions. We thank Rick Ryerson for valuable comments on the manuscript.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.


Mohammad Nikhian Sadiq
email: mohammad-nikhian.sadiq@univ-brest.fr

Marc Le Roy
email: marc.leroy@univ-brest.fr

André Pérennec

Paul Laurent

Noham Martin

Damien Passerieux

Aurelian Crunteanu

Rémi Boyer

Frédéric Dumas-Bouchiat

Marie-Blandine Martin

A Perennec

M-B Martin

L Divay

Q Levesque

G Tanné

Design and Characterisation of VO 2 Based Switches for Ultra-Fast Reconfigurable Devices

Keywords: Vanadium dioxide, ultra-fast switches, reconfigurable architectures, frequency shift keying

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The new military and civilian standards (e.g. massive MIMO and 5G Network purpose beamforming) and the allocation of the new frequency bands had increased the needs of reconfigurable and agile radiofrequency (RF) devices. In these circumstances, the switching speed of the RF devices is an essential and crucial parameter. Actual devices are, for the most part, based on mechanical technologies like the Micro-Electro-Mechanical Switches (MEMS) or are based on semiconductors components like PIN Diodes, varactors or Field Effect Transistors (FET) [START_REF] Rebeiz | RF MEMS: Theory, Design and Technology[END_REF]. The latter are compact and exhibit a switching time in the range of 50 ns down to 1 ns, but are less interesting over wide band applications where non-linear parasitic effects are hard to take into account. On the other hand, MEMS switching time ranges typically between 500 µs down to 1 µs [START_REF] Rebeiz | RF MEMS: Theory, Design and Technology[END_REF]. In that context, the Vanadium Dioxide (VO 2 ), which is a phase switching ultrawide band material [START_REF] Hariri | Very Wide-Band and Compact VO2 Based Switches[END_REF], has been identified as a promising solution with a switching time in the range of 0.1 µs down to 100 ps [START_REF] Rebeiz | RF MEMS: Theory, Design and Technology[END_REF].

This study focuses on the design of elementary agile RF devices using VO 2 based switches. The RF performances are experimentally characterized under an electrical control and the very low switching time using an optical excitation. The design and complete characterization of these devices are essential to have a set of reliable elementary blocks in order to build more advanced functions (SP2T, SP4T, switched-line phase shifters …).

II. VANADIUM DIOXIDE VO 2 : PROPERTIES AND ACTIVATIONS

The Vanadium Dioxide (VO 2 ) is a Metal-Insulator Transition (MIT) material also known as the Mott Transition [3]- [START_REF] Stefanovich | Electrical Switching and Mott Transition in VO2[END_REF]. When increasing the temperature from ambient temperature, the VO 2 is in its semiconductor state (insulator) until the insulator → metal transition temperature (T cm = 68°C), and becomes conductor above. Decreasing temperature, this reversible transition (metal → insulator) occurs at a slightly lower temperature (T ci = 65°C < T cm ) due to a well-known hysteresis properties. T cm , T ci and electrical resistivity can be adjusted according to material properties.

The VO 2 transition can be obtained either with a thermal [START_REF] Ben-Messaoud | High Contrast Optical Switching in Vanadium Dioxide Thin Films[END_REF], an optical [START_REF] Crunteanu | Microwave Switching Functions using Reversible Metal-Insulator Transition (MIT) in VO2 Thin Films[END_REF] or an electrical [START_REF] Mennai | High Cut-off Frequency RF Switches Integrating a Metal-Insulator Transition Material[END_REF] command.

The thermal command has the disadvantage of being too slow due to the classical thermal inertia needed to heat up and especially to cool down the VO 2 . Conversely, the electrical and optical commands are much faster (potentially under 1 ps for the optical command) [START_REF] Cavalleri | Femtosecond Structural Dynamics in VO2 During an UltraFast Solid-Solid Phase Transition[END_REF] and therefore particularly relevant for the foreseen RF applications. Those devices are implemented with different sizes (VO 2 gap length between 5 µs to 30 µs) to enhance the understanding and the characterization of the material properties and the impact on RF performances.

III. DESIGN AND IMPLEMENTATION

In a complementary approach, VO 2 material has also been integrated into a set of simple RF devices to make them agile. Several stub configurations have been designed: an open circuit one is detailed below (Fig. 1(d)). It consists of a VO 2 based switch in series configuration in order to increase the length of the stub and consequently shift the resonance frequency toward the low frequency, i.e. a shift from 15 GHz to 13 GHz.

In that case, a Molybdene biasing line has been added and optimized to control VO 2 state through an electrical command without degrading the RF performances.

The main process of devices implementations are illustrated in Fig. 2 together with the material characteristics. Moreover, four probes measurements showed a R OFF /R ON ratio of 2.5×10 4 after VO 2 deposition on Sapphire using PLD.

IV. RF RESULTS

The first measurements have been conducted at Lab-STICC using a 200 µm GSG pitched probe station and under voltage command applied between both ends of the VO 2 material. Only, the results of most relevant devices (series switches and frequency tunable stubs) for our final applications are shown here for brevity.

A. Series Switches Specifications

The evolutions of the S 11 and S 21 parameters of the 5 µm VO 2 gap based series switch are given in Fig. 3(a). In their insulator states (OFF state), all the VO 2 based switches understudied (gap of 5 µm, 10 µm, 20 µm and 30 µm) show good isolation level S 21 ≤ -20 dB in the 8 -18 GHz frequency band. In its metallic state (ON state), the VO 2 introduces variable losses with the gap length, i.e. VO 2 length. Those losses are consistent with VO 2 switch state-of-the-art, i.e. within -1.2 dB to -2.0 dB respectively for a gap varying from 5 µm to 30 µm, and with S 11 ≤ -14 dB over 8 -20 GHz frequency band.

Like the resistivity versus temperature curve, I-V curve displays a hysteresis (Fig. 3(b)) with two different commutation voltages: insulator, V ci, and metallic, V cm . The voltage pair (V ci , V cm ) is directly proportional to the length of the VO 2 gap: (5.2 V, 9.6 V) and (7.9 V, 20.2 V) respectively for a 5 µm and a 30 µm gap. Hence a higher voltage is needed to obtain the same current while increasing the VO 2 length (i.e. the corresponding resistance). At this stage on this wafer, 5 µm series switches seem to provide the best trade-off between isolation in OFF state, transmission losses in ON state and required switching power.

B. Stubs Specifications

To control the commutation of the VO 2 when integrated into the frequency tunable stub, an additional DC probe has been used and placed over the biasing pad (Fig. 1(d)).

S 11 and S 21 parameters curves versus frequency are given are given for both states in Fig. 3(c) (VO 2 gap of 5 µm). A slight shift in the resonance frequency has been noticed toward the low frequency (14.4 GHz instead of 15 GHz) in the insulator state of the VO 2 , primarily explained by the R-cut Sapphire anisotropic permittivity. The frequency shift occurred by switching the VO 2 in its metallic state is also slightly higher than expected (2.4 GHz instead of 2 GHz) and can be mainly explained by the lower value of the resistive line's resistivity than used in the simulation. In both state, good isolation parameters S 21 < -15 dB are obtained at both work frequencies A retro-simulation has allowed to obtain the effective permittivity of the R-cut Sapphire (ε r = 10.6) and Molybdene's resistivity (5.7×10 -6 Ω.m).

The measured characteristic curves I-V are given in Fig. 3(d) for stubs with 5 µm and 10 µm VO 2 gap. Voltage pairs (3.7 V, 12.3 V) and (4.8 V, 23.3 V) respectively for 5 µm and 10 µm gaps are noticed. Once more, 5 µm VO 2 gap provides the more relevant trade-off.

V. SWITCHING TIME AND OPTICAL COMMAND Experimentations based on an optical command have also been conducted on 10 µm series switches at Thales Research & Technology (TRT) on the same wafer by using a three wavelength laser (UV 355 nm, G 532 nm and NIR 1064 nm) with a 4 ns long pulse placed over a probe station with three optical lenses (×2, ×10 and ×50).

First, increasing the power of the excitation leads to a longer time in the metallic state (Fig. 4(a)) corresponding to the time the material needs to dissipate the additional energy provided by the laser. A deeper insight into the voltage at the VO 2 's terminals shows a step during the rising slope (Fig. 4(b)). The first slope corresponds to a fast reaction of the material but with an insufficient level of energy to maintain a long-lived conductivity. And the second one corresponds to a slower reaction where the energy is sufficient enough to reach conductivity as explained in [START_REF] Cocker | Phase Diagram of the UltraFast Photoinduced Insulator-Metal Transition in Vanadium Dioxide[END_REF].

Boltzmann sigmoid equation ( 1) has been used to estimate the rising time, where a and b are respectively the low and high horizontal asymptotes, τ the time constant and t 0 the central value.

𝑓(𝑡) = 𝑏 + 𝑎 -𝑏 1 + 𝑒 𝑡-𝑡 0 𝜏 (1)
(Model 1 ) considers only the fast reaction and estimates a 4.4 ns rising time whereas (Model 2 ) approximates both reactions and estimates a 33 ns rising time. Finally, a more complete model (Model 3 ), taking into account the response disruption in the reaction provides a 40.1 ns rising time. This approach is of prime interest to identify the different phenomena that rule of the switching time and to quantify the latter in a real RF switch implementation. From these results, complementary studies are being conducted in order to decrease the switching time.

Moreover, the VO 2 based series switches are observed under a microscope: before any optical excitation, the VO 2 appears blue (Fig. 4(c)). After several excitations using the NIR wavelength, degradations are visible (Fig. 4(d)) and electrical parameters change irreversibly (Fig. 4(a)). This phenomenon might be explained by a partial ablation of the material by the laser excitation.

VI. CONCLUSION

Promising results are shown about RF specifications and rising time of the VO 2 to be used in the design of RF devices for ultra-fast reconfigurable applications. Thanks to these experimental validations, VO 2 switches implementation and optical control are currently under improvement to design ultra-fast SP2T and SP4T optically controlled microwave switches dedicated to advanced functions such as 2 and 3-bits phase shifters.

Fig. 1 .

 1 Fig. 1. (a) Photo of the wafer with the devices, (b) series switch, (c) shunt switch and (d) frequency tunable stub based on VO2.Several devices are studied in this article (Fig.1(a)). Firstly, Coplanar Waveguide (CPW) switches are designed in series (Fig.1(b), i.e. VO 2 thin film between input and output RF ports) and shunt (Fig.1(c) with VO 2 between ground and hot lines of CPW) configurations.

Fig. 2 .

 2 Fig. 2. Cross-section (not to scale) of the different layer implementation.

Fig. 3 .

 3 Fig. 3. (a) S-parameters of the 5 µm VO2 gap series switch, (b) I-V characteristic curve of the 5 µm and 30 µm VO2 gap series switch (in series with a 2 kΩ resistor), (c) S-parameters of the 5 µm VO2 gap stub, (d) I-V characteristic curve of the 5 µm and 10 µm VO2 gap stub (with a 1.3 kΩ resistive line).

Fig. 4 .

 4 Fig. 4. (a) Voltage at the VO2 switches outputs while increasing NIR optical power and (b) model estimation of the switching time (10-90% rising time) with NIR wavelength, (c) Picture of the 10 µm gap VO2 series switch before NIR optical excitation, (d) and after several NIR optical excitations.

ACKNOWLEDGMENT

The authors would like to thank the French National Research Agency (ANR) for the MUFRED project funding and C. Prevot and B. Servet from TRT for optical experiments.


Andrea Mendizabal

Tatiana Fountoukidou

Jan Hermann

Raphael Sznitman

Stéphane Cotin

Stephane Cotin

A Combined Simulation & Machine Learning Approach for Image-based Force Classification during Robotized Intravitreal Injections

Introduction

Intravitreal injections are one of the most frequent surgical interventions in opthalmology with more than 4 million injections in 2014 alone [START_REF] Ullrich | Assistive device for efficient intravitreal injections[END_REF]. This procedure is used, for instance, in the treatment of age-related macular degeneration for injecting vascular endothelial growth factor inhibitors. Similarly, intravitreal therapy is also used in the treatment of diabetic maculopathy and retinopathy. With the increasing prevalence of diabetic patients and aging demographics, the demand for such intravitreal therapy is growing significantly.

At the same time, robotic assistance in ophthalmology offers the ability to improve manipulation dexterity, along with shorter and safer surgeries [START_REF] Meenink | Robot-assisted vitreoretinal surgery[END_REF]. Ullrich et al. [START_REF] Ullrich | Assistive device for efficient intravitreal injections[END_REF] has also proposed a robotized intravitreal injection device capable of assisting injections into the vitreous cavity, whereby faster injections in safer conditions are projected. Designing such robotic systems involves significant hurdles and challenges however. Among others, accurate force sensing represents an important topic, whereby the force required to puncture the sclera with a needle is very small and measuring this force is essential for the safety of the patient.

For many devices, force sensing plays a central role in the control loop and safety system of surgical robots [START_REF] Weber | Instrument flight to the inner ear[END_REF]. The use of force sensors and force-based control algorithms allows for higher quality human-machine interaction, more realistic sensory feedback and telepresence. Beyond this, it can facilitate the deployment of important safety features [START_REF] Haidegger | Force sensing and force control for surgical robots[END_REF]. Considerable amount of work relying on force sensors has previously been done, focusing on the development of miniaturized sensors to ease their integration with actual systems. In addition, they must be water resistant, sterilizable and insensitive to changes in temperature [START_REF] Haouchine | Vision-based Force Feedback Estimation for Robot-assisted Surgery using Instrument-constrained Biomechanical 3D Maps[END_REF]. The major limitation of standard force sensors is thus the associated cost, since most surgical tools are disposable [START_REF] Haidegger | Force sensing and force control for surgical robots[END_REF]. Alternatively, qualitative estimation of forces based on images has been proposed in the past [START_REF] Mura | Vision-based haptic feedback for capsule endoscopy navigation: a proof of concept[END_REF][START_REF] Haouchine | Vision-based Force Feedback Estimation for Robot-assisted Surgery using Instrument-constrained Biomechanical 3D Maps[END_REF]. Mura et al. presented a vision-based haptic feedback system to assist the movement of an endoscopic device using 3D maps. Haouchine et al. estimates force feedback in robot assisted surgery using a biomechanical model of the organ in addition to the 3D maps. In contrast, including detailed material characteristics allows precise quantitative assessment of forces but calculating the properties of the materials still remains highly complex [START_REF] Aviles | Sensorless force estimation using a neuro-vision-based approach for robotic-assisted surgery[END_REF]. Recent research have proposed the use of neural networks for force estimation, such as in [START_REF] Aviles | A recurrent neural network approach for 3d vision-based force estimation[END_REF] where interaction forces using recurrent neural networks are estimated from camera acquisitions, kinematic variables and the deformation mappings. In a follow up, [START_REF] Aviles | Sensorless force estimation using a neuro-vision-based approach for robotic-assisted surgery[END_REF] used a neuro-vision based approach for estimating applied forces. In this approach tissue surface deformation is reconstructed by minimizing an energy functional.

In contrast, our method consists in estimating the force, based only on an OCT image of the deformation of the eye sclera. The method relies on a deep learning classifier for estimating quantiles of forces during robotized intravitreal injections. Our contribution hence consists of two critical stages. We first build a biomechanical model of the sclera in order to generate virtually infinite force-OCT image pairs for training any supervised machine learning method. In particular, we will show that this allows us to avoid the need for large data sets of real OCT images. Second, we train end-to-end, a two-layer image classifier algorithm with the synthetic images and the corresponding forces. This allows a very simple estimation process to take place, without the need of a specific image feature extraction method beforehand [START_REF] Aviles | Sensorless force estimation using a neuro-vision-based approach for robotic-assisted surgery[END_REF][START_REF] Aviles | A recurrent neural network approach for 3d vision-based force estimation[END_REF]. The consequence of this two-stage process to produce synthetic data for training a NN is the ability to classify force ranges with 94% accuracy. We validate this claim on real OCT images collected on five ex vivo porcine eyes using a robotically-controlled needle.

Method

Deep learning has already been proposed to improve existent features of robotic assisted surgeries such as instrument detection and segmentation [START_REF] Pakhomov | Deep residual learning for instrument segmentation in robotic surgery[END_REF]. A key requirement for deep learning methods to work is the high amount of data to train on. Since for the moment, intravitreal injections are manually performed, there is no available information on the force applied by the needle. In this paper, we build a numerical model of the relevant part of the eye to generate images of the deformed sclera under needle-induced forces. The simulations are parametrized to match experimental results and compensate for the lack of real data.

Numerical simulation for fast data generation

Biomechanical model: The deformation of the sclera under needle-induced forces can be approximated by modeling the eye as an elastic material, shaped like a half-sphere subject to specific boundary conditions. Since applied forces and resulting deformations remain small, we choose a linear relationship between strain and stress σ, known as Hooke's law (1):

σ = 2µ + λtr( )I (1)
where λ and µ are the Lamé coefficients that can be determined from the Young's modulus E and Poisson's ratio ν of the material. The linearity of Hooke's law leads us to a simple relation σ = C where C is the linear elasticity constitutive matrix of homogeneous and isotropic materials. Boundary conditions are then added. Dirichlet boundary conditions are used to prevent rigid body motion of the sclera, while a constant pressure is applied on the inner domain boundary to simulate the intraocular pressure (IOP). The IOP plays an important role in the apparent stiffness of the eye and its variability is well studied, as high eye pressure can be an indication for glaucoma. The external force due to IOP is simply given by F = S × P where S is the surface area of the eye given in m 2 and P the IOP in P a. The external forces of the system are formed by the IOP, the force induced by the needle and gravity. Scleral thickness also plays an important role in the deformation of the sclera so it has to be taken into account.

Finite Element simulation: We solve the equation for the constitutive model using a finite element method. To discretize the eyeball, which we consider as nearly spherical, we generated a quadrilateral surface mesh of a half sphere of radius 12 mm, discretized using the Catmull-Clark subdivision method. This quadrilateral surface is then extruded according to the scleral thickness to generate almost regular hexahedra. The deformation is specified by the nodal displacements u and the nodal forces f, according to the following equation:

Ku = f p + f g + f n (2)
The matrix K is the stiffness matrix, and can be computed thanks to the elastic parameters of the material E and ν. E is Young's modulus and is a measure of the stiffness of the material while ν is Poisson's ratio and estimates the compressibility of the material. According to the values in the literature [START_REF] Asejczyk-Widlicka | The elasticity and rigidity of the outer coats of the eye[END_REF] we set E = 0.25 M P a and ν = 0.45. To model the needle pushing on the sclera, we apply a local force f n on a subset of nodes within a small region of interest near the virtual needle tip. f g is the force due to gravity, and f p the force normal to the surface due to the IOP. To compute accurately the deformation of the The simulated images correspond to a 2D cross-section of the entire 3D mesh.

Neural network image classification

Our goal is to create an artificial NN model to predict the force based solely on a single OCT image of the deformation. Using cross sectional OCT images allows one to visualize millimeters beneath the surface of the sclera and therefore provides information about the induced deformation. Using 2D OCT B-scans over 3D scans is preferred here as can be acquired at high-frame rates using low cost hardware [START_REF] Apostolopoulos | Efficient oct volume reconstruction from slitlamp microscopes[END_REF]. Instead of estimating the force as a scalar value, we opt to estimate an interval of forces in order to be less sensitive to exact force values. As such, we set up our inference problem as a classification task where a probability for a class label (i.e. a force interval) is desired.

Given our above simulation model, we can now train our prediction model with virtually an infinite amount of data. In our case, the forces applied in the experiments varied between 0.0 N to 0.08325 N . To ensure quality of the simulations (i.e. small displacements) forces going from 0.0 N to 0.06 N were applied on the Finite Element mesh of the sclera. The objective being the intime withdrawal of the needle to limit scleral damage, we decided to label the forces only using three classes. The classes are such that the first one translates no danger at damaging the sclera, the second means that a considerable force is being applied and the last class triggers the alarm to withdraw the needle. If we set the alarm threshold to 0.03 N , we have the following class ranges (see Fig. 1 for the deformation corresponding to each class):

-Class 0: force values smaller than 0.005 N -Class 1: force values from 0.005 N to 0.03 N -Class 2: force values bigger than 0.03 N Our NN then consists of four layers: an input layer having 5600 neurons, two fully connected hidden layers having 600 neurons each and an output layer of 3 neurons. The output layer is softmax activated. ReLu activation was used as well as a 0.8 dropout factor between layers. For training, we used Gradient descent optimization and the cross-entropy loss function. The network was implemented using Tensorflow and the Keras Python library and was trained from scratch. Validation was performed on unseen simulated images.

Results

In this section we first report the real data acquisition on ex vivo porcine eyes. Then, the virtually generated data set and the following training are presented. And last, we validate the NN on unseen real data.

Experimental set-up

To validate the ability of the NN on real data, we collected pairs of force-OCT images from ex vivo porcine eyes. Five porcine eyes were obtained from the local abattoir and transported in ice. Experiments began within 3 hours of death and were completed between 4 and 10 hours postmortem. During the experiments the eyes were moisturized with water. For the experiments, the eyes were fixed with super glue on a 3D-printed holder to ensure fixed boundary conditions on the lower half of the eyeball. The IOP was measured with a tonometer. For all the eyes the IOP was close to 2 mmHg (i.e. 266 P a). The IOP is low since it decreases dramatically after death [START_REF] Gnay | The importance of measuring intraocular pressure using a tonometer in order to estimate the postmortem interval[END_REF] and no fluid injection was made during the experiments. A medical robot [START_REF] Weber | Instrument flight to the inner ear[END_REF] was used to guide a needle while measuring the applied force. A 22G Fine-Ject needle (0.7 × 30 mm) was mounted at the tip of the end effector of the robot. The needle was placed as close as possible to the B-scan but without intersecting it to avoid shadows on OCT images. The margin between the needle and the B-scan was taken into account in the simulation. The robot moved towards the center of the eye (i.e. normally to the sclera) and forces in the direction of movement were continuously recorded with a sensor. To collect OCT images an imaging system was used to record B-scans over time. The OCT device used 840 nm ± 40 nm wavelength light source, with an A-scan rate of 50 kHz and 12 bits per pixel. The resulting 2D image had a resolution of 512 × 512 pixels, corresponding to roughly 15 × 4 mm. As the images were acquired at a lower rate than the forces, the corresponding forces were averaged over the imaging time for each frame (i.e. two seconds after the starting moment of the line move). Overall, 54 trials were performed. For each eye, one position near the corneal-scleral limbus was chosen and several line moves were performed in the same direction.

Data generation for neural network training

Measurements on ten porcine eyes are made to estimate the thickness of the sclera at the locations where the forces are applied. The values range from 460 µm to 650 µm with a mean value of 593 µm. From the literature [START_REF] Olsen | Porcine sclera: Thickness and surface area[END_REF], the thickness near corneal-scleral limbus is estimated to be between 630 µm and 1030 µm. Hence, we simulated scleras of five different thicknesses : 400 µm, 500 µm, 600 µm, 700 µm and 800 µm. The IOP is the same for all the eyes in the experiments so it is fixed at 266 P a for all simulations. For the smallest thickness, we generated 3200 images of deformed scleras undergoing the stated IOP and forces going from 0.0 N to 0.045 N . For each of the other four thicknesses, we generated 4000 images of deformed scleras where the forces vary from 0.0 N to 0.06 N at different random locations. For each thickness, the simulation took approximatively half an hour. Overall, a data set of 19200 synthetic images was generated within 2 hours and a half (see Fig. 2(b)). The generated images look like images in Fig. 1 below. This images are post-processed with a contour detection algorithm using OpenCV functions and binarized to obtain the images used to train the NN (see Fig. 2(b)).

To train the NN, the data set is split such that 90% of the images are for training and the remaining 10% for validating. Hence, the NN is trained on 17280 images and validated on the other 1920 images. In Fig. 2

Tests on unseen data

The aim of our work is to classify force ranges from real OCT images of the deformed scleras. All the OCT images obtained during the experiments (see Fig. 3(a)) are blurred and thresholded to obtain similar images to the synthetic ones (see Fig. 3(b) and Fig. 2(b)). Now that the OCT images look like the simulated ones they can be given as input to the NN for predictions. For each OCT acquisition, the force measured by the robot is converted into a class label (going from 0 to 2) and is considered as the ground truth (target class). The performance of the classification is given in the confusion matrix in table 1. Each raw of the table gives the instances in a target class and each column the instances in a predicted class. For each class, the correct decisions are highlighted in red. Overall the accuracy of the classifier is 94%.

A force class with label 0 is understood as a very small force meaning that there is no risk of damaging the sclera. On the contrary, a force class with label 1 means the sclera is being considerably deformed (forces in this range go from 0.005 N to 0.03 N ) and a force class with label 2 means that the sclera is potentially being damaged and a withdrawal of the needle is advised. With the eyes used during the experiments, the lowest precision of the NN was found for target class 0 (71%). For target class 2, the precision is of 100%. It is not dramatic if forces of range 0 and 1 are misclassified since the risk of damaging the sclera with forces so close to the alarm threshold is almost null (assuming the threshold is correctly set). On the other hand it is essential that the forces of range 2 are predicted correctly.

Conclusion and Discussion

In this paper, we have proposed a method aimed at improving the safety features of upcoming robotized intravitreal injections. Our vision-based method combining numerical simulation and Neural Networks, predicts the range of the force applied by a needle using only 2D images of the scleral deformation with high accuracy. It indicates in real-time the need to withdraw the needle as soon as a certain alarm threshold is reached. To avoid the need of large training data sets, the NN is trained on synthetic images from a simulated deformed sclera.

It is worth mentioning that more complex scenarios could be simulated such as different eye sizes, variable needle insertion angles and different intraocular pressures. We also propose to improve the simulated images to better match real surgical scenarios. In particular it seems important to add the shadows induced by the needle in the OCT image. From an imaging point of view, we know that the predictions of the NN are very sensitive to the image framing and scaling. To address this issue, we plan to randomly crop each simulated image, and enlarge the data set with it. We might also train a Convolutional Neural Network to reduce sensitivity to image properties and allow for more accurate predictions.

Fig. 1 .

 1 Fig. 1. Real OCT images of the sclera and their corresponding simulated images for the three different force ranges.

Fig. 2 .

 2 Fig. 2. (a) Loss and Accuracy curves for training and validation sets. (b) Fragment of the training data set generated by our numerical simulation.

  (a) the accuracy and loss of the model are shown on both training and validation data sets over each epoch. The validation accuracy curve shows 100% accuracy when classifying unseen synthetic images. This curve is above the training accuracy curve probably because of the high dropout applied during the training.

Fig. 3 .

 3 Fig. 3. (a) Unprocessed OCT image. (b) Blurred and thresholded OCT image.

Table 1 .

 1 Confusion matrix

	Target	Prediction Class 0 Class 1 Class 2 Precision
	Class 0	5	2	0	0.71
	Class 1	0	14	1	0.93
	Class 2	0	0	32	1.00
	Recall	1.00	0.88	0.97