J Philidet

C Gissinger

F Lignières

Magnetohydrodynamics of stably-stratified conducting regions in planets and stars

Magnetohydrodynamics of stably stratified regions in planets and stars

Introduction

Understanding the origin of the diversity of stellar magnetic fields is a research topic in which important advances have been recently made and much more are expected soon. This topic is pushed by new observations of groundbased instruments, in particular spectropolarimeters, as well as observations from space (with Corot and Kepler for instance). It also benefits from the results of simulations which, thanks to the numerical resources now available and the progress of astrophysical modeling, allow us to study the complex physics of stellar interiors taking into account many effects such as rotation, magnetism, turbulence. . . Until recently, only the Sun was accessible to us; now, stellar physics has become a privileged laboratory for deepening our knowledge of different aspects of fundamental physics. It is one of the few themes in astrophysics in which advanced theories can be directly confronted with observational constraints. However, although it is accepted that magnetism can play a major role in most stages of a star life, its effects are mainly ignored in stellar evolution codes (in particular unidimensional codes in which most physical effects are simply parametrised). Taking into account these magnetic effects is one of the major current issues in astrophysics.

Stellar magnetic fields can be classified into two broad categories. In the first category, magnetic fields are generated and maintained by convective motions in deep layers, through what is referred to as dynamo effect. There are countless different mechanisms through which dynamo can give rise either to stationnary fields, large-scale dipoles like those observed for rapidly rotating low-mass, M type stars, or to time-dependent fields. The periodic variation of the solar field has been observed since Galileo through the counting of solar dark spots. For these stars having masses less than or equal to those of the Sun, dynamo mechanisms seem to be very efficient in their thick convective zone, an idea supported by many studies on planetary magnetism [START_REF] Roberts | On the genesis of the Earth's magnetism[END_REF]King, 2013, Petitdemange, 2018]. In the second category, on the other hand, massive stars harbor strong dipolar fields with a stable configuration, commonly referred to as "fossil fields", for it has been established that it results from the formation of the star. Unlike dynamo-driven fields, the origin of magnetism in stars possessing a thick radiative zone located on the top is poorly explored by theoretical studies, and understanding it remains a major challenge.

Thanks to the improved sensitivity of a new generation of spectropolarimeters, the existence of a dichotomy in the magnetic field amplitudes for massive and intermediate mass stars is now clear [START_REF] Lignières | The dichotomy between strong and ultra-weak magnetic fields among intermediate-mass stars[END_REF]. Indeed, the observations show a large gap between strong fields (fossil fields with amplitudes above 300 G) and fields of ultra-weak, sub-gauss magnitude. This gap, referred to as the "magnetic desert", stretches over two orders of magnitude. Massive and intermediate mass stars have a thick radiative zone above a small convective zone (see [START_REF] Braithwaite | Magnetic fields in non-convective regions of stars[END_REF] for a review). Therefore, to understand these observations, we need to both determine how strong magnetic fields can be maintained on long time scales in their radiative zone, and at the same time study the processes of saturation of dynamos from shear flows in radiative zones, that is in stably stratified environments.

After its exponential growth phase, the dynamo-generated magnetic energy can indeed saturate to a level sufficiently high to affect the flow. It is well known that the magnetic tension significantly reduces the differential rotation. One of the complexities of magnetohydrodynamic (MHD) processes is the strongly nonlinear character of the retroaction, through the Lorentz force, of the magnetic field on the flow that initially amplified it. Differential rotation and magnetism thus appear as intimately linked phenomena. Differential rotation, by itself or by developing hydrodynamic instabilities, can generate dynamo-generated fields, and this magnetism can inhibit these instabilities or even greatly reduce the differential rotation by transporting angular momentum outward. This last process could explain the magnitude of the differential rotation observed by asteroseismology for different stars [START_REF] Deheuvels | Seismic constraints on the radial dependence of the internal rotation profiles of six Kepler subgiants and young red giants[END_REF][START_REF] Triana | The Internal Rotation Profile of the B-type Star KIC 10526294 from Frequency Inversion of its Dipole Gravity Modes[END_REF]. According to [START_REF] Fossat | Asymptotic g modes: Evidence for a rapid rotation of the solar core[END_REF], a strong shear layer could be present in the solar radiative zone.

MHD phenomena are crucial to the understanding of planetary interiors as well. The Earth magnetic field is known to be generated by dynamo effects taking place in its liquid outer core. The study of this region in the scope of MHD may therefore give some insight into the variations of the Earth magnetic field. Those take place on time scales ranging from less than a year to millions of years. In particular, the variations of the length of day on timescales ranging from 1 to 100 years are believed to be caused by geomagnetic fluctuations featuring these typical periods, so called secular variations (SV). However, while being observed for a long time and monitored with an ever increasing accuracy [START_REF] Chulliat | Geomagnetic secular acceleration, jerks, and a localized standing wave at the core surface from 2000 to 2010[END_REF][START_REF] Chulliat | Fast equatorial waves propagating at the top of the Earth's core[END_REF][START_REF] Finlay | Recent geomagnetic secular variation from swarm and ground observatories as estimated in the chaos-6 geomagnetic field model[END_REF], geomagnetic variations remain a fascinating problem as several features of the so-called magnetic jerks are still unclear. In particular, it was recently suggested [START_REF] Buffett | Evidence for mac waves at the top of earth's core and implications for variations in length of day[END_REF] that the presence of a certain type of MHD waves resulting from the interaction between magnetic tension, Archimedes force and Coriolis force in a 140 km thin stably-stratified layer located at the top of the Earth's liquid outer core [START_REF] Braginsky | MAC-Oscillations of the Hidden Ocean of the Core[END_REF] may account for unexplained secular fluctuations of the geomagnetic dipole field. The presence of this layer at the core-mantle boundary (CMB) is believed to be due to the diffusion of light elements from the mantle to the core driven by mere pressure gradients [START_REF] Gubbins | The stratified layer at the core-mantle boundary caused by barodiffusion of oxygen, sulphur and silicon[END_REF].

Among the impacts a magnetic field can have on a flow embedded in it, one is the development of MHD instabilities. Under the conditions existing in a stellar radiative region, three main types of instabilities are likely to occur, which differ by the source of free energy they derive from. Gravitational energy can be released through Parker-like instabilities [START_REF] Parker | The Dynamical State of the Interstellar Gas and Field[END_REF], magnetic free energy through purely magnetic instabilities (like the Tayler instability [START_REF] Tayler | The adiabatic stability of stars containing magnetic fields-I.Toroidal fields[END_REF] or induced by a gradient of toroidal magnetic field lines [START_REF] Acheson | Hydromagnetics of rotating fluids[END_REF]) while shear instabilities (like the magneto-rotational instability MRI or stratorotational instability SRI) depend on the free kinetic energy in the differential rotation [START_REF] Balbus | A powerful local shear instability in weakly magnetized disks. I -Linear analysis. II -Nonlinear evolution[END_REF], 1998[START_REF] Menou | Local Axisymmetric Diffusive Stability of Weakly Magnetized, Differentially Rotating, Stratified Fluids[END_REF] for the MRI and Shalybkov and Rüdiger [2005] for the SRI). Differential rotation can also result from baroclinicity as shown by [START_REF] Rieutord | The dynamics of the radiative envelope of rapidly rotating stars. I. A spherical Boussinesq model[END_REF] who considered a Boussinesq model of a radiative envelope. In the simulations presented in this paper, these hydrodynamic or MHD instabilities are not taken into account as we focus here on the low Rossby number regime in which an axisymmetric steady state configuration exists. We postpone the numerical study of three dimensional structures induced by a strong differential rotation to a future paper. We only argue on the possible development of MHD instabilities by means of analytical developments.

In this paper, we therefore propose a direct, global numerical modeling of stably stratified layers, by taking into account both the rotational influence of neighboring layers and the effect of the magnetic field generated in the conducting convective zone, and use this model to describe both the stably-stratified layer at the Earth CMB and stellar radiative zones. We present the numerical model in section 2; in section 3 and 4, we discuss the steady-state flow configuration yielded by these simulations. We show the influence of rotation and stratification on the geometry of the flow in purely hydrodynamical models in section 3 and we focus on the development of super-rotation induced by the combined influences of magnetism and stratification in section 4. In section 5, we perform a local, linear stability analysis to show that the combined effects of shear, magnetic field and stratification could yield MHD instabilities able to account for geomagnetic secular variations in the planetary context, or that they could affect the stability of dipolar fields observed for massive stars in the stellar context. We argue that this latter analytical result could explain the observed "magnetic desert" [START_REF] Lignières | The dichotomy between strong and ultra-weak magnetic fields among intermediate-mass stars[END_REF].

Numerical model

In the present paper, we aim at modeling stably-stratified layers in the geophysical and astrophsyical contexts. In the planetary context, this corresponds to the 140 km thin shell located at the Core-Mantle Boundary (CMB), between the fluid, convective dynamo region and the solid mantle. In the stellar context, it corresponds to a radiative zone (more precisely a radiative envelope in the case of a massive or intermediate mass star). In both cases, the region that is being modeled is located between two spheres, and the aspect ratio in these two situations is drastically different (the planetary case is in the thin gap limit, whereas the stellar context is in the thick gap limit). It is important to stress out that a 3D modeling of such a flow would be extremely time-consuming when it comes to use parameter values relevant to both regimes. In particular, the Ekman number E, which measures the ratio between the viscous force and the Coriolis force, is close to E ∼ 10 -15 in both cases, whereas classical DNS can only reach E 10 -7 -10 -8 [START_REF] Schaeffer | Turbulent geodynamo simulations: a leap towards Earth's core[END_REF]. Furthermore, while classical models are essential to understand the evolution of magnetic fields on the magnetic dissipation time scale (about 10000 years), they are less suitable to the study of much shorter timescales (less than a century), such as we are here. For these different reasons, we use a more simplified approach, in which inertia (i.e advection terms) is neglected and axisymmetric calculations are performed, as described below.

The stably-stratified layer is therefore modeled as an axisymmetric spherical Couette flow embedded in a dipolar magnetic field. The inner sphere of radius r i in which this magnetic field is generated represents the transition to the convective region in both contexts, and has the same conductivity as the modeled fluid. The outer sphere of radius r o , however, is insulating. In the planetary context, magnetic torques acting on the convective (inner) region make it rotate slightly faster than the outer, insulating sphere, such that a weak differential rotation (low Rossby limit) can develop in the stably stratified layer. In the stellar context, the rotation rate of the deeper convective region can increase during different phases of the star's likfe (see introduction). We model this by applying a slightly different rotation rate to the two shells, which will propagate to the fluid in between through no-slip mechanical boundary conditions. The calculations being performed in the frame of reference of the outer shell, the latter is motionless, whereas the inner shell has a slow rotation (compared to the global rotation of the frame of reference). Finally, a stable stratification is applied through a temperature difference between the inner and outer shells.

We consider the Boussinesq approximation, allowing us to neglect variations of the fluid density except in the buoyancy term, leading to the following dimensionless MHD equations :

                 ∂v ∂t = -(v.∇)v + E∆v -2e z × v -∇P + EΛ P m (∇ × B) × B + E RaΘe r ∂B ∂t = E P m ∆B + ∇ × (v × B) ∂T ∂t = -(v.∇)T + E P r ∆T (1)
where v, P , B and T are the dimensionless velocity, pressure, magnetic field and temperature, Θ is the temperature fluctuation as a function of which we express the density fluctuations in the buoyancy term, e z the unit vector representing the rotation axis and e r the local radial unit vector. Time scales have been normalized by Ω -1 o (rotation rate of the outer shell), spatial variables by the radius r o of the outer sphere, the magnetic field by its amplitude B o at r = r o in the equatorial plane and temperatures by the temperature difference ∆T between the shells. We introduce the following dimensionless numbers : the Ekman number E = ν/Ω o r 2 o , the Elsasser number Λ = B 2 o /µρηΩ o (ratio between the Lorentz force and the Coriolis force), the magnetic Prandtl number P m = ν/η, the (modified) Rayleigh number Ra = αg∆T r o /νΩ o , the thermal Prandtl number P r = ν/κ and the Rossby number Ro = (Ω i -Ω o )/Ω o . Following traditional notations, we denote the momentum diffusivity of the fluid as ν, its magnetic permeability as µ, its mean density as ρ, its magnetic resistivity as η, its thermal expansion coefficient as α, its thermal diffusivity as κ and the gravitational acceleration at the outer shell as g.

In order to carry out the study of this system, we follow what has been done in previous studies [START_REF] Dormy | MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field[END_REF][START_REF] Dormy | A super-rotating shear layer in magnetohydrodynamic spherical Couette flow[END_REF] and consider a weak differential rotation, which corresponds to negligible non-linear terms. Thus, we have linearised Eq. 1 by expressing the magnetic field as B = B o + b (B o being a dipolar magnetic field, whose amplitude is controlled by Λ), the temperature as T = T s + Θ (where T s is the stationary solution to the heat equation without source ∆T s = 0), and by neglecting the advection term in the Navier-Stokes equation. We have then integrated these linearised equations using the PaRoDy code [START_REF] Dormy | MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field[END_REF] and later collaborations), over a time long enough to reach a stationary state. Note that compared to previous studies [START_REF] Dormy | MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field[END_REF][START_REF] Dormy | A super-rotating shear layer in magnetohydrodynamic spherical Couette flow[END_REF], the novelty here lies in the introduction of stable stratification. Otherwise, we have followed their numerical method.

P r N Ω 2 ≤ E 2/3 P r N Ω 2 > 1 E = 10 -5 , P r = 1 (a) Q = 10 -4 (b) Q = 0.15 (c) Q = 0.3 (d) Q = 1.5 (e) Q = 3 E = 10 -6 , P r = 0.05 (f) Q = 8 • 10 -5 (g) Q = 0.015 (h) Q = 0.15 (i) Q = 1.5 (j) Q = 3 E = 10 -5 , P r = 0.05 (k) Q = 5 • 10 -4 (l) Q = 0.1 (m) Q = 1 (n) Q = 10 (o) Q = 20 E = 10 -6 , P r = 0.05 (p) Q = 10 -4 (q) Q = 0.1 (r) Q = 1 (s) Q = 10 (t) Q = 20 ---------------------------------→ P r N Ω 2 ≡ P r E Ra r o H ≡ Q
3 Hydrodynamical results and symmetry of the flow

We will discuss in detail the full MHD results in the next section; however, it is instructive to first describe the effect of a stable stratification on a spherical Couette flow in the purely hydrodynamical regime, that is with Λ = 0. Previous studies [Barcilon and Pedlosky, 1967a,b] investigated such an effect in a cylindrical geometry, and found that although the Rayleigh number and the Ekman number are independent control parameters, the structure of the flow is merely controlled by the value

Q ≡ P r N Ωo 2
where N is the Brünt-Väisälä frequency. Using the dimensionless parameters introduced above, we can write

Q = P r × E × Ra × r o r o -r i
and by estimating N 2 = αg∆T /(r o -r i ). In spherical shells, [START_REF] Garaud | On rotationally driven meridional flows in stars[END_REF] has studied the nonlinear problem by combining analytical and numerical developments, assuming axial symmetry of the problem and within some limits that are not considered here. [START_REF] Garaud | On rotationally driven meridional flows in stars[END_REF] has also noticed the important role of Q. We found a similar dependance with Q in our data. More precisely, 3 different regimes can be observed depending on the value of Q, as illustrated by Fig. 1 [ Barcilon and Pedlosky, 1967a,b]. If Q E 2/3 (left panels in Fig. 1), the flow is driven by Coriolis forces and exhibits a cylindrical geometry. On the other hand, if Q 1 (right panels in fig 1), the flow is driven by buoyancy, and a spherical symmetry is observed for the flow. In between, neither one dominates the other. To draw fig 1, we have varied the value of Q by changing E, P r and the aspect ratio. It is important to note that additional simulations (not shown) have been run with different parameters and also show that Q is the key parameter.

In the solar case, [START_REF] Spiegel | The solar tachocline[END_REF] have also noticed that Q is an important parameter in their study on the thickness of the tachocline. Due to the smallness of P r (less than 10 -5 according to Dormy, E. and Soward, A. M.

[2007] and [START_REF] Thomas | Sunspots and starspots[END_REF]), Q would be of the order of 1 in the solar radiative zone. Observations of stellar rotation periods show that Ω can vary by several orders of magnitude [START_REF] Donati | Magnetic Fields of Nondegenerate Stars[END_REF] and the different regimes highlighted here are all relevant in the stellar context.

For the Earth case, [START_REF] Roberts | On the genesis of the Earth's magnetism[END_REF] have provided an estimate for the Prandtl number at P r = 0.1. However, in the outer core, the kinematic viscosity and the thermal diffusivity are so low compared to inertia that turbulence can significantly increase their effective values (turbulent diffusivities) such that P r may approach unity. [START_REF] Buffett | Evidence for mac waves at the top of earth's core and implications for variations in length of day[END_REF] have argued on the possibility that the buoyancy frequency N = αg∆T /H may be comparable to the Earth rotation rate Ω o . Consequently, the Rayleigh number would be linked to the Ekman number by the relationship Ra E ∼ (r o /H). In this case, the Earth would lie in the intermediate regime (0.1 Q Earth 1), meaning that the rotation profile is unlikely to present a simple geometry ; furthermore, the Elsasser number for the Earth has a value close to Λ Earth ∼ 1, which means that the effects of the magnetic field are not negligible and can drastically change this geometry. (How about stars ?)

4 Super-Rotation induced by the combined influences of magnetism and stratification

We carry out the same direct numerical simulations (DNS) as presented in the previous section, but this time we add a dipolar magnetic field (thus flows are no longer purely hydrodynamic). For the reasons we have developed above, non-linear terms in the set of equation 1 have been dropped altogether. We first show in the left panel of Fig. 2 the rotation profile and corresponding azimuthal magnetic field for a flow where no stratification is applied. Note that the linearity of the equations implies that only the relative values for different regions make sense.

As illustrated by Fig. 2-a, the structure of the flow is quite different from the hydrodynamic case, in which the rotation rate shows values halfway between that of inner and outer spheres. Here, the magnetic field lines corotate with the conducting inner core, and the magnetic torque naturally produces a co-rotation of the majority of the flow with the inner core, except at the outer boundary layer where viscosity effects are no longer negligible, and the no-slip condition makes the fluid corotate with the outer boundary. However, Fig. 2-b shows that for sufficiently large magnetic field, part of the fluid located in the equatorial region rotates faster than the inner sphere. This phenomenon is commonly referred to as super-rotation [START_REF] Dormy | A super-rotating shear layer in magnetohydrodynamic spherical Couette flow[END_REF] and has been observed in both numerical simulations [START_REF] Hollerbach | Spherical Couette flow in a dipolar magnetic field[END_REF] and experiments [START_REF] Cardin | The influence of toroidal magnetic field on thermal convection in the core[END_REF] of magnetized spherical Couette flow. The mechanism involved is the following : since the magnetic field lines corotate with the conducting inner core, the bulk flow will therefore also corotate with it. So does the inner boundary layer, because the current can recirculate inside the conducting inner core there ; however, the insulating outer boundary layer can not, which means that the current cannot align with the dipolar magnetic field lines. It is therefore subjected to azimuthal Lorentz forces inducing an extra azimuthal velocity, that can be so high as to overcome the rotation rate of the core. This extra rotation rate propagates to the entire field line, thus explaining the characteristic crescent shape of the super-rotating regions.

Unfortunately, this fascinating phenomenon is known to vanish in presence of a strong global rotation [START_REF] Dormy | A super-rotating shear layer in magnetohydrodynamic spherical Couette flow[END_REF]. As shown in the following, it can however be restored in the presence of sufficiently large stable stratification, allowing super-rotation to play a role in the magnetostrophic regime relevant to the planetary context as well as the stellar context (Λ ∼ 1, R o 1 and E 1). Fig. 2-c shows rotation profiles at the equator for Λ = 1 and different values of the Ekman number E. Regardless of the value of E, super-rotation increases as Ra increases and the super-rotating region migrates outward. As a result, we observe the emergence of a localized shear layer close to the CMB when the stratification is sufficiently high. Note that the position of the super-rotating region on the equator remains sensibly the same when Ra is increased, while its magnitude becomes larger.

Interestingly, the presence of stratification also strongly modifies the behavior of super-rotation in the limit of small Ekman number. Here we take the Earth context as an example : Fig. 3-left shows the evolution of the rotation profiles with the global rotation for a given value of Q and Λ, which have been both fixed to their value expected for the Earth outer core, Q ∼ 2 and Λ ∼ 1.

As the Ekman number is decreased, the super-rotating region is enhanced, suggesting a possible magnetostrophic super-rotation close to the CMB. Precise values of P r and Ra in the thin stably-stratified region of the Earth's outer core are still difficult to constrain, making very conjectural any prediction on the exact structure of such a shear layer. Still, the right panel of Fig. 3 strengthens the idea that a super-rotating region similar to the one described here may be generated in this layer. Indeed, it shows that for sufficiently large stratification, the magnitude of the super-rotation seems to tend to a constant value in the magnetostrophic regime, independent of both the Ekman number and the dimensionless number Q.

These numerical simulations therefore suggest that the stably-stratified layer in the Earth outer core might be characterized by a super-rotating region in which the fluid could rotate up to 30% faster than the inner core. Because this mechanism is associated with a strong differential rotation, it provides a simple source for the generation of the MAC waves generally invoked to account for some of the secular variations of the geomagnetic field. Below, we show that in addition to such waves, the stably-stratified layer can also be prone to MHD instabilities evolving on time scales similar to the Earth's magnetic secular variations.

Local description of unstable MHD modes

As pointed out above, the presence of stratification may strongly enhance the differential rotation, and leads to very localized shear through the generation of so-called super-rotation. The purpose of the present section is precisely to investigate the waves and instabilities prone to be generated in the stablystratified layer whereas previous linear studies have ignored either the effects of differential rotation [START_REF] Buffett | Evidence for mac waves at the top of earth's core and implications for variations in length of day[END_REF] or the effects of stratification [START_REF] Petitdemange | Axisymmetric and nonaxisymmetric magnetostrophic MRI modes[END_REF] in the planetary context. In the stellar context, our stability analysis is strongly inspired by [START_REF] Menou | Local Axisymmetric Diffusive Stability of Weakly Magnetized, Differentially Rotating, Stratified Fluids[END_REF]. By using recent observational constraints and a simpler linear analysis, we argue on the generation of unstable modes and their possible influences.

Because of the force balance achieved between the Lorentz, Coriolis and Archimedes forces, one expects so-called MAC waves to propagate in the stablystratified layer, in addition to Alfvèn waves, purely hydrodynamic waves or torsionnal oscillations. Furthermore, in the presence of a weak magnetic field, a differentially rotating, conducting fluid can enter an unstable regime and develop what is refered to as the MagnetoRotationnal Instability (MRI). It is essential, for instance, to the understanding of the dynamics of accretion disks owing to the presence of a strong shear in these objects [START_REF] Balbus | Instability, turbulence, and enhanced transport in accretion disks[END_REF]]. Although rotation properties of planetary interiors obviously differ from those of an accretion disk, it was also shown to be impacted by this instability [START_REF] Petitdemange | Magnetostrophic MRI in the Earth's outer core[END_REF]. In this section, we study the possibility to see the MRI arise, and especially the impact of stratification and differential rotation on this instability, in two different contexts : planetary CMB and stellar radiative envelopes. In particular, this study stands out in that it extends earlier studies [START_REF] Petitdemange | Magnetostrophic MRI in the Earth's outer core[END_REF] to account for the influence of stratification, which had not been accounted for before. Linear calculations show that stratification increases the threshold of the centrifugal instability. However, stable stratification acting on a differentially rotating fluid allows the development of the so-called StratoRotational Instability (SRI), which amplifies non-axisymmetric modes [START_REF] Molemaker | Instability and equilibration of centrifugally stable stratified taylor-couette flow[END_REF]. It was studied numerically [START_REF] Shalybkov | Non-axisymmetric instability of density-stratified taylor-couette flow[END_REF] and experimentally [START_REF] Bars | Experimental analysis of the stratorotational instability in a cylindrical couette flow[END_REF]. The SRI has no effect on axisymmetric modes.

We follow the linearisation of Eq. 1 we have used in section 3, but this time we restrain ourselves to a local analysis and use the WKB approximation to extract a dispersion relation. We place ourselves near the equator where the magnetic field has no radial component. Ignoring the azimuthal component of the background magnetic field is not an important restriction as it is known that the stability of axisymmetric MRI modes are not influenced by this component [START_REF] Balbus | Instability, turbulence, and enhanced transport in accretion disks[END_REF]]. In addition, this azimuthal component has to match zero at r = r o as the outer shell is insulating in our model. Consequently we are left with the vertical component of B 0 only. Consequently, the state around which we perform the linearisation is described by an azimuthal velocity field V 0 = Ω(s, z)se φ (with a rotation rate whose gradient we leave free), a uniform magnetic field B 0 e z , and the stationary temperature T s (r) described in section 2. Note that in the Earth case, the thermal, momentum and magnetic diffusivities follow κ ≤ ν η, so that we may drop the thermal diffusivity κ (estce que c'est le cas dans les étoiles d'ailleurs ?). We find the following dispersion relation : a 4 σ 4 + a 3 σ 3 + a 2 σ 2 + a 1 σ + a 0 = 0 with :

a 4 = k 2 k 2 z a 3 = 2 k 4 k 2 z E η + RaE k 2 a 2 = k 2 k 2 z (E 2 η k 4 + 2Λ m k 2 z ) + 2 RaEE η -a a 1 = 2E η Λ m k 4 + E 2 η E Rak 2 + RaE k 2 Λ m k 2 z -2aE η k 2 a 0 = Λ 2 m k 2 z k 2 + RaEE η Λ m k 2 z -aE 2 η k 4 -Λ m k 2 z (a + 4) (2)
where

a ≡ 1 s 3 Ω 2 k s k z ∂ z (s 4 Ω 2 ) -∂ s (s 4 Ω 2 ) (3)
and Λ m = EΛ/P m , E η = E/P m , k s and k z are the nondimensionnalised radial and vertical component of the wavevector, and σ the nondimensionnalised complex growth rate of the spatial mode k. This dispersion relation is compatible with the one obtained by [START_REF] Menou | Local Axisymmetric Diffusive Stability of Weakly Magnetized, Differentially Rotating, Stratified Fluids[END_REF], where both their entropy dependent terms and the thermal diffusivity κ are neglected here. Note that, in the case of a purely radial differential rotation (as is the case locally near the equator), the parameter which we denote as a simplifies to a = -4 + 2Ro where the definition of the shear rate Ro is based on the angular velocity jump δΩ over a characteristic length scale δs: Ro = (s/Ω)δΩ/(δs). We warn the reader that despite the fact they both measure differential rotation, the shear rate Ro (which is a local quantity) is different from the Rossby number Ro introduced in section 2 (which is a global quantity).

The Earth case

The relevant force balance in the Earth core is the so-called MagnetoStrophic regime, which allows for the development of a slightly modified version of the classical MRI. We refer to this as the MagnetoStrophic-MRI, or MS-MRI [START_REF] Petitdemange | Magnetostrophic MRI in the Earth's outer core[END_REF]. In particular, the differential rotation and the stable stratification only play minor roles in this force balance. In this regime, a simpler (second order) dispersion relation can be obtained. Note that a similar dispersion relation can also be derived outside the magnetostrophic regime for a stably-stratified fluid in the special case where the eddy-diffusivities η and κ are equal :

4Ω 2 (σ + ηk 2 ) 2 + k 2 k 2 z (k z V Az ) 4 + αg∆T (k z V Az ) 2 + 2Ω(k z V Az ) 2 s dΩ ds = 0 (4)
where V Az is the Alfvèn speed V Az = B z / √ µρ. The most unstable mode corresponds to a purely vertical one (k s = 0), and its growth rate is :

σ r Ω = |Ro | - H r o RaE Λ/2 1 + √ 1 + Λ 2 (5)
This dispersion relation allows us to study the stability of a given mode, we do in the following. In the Earth case, among the dimensionless parameters which play a role in the dispersion relation, several are fairly well constrained : we take E ∼ 10 -15 , Λ ∼ 1 and P m ∼ 10 -5 (see Dormy, E. and Soward, A. M. [2007]). Both differential rotation and stratification however are poorly constrained, and as a result we keep Ro and Ra free in the following. In Fig. 4, we have selected one spatial mode and studied both its stability and its oscillatory behaviour when differential rotation and stratification are changed.

When there is no stratification (left panel in Fig. 4), two different regimes can be extracted : if the shear rate is very low, the mode oscillates with a pulsation that is independent on Ro . At higher Ro , however, the mode becomes unstable, and no longer oscillating. This is not surprising, as the differential rotation acts as a reservoir of energy for the MRI. Adding stable stratification (right panel in Fig. 4) decreases the growth rate of the instability in the unstable regime, especially for strong differential rotation. The purely oscillating, low-Ro regime is replaced by a purely evanescent, non oscillating regime. Here, we see what equation 5 already taught us, that is stable stratification plays a stabilizing role for MS-MRI modes, in the sense that it suppress stable oscillations.

This local study also allows us to search for the most unstable spatial mode for a given set of dimensionless parameters. Figure 5 shows the evolution of the growth rate of the most unstable mode with the shear rate. In the left panel, we present the non-stratified case : the blue curve shows the numerical computation, with a low-Ro regime where σ r varies linearly with the shear rate, and a high-Ro regime where σ r varies with the square-root of the shear rate. We show the prediction from the magnetostrophic regime in orange, which we predictably retrieve for very low shear rates.

We show the stratified regime in the right panel. As we previously mentionned, stable stratification prevents the development of unstable modes at low shear rates. However, it leaves it unscathed when the differential rotation is stronger.

On a final note for the Earth case, we point out that for a shear rate of the order of several percents (which is the typical value obtained if the shear is Figure 4: Left : Evolution of the growth rate σ of a given spatial mode with the shear rate Ro (E = 10 -15 , P m = 10 -5 , Λ = 1 and R a = 0). The blue curve shows the real part of σ when it is positive (i.e when the mode is unstable), the red curve shows the imaginary part of σ when it is non-vanishing. Right : Same as left panel, but with R a = 10 12 . The graph solid blue line shows the real part of σ when it is positive, and the dashed blue line shows its opposite when σ r is negative.

Figure 5: Left : Evolution of the growth rate σ of the most unstable mode with the shear rate Ro (E = 10 -6 , P m = 10 -5 , Λ = 1 and R a = 0). The orange curve shows the analytical solution that can be derived in the manetostrophic regime (equation 5). The blue curve shows the solution computed numerically.Right : Same as in the left panel, but with R a = 10 12 . due to super-rotation), and if we consider that stratification is not high enough to suppress the instability, the numerical study yields σ r ∼ 10 -3 in units of Ω o , corresponding to a typical growth time τ ∼ 3 years, in agreement with the observed typical time scale for the Earth magnetic jerks.

The stellar case

The same dispersion relation can be used regardless of the context. However, the observational constraints are very different in the stellar case and in the Earth case, and we will therefore perform the analysis in the stellar case in a different fashion. In particular, following the stellar community, we prefer to describe the stratification with the Brunt-Väisälä frequency N rather than with the modified Rayleigh number, the two being linked via :

N Ω 2 = r o r o -r i RaE (6)
As an example, we consider here the case of a star having just begun to climb the red giant branch (RGB). Due to the contraction of the core and the swelling of the enveloppe, we consider that the differential rotation is quite strong. Observational constraints indicate that the angular momentum is then redistributed throughout the star (see for instance Deheuvels 2014), thus damping the degree of differential rotation ; however, we place ourselves before the redistribution mechanism kicks in, and for illustrative purposes take a ∼ 1 (which corresponds, for instance, to a keplerian rotation field). Furthermore, asteroseismic constraints also allow to probe the Brunt-Väisälä frequency in the interior of red giants, and give an estimate for the stratification. In particular, for a rapidly rotating red giant, we may take N ∼ Ω. Taking the Ekman number as E ∼ 10 -15 (see Rincon 2019 again) then gives us an estimate for the Rayleigh number. Note to self : this entire paragraph must be burnt at all cost, and replaced by something that might possibly make sense instead.

We present in Figure 6 the results yielded by the dispersion relation 2 in this stellar case. The first three panels show the phase diagram of the system stability in the (P m, Λ) plane for different amplitudes of the stratification. The non-stratified case (top-left panel in figure 6) yields a phase diagram that is well-known : there is a P m min below which the instability does not arise, and above that critical value one can see that the system is only unstable if Λ is in a certain window, whose width increases with P m. The other panels show that, as we previously encountered with super-rotation, stable stratification has a stabilizing effect on the flow, in the sense that the unstable domain shrinks when Ra is increased. Note that in order to find the most unstable modes, we have arbitrarily restrained the search to spatial modes comprised between 10 -7 and 10 7 in units of the radius gap H between the two shells, for both the radial and vertical component. This was done with the understanding that wave vectors whose norm is below 10 -7 are bigger than the system itself, rendering this local analysis invalid, and that wave vectors whose norm is above 10 7 are killed by dissipative effects. Therefore, while the dispersion relation always formally produce unstable modes regardless of the value taken by the dimensionless parameters, there is a domain, drawn in white in Figure 6, in which the only unstable spatial modes are not physically meaningful.

In the last panel of Figure 6 we present a vertical cut of the other panels at P m = 10 -2 , which corresponds to the value of the magnetic Prandtl number at the top of the Sun radiative core for instance. One can clearly see that the instability window shrinks with increasing stratification. The typical growth rate of the instability within the window remains, however, independent from the stratification.

The results presented above support the idea that, in the parameter range relevant to the context of stellar radiative zones, the MRI is indeed likely to arise, but only if the magnetic field lies in a certain range. Therefore, the MRI may provide a viable mechanism to explain magnetic desert among intermediate mass stars. Indeed, if the magnetic field harbored by the star during the early phases of its life has a high enough amplitude, it will prevent this instability from developing, and the star will continue to feature this fossil field. Otherwise, the MRI can develop, and its saturation will destroy the magnetic field that generated it. Alone would survive the magnetic field generated by dynamo effect in the convective zone of the star; given the small size of said region, the surviving field will be much lower than the original field, and may not be subjected to the MRI. As a result, we propose the following picture of intermediate mass star magnetism : A p /B p magnetism is a result of fossil fields strong enough to prevent the MRI from developing; and in Vega-like magnetism the fossil field has been destroyed by MRI, thus letting dynamo-generated fields dominate, which are too small to trigger the MRI. Note that, in regards to the magnetic instability windows shown in Fig. 6, its width in terms of Λ is twice that in terms of B o , since Λ ∝ B 2 o . There is another type of stars in which differential rotation is of primary importance, and those are red giant stars. It is well known, indeed, that these stars feature a strong, localised shear layer near their H-burning shell [START_REF] Deheuvels | Seismic constraints on the radial dependence of the internal rotation profiles of six Kepler subgiants and young red giants[END_REF], which is however much weaker than predicted from theory or simulations (the shear rate Ro inferred from observations is of order unity, while the predicted value would be of order 10 2-3 ). The MRI is likely to play a crucial role in the dynamics of such region; if that is the case, the saturation of this instability would, as we mentionned in the case of main-sequence, intermediate mass stars, destroy the flow configuration that allowed it to develop, thus smoothing differential rotation. Indeed, since the convective region in giant stars is much deeper than in intermediate mass stars, dynamo-driven magnetic fields may in this case resist the development of the MRI, and the flow would then be affected not through its magnetic field, but through its differential rotation. Therefore, it is possible that the MRI may play a role in the angular momentum transport needed to account for the observed rotation rate difference between the core and the enveloppe.

Summary and perspectives

In this study, we have modeled a stably-stratified region with an axisymmetric spherical Couette flow embedded in a dipolar magnetic field. We have split this study two ways, by applying this model to two different astrophysical objects : the thin, stably-stratified layer at the Earth CMB; and stellar radiative zones.

In the Earth context, the most striking feature is the generation of superrotation in the limit E 1, Λ ∼ 1. Combined with MAC waves which are known to transport angular momentum in such a stably-stratified layer [START_REF] Buffett | Evidence for mac waves at the top of earth's core and implications for variations in length of day[END_REF], the magnetostrophic super-rotation can exert strong shear at the CMB. This new source of shear might be crucial for the understanding of some of the behavior of the Earth's magnetic field, such as the variations in the length of day, or the geomagnetic jerks. Indeed, our local analysis shows that the presence of strong enough shear at the Earth's CMB might trigger MHD instabilities. We showed that the stratification has little direct impact on the growth rate of such instabilities, which arises as long as Ra remains reasonably low. This suggests that such a flow is likely to be unstable as long as the amplitude R o of the differential rotation due to super-rotation remains close to the rough estimate ∼ 10 -2 . The stratification enables the presence of such a localized shear by restoring super-rotation (in the magnetostrophic regime) close to the CMB. In addition, the localization of super-rotation at the equator may explain why several magnetic features are mainly observed at very low latitudes (see for instance [START_REF] Chulliat | Geomagnetic secular acceleration, jerks, and a localized standing wave at the core surface from 2000 to 2010[END_REF] and [START_REF] Finlay | Recent geomagnetic secular variation from swarm and ground observatories as estimated in the chaos-6 geomagnetic field model[END_REF]). Furthermore, the observed variations in the length of day might be caused by an acceleration of the Earth's mantle due to the non-linear part of the instability. Finally, note that a non-axisymmetric model would be required to explain periodic local magnetic jerks [START_REF] Petitdemange | Axisymmetric and nonaxisymmetric magnetostrophic MRI modes[END_REF]. Further investigations on this matter would require a more in-depth study of non-axisymmetric modes in a stratified flow.

In the stellar context, we focused not so much on the possibility to see super-rotation arise from the interaction between differential rotation and magnetic field (such a shear structure may in this context be subject to large-scale flows, like meridional recirculation, which we have not taken into account) as we have on the potential instabilities which might dramatically affect the flow, in particular the MRI. We show here through a local linear analysis that in the parameter range relevant to the context of stellar radiative zones, the MRI is indeed likely to arise. One of the recognisable characteristic of this intability is that it requires the magnetic field amplitude to lie in a certain instability window. We show that stable stratification plays a stabilizing role with respect to the MRI, in the sense that it reduces this unstable window. While it may seem far-fetched to draw quantitative conclusions at this stage, we nonetheless propose a mechanism through which the MRI may explain the magnetic desert observed among intermediate mass stars.

The model we propose here is a simple one, and many aspects of the situation have been overlooked in the process. In particular, we have overly simplified the behaviour of the transition region between the convective interior and the radiative enveloppe at radius r i . While we have imposed a mechanical no-slip boundary condition there, it is well known that, because of inertia, the convective flow in the former overshoots in the latter. The radial velocity entailed by this mechanism thus have on the flow in the stably-stratified layer an impact that may be global; further studies will have to adress this point. Furthermore, it is natural that our axisymmetric study be extended to non-axisymmetric, tridimensionnal simulations, for a lot of instability mechanism in the studied contexts do not arise in axisymmetric environments. Finally, we have focused in this paper on the linear part of the instability, before its saturation. While we have strived to draw conclusions about the impact of the saturated phase on the flow, it will be necessary to study it numerically to confirm said conclusions. In particular, only then may the results obtained be quantitatively compared to observations. Observations are indeed plentiful as far as this study goes. The Kepler mission has in particular given us access to extensive data on giant stars, thanks to which the issue of angular momentum redistribution can be studied in excruciating detail. Furthermore, concerning the intermediate mass stars magnetic desert, new observations of very young stellar objects may help us shed some light on the origin of magnetism in these stars, and in particular may help us to support the validity of the mechanism proposed in this paper. Tridimensionnal simulations are also crucial when it comes to understanding the origin of secular variations (SV) in geomagnetism, which are inherently non-axisymmetric. citer Olsen have investigated on the matter already, but in a parameter range that did not allow for the developement of the MRI. Performing a similar study with much lower Ekman number will allow for the study of the impact of this instability on the geomagnetic SV. The vastness of issues still open on the subject, as well as the innumerable observationnal data available to help solve them, mean that there is still a lot of work to be done.

Figure 1 :

 1 Figure 1: Rotation profile for purely hydrodynamic Couette flows with a weak differential rotation (Ro = 0.05) for different values of E and P r, two different aspect ratios (r i /r o = 0.35 for the first two rows, 0.9 for the last two), as well as increasing values of Q (from left to right). As the amplitude of stratification is increased, the flow undergoes a continuous transition from cylindrical to spherical symmetry. The Earth's stratified layer lies in the intermediate regime, while the case of stellar radiative zones clearly lies in the spherically symmetric, high-Q regime. The solid black lines shown in non-trivial geometries are contours of equal rotation rate.

  Figure 2: On the left (panel a, d) : Azimuthal magnetic field b φ (left) and rotation profile (right) for Λ = 1, R o = 0.01, P r = 0.05 and E = 10 -5 , for two different aspect ratios (r i /r o = 0.9 for the top row, 0.35 for the bottom row). Otherwise, equatorial rotation profile normalized by the differential rotation ∆Ω, shown for different values of Λ at E = 10 -6 , Ra = 2 • 10 5 , P r = 0.05, (panel (b,e)), and for different values of Ra at Λ = 1, P r = 1, E = 10 -7 (panel (c)). Regardless of the aspect ratio (thin gap on the top and thick gap on the bottom), super-rotation increases with both the magnetic field (Λ) and the stratification ( Ra).

Figure 3 :

 3 Figure 3: Left : equatorial rotation profiles for different values of E ( Ra×E = 2, all other parameters same as Fig.1), showing the generation of a strong shear layer at low Ekman number. Right : Maximum amplitude of the equatorial angular velocity Ω max normalized by the applied differential rotation ∆Ω = Ω i -Ω 0 for different values of P r and E × Ra, all other parameters same as Fig.1. Sufficiently large stratification leads to super-rotation independent of E at low Ekman number.

Figure 6 :

 6 Figure 6: Top left : Phase diagram of the flow stability in the (log Λ, log P m ) plane, for E = 10 -15 , Ro = 0.01 and E × Ra = 10 -4 . The color bar shows the growth rate of the most unstable spatial mode (logarithmic scale) when it exists. The white zones on the left hand side of the phase diagram corresponds to the stable domain. Top right : Same with E × Ra = 1. Bottom left : Same with E × Ra = 10 4 . Bottom right : Evolution of σ r with Λ, for the same parameters as the three other panels, P m = 10 -2 and different values of Ra, showing that the unstable Λ-window shrinks as the amplitude of stratification is increased.

Acknowledgments

This study was granted access to the HPC resources ofMesoPSL financed by the Région île-de-France and the project EquipMeso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. Numerical simulations were also carried out at the CINES Occigen computing centers (GENCI project A001046698). LP acknowledges financial support from "Programme National de Physique Stellaire" (PNPS) of CNRS/INSU, France.


D Bouchet

J Scholler

G Blanquer

Y De Wilde

I Izeddin

V Krachmalnicoff
email: valentina.krachmalnicoff@espci.fr

Probing near-field light-matter interactions with single-molecule lifetime imaging

Nanophotonics offers a promising range of applications spanning from the development of efficient solar cells to quantum communications and biosensing. However, the ability to efficiently couple fluorescent emitters with nanostructured materials requires to probe light-matter interactions at subwavelength resolution, which remains experimentally challenging. Here, we introduce an approach to perform super-resolved fluorescence lifetime measurements on samples that are densely labelled with photo-activatable fluorescent molecules. The simultaneous measurement of the position and the decay rate of the molecules provides a direct access to the local density of states (LDOS) at the nanoscale. We experimentally demonstrate the performance of the technique by studying the LDOS variations induced in the near field of a silver nanowire, and we show via a Cramér-Rao analysis that the proposed experimental setup enables a single-molecule localisation precision of 6 nm.

Single fluorescent emitters constitute an excellent probe to access the evanescent near-field of a nanostructure with far-field measurements. Indeed, the advent of super-resolution microscopy in the field of biophotonics has uncapped an unprecedented detail of observation of subcellular structures revealing structural features of tens of nanometres [1][2][3], one order of magnitude below the resolution limit imposed by the diffraction of light. While the main super-resolution approaches are based on fluorescence intensity measurements, there exists a strong interest in developing techniques capable of probing lifetime variations at the nanoscale by associating fluorescence lifetime imaging microscopy (FLIM) with subwavelength spatial information. The far-reaching potential of fluorescence lifetime imaging with nanometre resolution is straightforward not only for biological studies [4,5] but also for nanophotonics applications [6,7], as the lifetime of fluorescent emitters is inversely proportional to the LDOS [8].

In the last few years, different experimental approaches have been proposed to achieve lifetime measurements at the nanoscale. Super-resolution lifetime imaging was first demonstrated in combination with stimulated emissiondepletion (STED) microscopy [9], mostly used for biological applications, and more recently by making use of scanning-probe microscopy to characterise the response of nanostructured plasmonic [10][11][12][13][14] or dielectric [15] materials to light. Despite the contribution of these methods to nanoscale imaging, a wide-field scheme rather than a scanning approach is essential in order to study dynamic phenomena and to reach molecular resolution. Several groups have recently proposed wide-field approaches to obtain super-resolved LDOS measurements. The association of wide-field localisation with a scanning scheme was used to probe lifetime variations induced by periodic structures [16]. Elegant techniques (although arduous to master) were implemented to measure the lifetime of single quantum dots positioned with microfluidic flow control [17,18] or using surface-bound motor proteins [19], allowing one to image LDOS variations induced by plasmonic nanostructures. Other methods based on point accumulation for imaging in nanoscale topography (PAINT) [20,21] and photo-activated localisation microscopy (PALM) [22] need numerical simulations to estimate the LDOS from intensity-based measurements.

In this Letter, we introduce a novel approach that overpasses these limitations and combines lifetime and superresolved spatial information based on stochastic optical reconstruction microscopy (STORM) [3], a stochastic imaging technique widely used in biological imaging [23]. This method allows to map the lifetime τ of stochastically photo-activated single molecules in close vicinity of a densely-labelled nanostructure. It can be readily implemented with a standard microscope and can be applied to biological samples or artificially fabricated nanostructures, either dielectric, metallic, or hybrid metallo/dielectric. Here, we demonstrate the performance of the technique by mapping the LDOS variations induced by a silver nanowire on single molecules located a few nanometers apart. Plasmonic nanowires are an ideal playground to demonstrate the ability of a superresolved technique to measure light-matter interactions on the nanometer range. They induce strong variations of the lifetime of nearby emitters on the nanometre scale, highlighting the large dynamic range in terms of lifetime modification explorable with our technique. Moreover, due to their geometric simplicity, they enable handleable theoretical studies easily comparable to experimental results.

The sample consists of silver nanowires on a glass coverslip, the whole covered with photo-activatable fluorescent molecules, and is illuminated in wide field with a pulsed laser through an oil immersion objective mounted on an inverted microscope. The studied nanowires have a diameter of ∼115 nm and a length of several tens of microns. Their large longitudinal dimension ensures that they weakly radiate to the far-field, therefore strongly limiting the shift in the apparent position of the emitters that has been observed for resonant nanostructures [18,20,24]. The specificity of our method relies on the simultaneous detection of fluorescence photons, through the same microscope objective as the one used for the excitation, on an electron-multiplying charge-coupled device (EM-CCD) camera for super-localisation and on a single-photon avalanche diode (SPAD) coupled to a timecorrelated single-photon counting (TCSPC) system for lifetime measurements (see Supplementary Section 1). The EM-CCD camera records wide-field images of the sample with a field of view of tens of micrometers on the sample plane. In contrast, the SPAD, which is a single-channel detector, is conjugated with the center of the camera image via a 50 µm confocal pinhole and covers an area on the sample plane of ∼1 µm 2 . By setting the excitation and photoactivation laser power so that no more than one molecule is active at a given time on the area conjugated to the SPAD, the decay rate Γ = 1/τ can be properly estimated for each individual molecule and can be associated to its position.

This approach is illustrated in Fig. 1. A single fluorescent molecule is identified on a sequence of wide-field images (Fig. 1a), and the position of this molecule is estimated by fitting a two-dimensional Gaussian function to the measured point spread function (Fig. 1b). At the same time, the detection of a fluorescent molecule appears as a burst on the signal of the SPAD time trace (Fig. 1c). For each SPAD burst, we build the associated decay histogram with a time resolution of 16 ps ( Fig. 1d,e). To estimate the decay rate, the convolution of the instrument response function (IRF) and a decreasing mono-exponential function is fitted to the decay histogram. Based on the time correlation between the events detected by the camera and the SPAD (see Supplementary Section 2), we can associate position and decay rate for a large number of photo-activated molecules detected in a single experiment and obtain the superresolved decay rate map shown in Fig. 1f. This map is reconstructed from simultaneous position and decay rate measurements of 3119 molecules, located in a sample region of 1 µm 2 containing one silver nanowire. The typical localisation precision, calculated via a Cramer-Rao lower bound analysis as explained below, is of the order of 6 nm. Spatial variations of the decay rate are observed well below the diffraction limit, demonstrating the ability of the technique to obtain super-resolved LDOS images in a wide-field optical configuration.

A unique insight as allowed by this new approach is revealed by the study of the density of detected molecules along the center of the nanowire axis (see also Supplementary Section 3). Fig. 2a shows that, on average, twice as many molecules are detected for a distance to the nanowire axis d = ±50 nm than for d = 0 nm. Indeed, the interaction between the excitation field and the nanowire results in a non-uniform excitation intensity distribution, as shown in Fig. 2b by finite-difference time-domain (FDTD) simulations (see Supplementary Section 4). A local enhancement of the excitation intensity is observed on the sides of the nanowire, with a lateral extension of about 20 nm, as well as an extinction of the excitation intensity on the top of the wire. Therefore, the molecules located in the higher excitation intensity regions have a larger probability to be detected, supporting the observed variations of the density of detected molecules. Furthermore, the image is formed by a two-dimensional projection of fluorescent events around a cylindrical nano-object. This also affects the apparent density of detected molecules. In order to get a deeper insight into the observed decay rate variations, we further studied the dependence of the decay rate on the distance d to the nanowire axis (Fig. 2c). Molecules detected far from the nanowire axis (d > 200 nm) show an average value of the decay rate of 0.68 ns -1 with a standard deviation of 0.17 ns -1 . In contrast, the decay rate is higher than 10 ns -1 for many molecules detected at distances d < 60 nm from the nanowire axis. This leads to a decay rate enhancement of a factor 15, only limited by the IRF of the setup (see Supplementary Section 5). This measurement confirms that molecules with the largest decay rates are those attached to the nanowire or in its closest vicinity. We further numerically simulated the enhancement of the decay rate induced by the presence of the nanowire for three orthogonal dipole moment orientations (Fig. 2d). Experimental and numerical results are in good qualitative agreement, supporting the validity of the experimental technique. Different dipole moment orientations can explain the lifetime dispersion observed in the vicinity of the nanowire.

The performance of the proposed method ultimately relies on the precision at which we can estimate both the position and the decay rate of the detected single fluorophores. We can assess a lower bound on these parameters by calculating the Cramér-Rao lower bound [START_REF] Kay | Fundamentals of Statistical Processing, Volume I: Estimation Theory[END_REF] on the standard error of the position and lifetime estimators, respectively noted σ x,y and σ Γ (see Supplementary Sections 6 and 7). Such analysis is standard in localisation microscopy to assess the localisation precision [START_REF] Deschout | [END_REF]27]. Fig. 3a shows the dependence of σ x,y upon the number of fluorescence photons detected by the EM-CCD camera. The fundamental limit (red curve) is set by the shot noise and the finite pixel size as sources of error on the measurement. The instrumental limit (green curve) also accounts for the readout noise of the camera and the noise introduced by the electron multiplying process. The actual limit of our experiment (blue curve) is calculated by considering additional sources of noise such as substrate luminescence. The number of fluorescence photons experimentally detected by the camera from each molecule ranges from 150 to more than 10 4 fluorescence photons with a median value of 1228 photons (3a, bottom). With this value the Cramér-Rao bound for position estimations is 6 nm.

A similar analysis can also be performed for lifetime estimations [28,29]. The fundamental limit on the relative standard error of decay rate estimators σ Γ /Γ is simply given by 1/ √ N , where N is the number of detected photons (Fig. 3b, red curve). We calculated the Cramér-Rao bound for Γ = 0.7 ns -1 (molecules on glass) and Γ = 7 ns -1 (molecules close to the nanowire), which corresponds to a lifetime of 140 ps, comparable to the FWHM of the IRF (240 ps). As expected, σ Γ /Γ deviates from the fundamental limit when the number of measured fluorescence photons is smaller than 1000 due to the influence of background noise. In the experiment, the median value of detected photons is 367 photons (Fig. 3b, bottom). For this value, σ Γ /Γ ranges from 8% to 10% depending on the value of Γ. The Cramér-Rao analysis thus demonstrates that the proposed experimental setup enables state-of-the-art measurements of light-matter interactions with a localisation precision of 6 nm together with a relative error of 10% for lifetime estimations. Future prospects will include accessing the axial position of the detected molecules with the implementation of three-dimensional localisation methods [30][31][32][33]. The technique can notably be adapted to perform three-dimensional imaging using metal-induced energy transfer, as suggested by a recent article reporting three-dimensional localisation for sparsely distributed single molecules [34]. Additionally, by taking advantage of SPAD arrays constituted of several independent channels [35], a field of view of tens of micrometres in the sample plane could be reached, opening a wide range of interesting opportunities for imaging and sensing applications. The acquisition time could also be reduced by actively optimising the number of molecules simultaneously photo-activated on the region conjugated to the SPAD. Multiple simultaneous detections could be treated with an improved detection, fitting, and reconstruction algorithm.

The readiness of the technique to be implemented with a standard microscope suggests a great potential to rapidly expand into a wide variety of applications, ranging from nanophotonics and plasmonics to biophotonics. Topical applications in nanophotonics include the direct characterisation of samples presenting rich LDOS patterns and strongly-confined electromagnetic fields, with concrete perspectives for the study of light localisation in strongly scattering media [36,37]. The technique is not constrained to the fluorphores used in the present realisation, but can be extended to photoactivatable fluorophores of different wavelengths and to DNA-PAINT for an a priori knowledge of the fluorophore position. Thanks to these extensions, it will be possible to characterise the resonant and non-resonant behavior of a nanostructure and to tackle the mislocalisation of resonant fluorophores with fluorescence lifetime measurements. In the field of biophotonics, wide-field FLIM images with nanometre resolution will allow to probe local dynamic phenomena in living cells. We also foresee that, by associating our approach with techniques based on Förster resonance energy transfer (FRET), in cellulo nanoscale imaging of molecule-molecule interactions will soon become within reach.

Probing near-field light-matter interactions with single-molecule lifetime imaging: supplementary material To prepare the sample, we spin-coat a dilute solution of silver nanowires in isopropyl alcohol on a glass coverslip. A microfluidic chamber is then prepared as follows [1]: we cover the sample with a ring made of parafilm, we place two micro-pipettes on opposite sides of the parafilm ring and we cover them with another glass coverslip before heating the sample up to 70 • C in order to melt the parafilm. We let the microfluidic chamber cool down for a few minutes before using the micro-pipettes to inject biotin diluted in a phosphate-buffered saline (PBS) solution at a concentration of 1 g/L. We leave this solution incubate for 2 hours. Then, we inject streptavidin-conjugated fluorescent molecules (Alexa 647) diluted in a PBS solution at a concentration of 0.005 g/L, and we leave this new solution incubate for 2 hours. We add a PBS solution containing a few polystyrene fluorescent beads 100 nm in diameter (Red FluoSpheres, ThermoFisher Scientific) which we use as fiducial markers, and we then fill the chamber with an oxygen-reducing buffer [2]. This buffer is prepared according to the following protocol [3]: we use a PBS solution in which we dilute dextrose (100 mg/mL), cysteamine (3.86 mg/mL), glucose oxidase (0.5 mg/mL) and catalase (1.18 µL of an aqueous solution concentrated at 20-50 mg/mL).

B. Optical setup

Before the experiment, we select an area on the sample in which a silver nanowire can be identified by basic transmission imaging, thus ensuring that only one nanowire is present in the detection volume. Then, we place the area of interest in the middle of the field of view of the camera by using a piezoelectric stage (PXY 200SG, Piezosystem Jena). Photo-activatable molecules (Alexa Fluor 647) are excited by a pulsed laser diode emitting at λ = 640 nm (LDH Series P-C-640B, Pico-Quant) at a repetition rate of 80 MHz. The intensity incident on the sample averaged over a repetition period is 10 µW/µm 2 . The laser polarisation is set perpendicular to the nanowire axis in order to minimise the backscattering of the laser light by the nanowire. The molecules are photo-activated with a laser diode emitting at λ = 405 nm (LDH Series P-C-640B, Picoquant). During the acquisition, the density of activat-able molecules decreases in time since several molecules are photobleached by the excitation laser. To compensate for this effect, we progressively turn on the photoactivation laser, with an average intensity on the sample up to 50 nW/µm 2 . A third laser (Fianium SC450) filtered at λ = 568 nm is required for the excitation of fiducial markers that are used for real-time drift correction. These three lasers illuminate the sample through an oil immersion objective (UPLSAPO 100XO, NA=1.4, Olympus) mounted on an inverted microscope (Fig. S1). Wide-field illumination over an area of approximately 200 µm 2 is achieved by placing a lens (f = 300 mm) before the objective. Fluorescence from the sample is then collected by the objective and filtered by a dichroic mirror as well as two long-pass filters. Then, a 50:50 beamsplitter splits the signal towards two paths. On the first path, fluorescence photons are directed towards an EM-CCD camera (iXon 897, Andor). On the second path, a SPAD (PDM-R, Micro Photon Devices [4]) is connected to a time-correlated single-photon counting (TCSPC) system (HydraHarp400, Picoquant). together with the photo-activation laser (λ = 405 nm) and the laser used for sample stabilisation (λ = 568 nm), illuminate the sample via a high numerical aperture oil objective (NA=1.4). A lens (f=300 mm) is located on the excitation path to ensure wide-field illumination. Fluorescence from the sample is filtered with a dichroic mirror (DM) and passes through a tube lens (TL). A 50:50 beamsplitter (BS) splits the light towards an EM-CCD camera and a SPAD. The sample under study contains photo-activated single molecules in the near-field of a silver nanowire (NW).

C. Drift correction

To determine and correct the drift in the sample plane, we estimate the position of a fiducial marker from the wide-field images acquired by the camera and we use a feedback loop to maintain the marker at a fixed position. Every 5 s, the drift is estimated by fitting a two-dimensional Gaussian function to the image of the marker. A feedback signal is then applied on the piezoelectric stage (PXY 200SG, Piezosystem Jena) controlling the in-plane position of the sample in order to compensate for the drift.

In order to estimate the drift of the sample in the axial direction with respect to the focal plane, we analyse images of the fiducial marker accumulated over several seconds. The defocus-correction system is based on a real-time maximisation of the power spectral density of the measured images, and the axial position of the objective with respect to the sample is corrected in real time with a piezoelectric positioning system (MIPOS 20SG, Piezosystem Jena) located between the objective and the microscope turret.

POSITION AND DECAY RATE ASSOCIATION A. Position estimations

The EM-CCD camera acquires 31 frames per second with an acquisition time of 30 ms per frame. The full sequence of wide-field images saved by the camera (over a subset of 13×13 pixels, pixel size = 160 nm) is imported by ImageJ [5] and the positions of the photo-activated molecules are estimated using ThunderSTORM [6]. First of all, each frame is filtered using a wavelet filter, as proposed by Izeddin et al. [7]. For each frame, approximate localisation of the molecules is then performed by applying a threshold that depends on the signal-to-noise ratio of the camera data. For this acquisition, we set it to 2.7 times the standard deviation of the intensity values obtained in the filtered image. Finally, sub-pixel localisation of the molecules is performed by fitting a two-dimensional Gaussian function to the data using the weighted least squares method on a restricted domain around the molecule (7×7 pixels). As some molecules can be identified over consecutive frames, we perform a merging of the data acquired by the camera if the estimated distance between successive detections is less than 40 nm. Then, the position of the molecule is determined by using the average value of the positions estimated from the different frames. Using this strategy, we obtain approximately 24,000 different detections for the whole experiment. This number is limited by the weak activation power required to ensure that no more than one molecule is typically active at a single time on the area conjugated to the SPAD.

B. Decay rate estimations

In addition to EM-CCD images, we also record the arrival time of each photon detected by the SPAD. To deal with the large size of the resulting file (∼15 GB), the 10hour-long acquisition is split into several sequences of approximately 50 minutes. Then, we compute the number of detected photons as a function of time with a resolution of 500 µs. The intensity of background noise associated with this signal usually decreases during the experiment due to a decreasing number of activated molecules in the periphery of the detection area. Hence, the intensity time trace is Fourier filtered in order to remove low frequency components associated with temporal fluctuations longer than 30 s. Then, we consider that a molecule is potentially detected for each burst surpassing a given threshold that depends on the signal-to-noise ratio of the SPAD data. For each 50-minutes-long sequence, we set it to 2.6 times the standard deviation of the filtered signal. If another burst occurs within the typical blinking time scale (20 ms), it is attributed to the same molecule. In total, we identify approximately 14,000 events over the 10hour-long acquisition. This value is small in comparison to the number of detections obtained from camera data. Indeed, the area of the sample conjugated to the SPAD (see Fig. S2b) is smaller than the area over which the localisation is performed (∼ 1100 × 1100 nm). For each SPAD event, we build the associated decay histogram with a resolution of 16 ps in order to estimate the decay rate. To do so, the contribution of background noise is estimated by using close-by time intervals in which no burst can be identified. Then, the convolution of the instrument response function (IRF) and a decreasing mono-exponential function is fitted to the decay histogram using the least-squares method. The value of the decay rate is set to 10 ns -1 if the fit yields a value higher than this limit. Indeed, the IRF of the setup is characterised by a FWHM of approximately 240 ps (corresponding to 4 ns -1 ) and we consider that estimates above 10 ns -1 are not meaningful even after the deconvolution process. While sample heterogeneities could induce multi-exponential decays, the small number of photons detected by the SPAD from each molecule does not allow to resolve different lifetimes. For this reason, we restrict the analysis to a mono-exponential decay, which would therefore correspond to an average over different decays.

C. Temporal and spatial correlations

It is important to keep in mind that the SPAD does not include information about the position of the molecules. We therefore need to ensure that the lifetime information provided by the SPAD is properly associated with the position of the molecules provided by the EM-CCD camera. At the beginning of the experiment, the acquisition of both camera and SPAD data is started by using an in-house software, and we can expect a time offset of several milliseconds between the two different channels. In order to precisely determine this time offset, we build two binary representations respectively associated with the SPAD events and the camera detections (1 for a SPAD event or a camera detection, 0 otherwise). We then calculate the time correlation of these binary representations with a resolution of 500 µs, as shown in Fig. S2a for a typical sequence of 50 minutes. The maximum of this correlation coefficient gives an accurate estimate of the time offset between the camera and the SPAD. This delay is typically around 20 ms, which is consistent with the data acquisition procedure. Note that the correlation coefficient does not reach unity but is typically between 0.3 and 0.5. Indeed, the conditions required for the detection of a molecule by the camera and by the SPAD are different. In comparison to the SPAD, the camera is characterised by a larger field of view and a larger quantum efficiency. However, its lower temporal resolution makes the identification process less efficient for molecules characterised by fast temporal fluctuations. Hence we can expect some molecules to be detected by only one of the two detectors, resulting in a value smaller than unity for the maximum of the correlation coefficient. In order to characterise the spatial correlation between SPAD events and camera detections, we must identify the pixels of the camera that are conjugated to the area of the sample seen by the SPAD. Hence, we measure the response of the SPAD by scanning a fluorescent bead with a diameter of 100 nm over a large area in the sample plane. Figure S2b shows the number of photons detected by the SPAD as a function of the bead position. The FWHM value of the measured profile is of the order of 500 nm, as expected from the diameter of the confocal pinhole (50 µm) and the magnification of the optical system (×100). This response can be modelled by a function h(x, y) which is the convolution of a 500 nm gate and a two-dimensional Gaussian function.

D. Association conditions

Once the time offset between the camera and the SPAD is estimated and compensated, we can quantify the time overlap between a camera detection and a SPAD event.

To do so, we simply calculate the ratio of the time overlap ∆t ij to the time interval ∆t j corresponding to the SPAD event. The camera detection and the SPAD event are likely to be associated to the same molecule whenever this ratio is close to unity. We can then associate position and decay rate in the following situations:

• In 77% of the cases, the association between position and decay rate is straightforward. In such cases, only one camera detection is identified in the emission time ∆t j corresponding to a SPAD event.

In addition, this SPAD event is the only one identified in the emission time ∆t i corresponding to the camera detection. Therefore, the camera detection i and the SPAD event j can be associated.

• In 18% of the cases, several camera detections at different positions are identified in ∆t j . In such cases, we can estimate the number of photons to be detected by the SPAD from a given camera detection. Let x i and y i be the coordinates in the sample plane corresponding to a detection and N i the number of fluorescence photons measured by the camera, we can simply assume that the number of photons to be detected by the SPAD is proportional to N i h(x i , y i ). An association condition can thus be set on the base of the value taken by T ij = N i h(x i , y i )∆t ij /∆t j . After the identification of the detection k on the camera associated with the maximum value of T ij , we consider that the association between position and decay rate can be performed only if T kj > α a n i=1 T ij where n is the number of camera detections in ∆t j and α a is a threshold characterising the association condition. If α a is low, camera detections are more frequently associated to SPAD events. However, this increases the number of cases in which the sured decay histograms the sum of different cay histograms that cannot be properly separated by a post-processing analysis. As a trade-off, we use α a = 80% in the experiment.

• In 5% of the cases, several SPAD events are identified in ∆t i . Then, if the difference between these decay rates is smaller than 30%, we merge the SPAD events and we calculate the average decay rate. Otherwise, we evaluate the likelihood of each event to be the one corresponding to the camera detection, based on the number of fluorescence photons measured by the SPAD. To do so, we identify the event k associated with the highest number of photons N k and we perform the association between position and decay rate only if N k > α a n i=1 N i where N i is the number of photons associated with the overlapping SPAD events and α a is the threshold previously mentioned (α a = 80%).

a. Post-process filtering Two additional conditions are required in order to correctly perform the association between position and decay rate. For each molecule, at least 150 fluorescence photons must be detected on each detector. Moreover, the standard deviation of the Gaussian function fitted to the camera data must be smaller than 190 nm. These two conditions avoid the occurrence of false detections that would be due to noise. Using this procedure, we associate the position of 3,581 camera detections with their decay rate. We then perform post-processing filtering to account for the few remaining loopholes of the procedure. To do so, we compare each decay rate to the decay rate of the 10 closest detections. On average, this corresponds to a of 19 between the detection and its neighbours. Then, we perform an outlier identification based on the median absolute deviation (MAD). A decay rate Γ is rejected if the decay rates Γ k of the closest neighbours satisfy the following condition:

|Γ -Med(Γ k )| > α r Med |Γ k -Med(Γ k )| 0.675 , ( S1
)
where Med is the median operator and α r is a rejection threshold. The factor 0.675 is used so that MAD and standard deviation are approximately equal for large normal samples [8]. It should be noted that no outlier identification is performed if more than 50% of the neighbours have a decay rate equal to the upper limit previously mentioned (10 ns -1 ) since the right-hand side of Eq. (S1) equals zero in this case. With the approach expressed by Eq. (S1), using a small threshold α r allows the identification of many outliers but may also identify actual detections as outliers. As a trade-off, we use α r = 5 resulting in the identification of 6% of outliers. By removing them, the number of actual detections reduces to 3,352.

DENSITY AND INTENSITY MAPS

From data acquired by the EM-CCD camera, we can render a density map of the detected molecules (Fig. S3a), as for a usual single-molecule localisationbased super-resolution image reconstruction. In Fig. S3a, we observe strong density fluctuations due to an inhomogeneous labelling of our sample. However, note that for the purpose of obtaining a map of the LDOS, inhomogeneous labeling is not a limitation given a high enough spatial sampling, which underlines the robustness of our fluorescence lifetime measuring technique. It is important to underline that, in the image reconstruction in Fig. S3a, the strong density differences renders an image where black regions do not necessarily represent a lack of detections. In the case of biological applications, the labelling is specific to the protein of interest and thus density fluctuations represent structural changes of the sample which is not the case in our LDOS nanocartography. Additionally, we can also reconstruct a color map coding the measured fluorescence intensity for each detection (Fig. S3b). Note that, if several molecules are detected within the same area, we plot the average intensity. While the density of detected molecules is higher along the sides of the nanowire than on the substrate, we observe that the collected intensity is lower for the molecules on the nanowire. Indeed, although the excitation field is larger for the molecules on the sides of the nanowire, their radiative quantum yield is reduced due to coupling to non-radiative modes (surface plasmon modes and quenching).

NUMERICAL SIMULATIONS

Simulations are performed using the FDTD simulation software MEEP [9]. The relative permittivity of silver is modelled with a Lorentz-Drude model, the relative permittivity of the buffer solution is set to 1.77 and the relative permittivity of glass is set to 2.25. In order to estimate the influence of the excitation field on the observed density variations, we model the system in two dimensions, with a mesh resolution of 0.5 nm. The nanowire, located on a glass substrate, is illuminated by a plane wave at λ = 640 nm polarised perpendicularly to the nanowire, as in the experiment. In this configuration, a two-dimensional simulation gives the exact solution due to the invariance of the structure and the source along the longitudinal dimension. In contrast, in order to study the decay rate enhancement due to the nanowire, we model the system in three dimensions, with a mesh resolution of 1 nm. As the effect of the substrate on the decay rate is small due to the low contrast between the relative permittivities of the buffer solution and the glass coverslip, we perform the simulations without the substrate to limit the computational time. In each simulation, the emitter is modelled as an electric dipole source that generates a Gaussian pulse at λ = 670 nm, and the decay rate is estimated from the value of the electric field at the source position. We assume that the intrinsic quantum yield of Alexa Fluor 647 dyes is 0.33, as specified by the provider, in order to calculate the total decay rate enhancement.

DECAY HISTOGRAMS OF SINGLE MOLECULES

In this section, we show decay histograms for different molecules far from the nanowire and in its close vicinity, providing clear evidence of the decay rate enhancement. Figure S4 (a) to (f) shows the signal measured by the camera and by the SPAD during the experiment for three molecules characterised by different decay rate. Figure S4 (g) shows the associated decay histograms, together with mono-exponential fits. While the decay rate of the molecule far from the nanowire is not enhanced (molecule 1), the decay rate of the two molecules in the close vicinity of the nanowire show a strong decay rate enhancement (molecules 2 and 3). For the third molecule considered, the decay rate cannot be resolved by the current experimental setup as the decay histogram and the IRF are superimposed.

CRAMÉR-RAO ANALYSIS: POSITION ESTIMATIONS

To estimate the Cramér-Rao lower bound on the standard error of position estimators σ x,ȳ , we follow the approach described in [10]. The data acquired by the EM-CCD camera are modeled using the Airy function to describe the fluorescence signal, as well as a uniform background noise originating from the luminescence of the substrate. Then, we consider that the probability density function describing the number of photoelectrons per pixel is given by the convolution of the amplified signal and the Gaussian readout noise. The information matrix is calculated from this probability density function, and numerically inverted in order to compute the Cramér-Rao bound.

a. Point spread function We consider the simple situation in which a far-field microscope is used to collect the photons emitted by a single molecule located in the object plane. We assume that the 2-dimensional probability density function (PDF) describing the intensity distribution in the image plane can be expressed from the coordinates in the image plane noted (x , y ) and the coordinates of the molecule in the object plane noted (x 0 , y 0 ) as follows:

q(x , y ) = J 2 1 2πNA (x -M x 0 ) 2 + (y -M y 0 ) 2 M λ 0 π [(x -M x 0 ) 2 + (y -M y 0 ) 2 ] ,
(S2) where J 1 is the first-order Bessel function of the first kind, NA is the numerical aperture of the objective, M is the magnification and λ 0 is the free-space emission wavelength. The expectation of each data item -that is, the expectation of the value measured on each pixel by the camera -is then expressed as follows:

f i =N (x ,y )∈pixel q(x , y ) d x d y + N b (x ,y )∈pixel q b (x , y ) d x d y , (S3)
where N is the total number of photons emitted by the molecule and detected by the camera and N b is the number of photons due to background noise which follows a PDF noted q b (x , y ). In Eq. (S3), the integration is performed over the area that defines the considered pixel. Note that a dedicated study of the point spread function in our geometry could improve the prediction of position of the dipoles along the nanowire [11].

b. EM-CCD data model We can now derive a functional form for the likelihood function that describes the number of events measured on each pixel by the camera. Assuming that fluorescence photons detected by the camera are statistically independent, the number of photons impinging on each pixel during a given time interval follows a Poisson distribution of expectation f i . If we do not consider the additional noise arising from the detection process, the PDF associated with the observation of X photoelectrons on a given pixel is

p p i (X; θ) = f X i X! e -fi , (S4)
where θ are the parameters that must be estimated from the data (here, the parameters are the coordinates of the molecule). This sets the fundamental limit achievable by a perfect camera. However, the multiplication register of an EM-CCD camera enhances the number of generated photoelectrons in order to beat the readout noise of the camera, and the PDF followed by the number of photoelectrons generated by the process depends on the gain g. As shown in Ref. [12], this PDF noted p e i (X; θ) can be approximated, for large gain values, by

p e i (X; θ) =            e -f i , for X = 0, e (-X/g-f i ) fiX g I1   2 fiX g   X , for X > 0, (S5)
where I 1 is the first-order modified Bessel function of the first kind. In addition, the readout process induces a Gaussian noise on each pixel characterised by an expectation η g and a standard deviation σ g . This Gaussian noise can be described by the following PDF:

p g (X; θ) = 1 σ g √ 2π exp - (X -η g ) 2 2σ 2 g . ( S6
)
The PDF describing the readout noise of the camera is the same for all the pixels. Therefore, we can consider that the PDF describing the number of photoelectrons per pixel for a real EM-CCD camera is given by

p i (X; θ) = [p e i (X; θ)] * [p g (X; θ)] , (S7)
where the asterisk ( * ) represents the convolution product. Then, the information matrix can be numerically evaluated from its general expression given by [13] [I(θ)] jk = n i=1 E 1 [p i (X; θ)] 2 × ∂p i (X; θ) ∂θ j ∂p i (X; θ) ∂θ k .

(S8) c. Cramér-Rao bound After having experimentally measured the value of the parameters involved in the model, we can compute the Cramér-Rao bound on the variance of position estimators in order to evaluate a lower bound on the standard error σ x,y on the position estimates performed using one frame. Assuming that there is no preferred direction in space -this is not exactly true because of the shape of the pixels, but is a good approximation for squared pixels -the Cramér-Rao inequality reads

σ x,y ≥ 1 I xx = 1 I yy . ( S9
)
In the experiment, a molecule is typically detected on two successive frames. Its position is then estimated by the mean of the individual estimates, so that the standard error on the resulting position estimate is σ x,ȳ = σ x,y / √ 2. Different situations can then be compared: the fundamental is calculated using Eq. (S4) with N b = 0, the instrumental limit is calculated using Eq. (S7) with N b = 0, and the experimental limit is calculated using Eq. (S7) with the value of N b measured in the experiment.

CRAMÉR-RAO ANALYSIS: DECAY RATE ESTIMATIONS

To estimate the Cramér-Rao lower bound on the relative standard error of decay rate estimators σ Γ /Γ, we adopt a similar approach, described in [14]. In order to estimate the Cramér-Rao bound for our experiment, we model the fluorescence decay by the convolution of the IRF and an exponential distribution. After proper inclusion of time-dependent background noise in the model, it can then be considered that each point of the decay histogram follows a Poisson distribution. Thus, we can compute the information matrix from this distribution and numerically invert it in order to obtain the Cramér-Rao bound.

a. SPAD data model By modelling a molecule by a two-level system, the PDF that describes the photon emission time t is given by an exponential distribution. Since the agreement between experimental data and the mono-exponential model is satisfactory, we consider here that this model is relevant. Then, the PDF followed by the photon detection time measured by the experimental system is q(t) = q irf (t) * Γe -Γt , (S10)

where q irf (t) is the PDF describing the IRF of the setup.

From this expression, we can find the expectation of each data item; that is, the expectation of each data point of the decay histogram. We obtain

f i = N +∞ l=0 ti+1+lT ti+lT q(t) d t + N b ti+1 ti q b (t) d t , ( S11
)
where N is the number of photons emitted by the molecule and detected by the system, N b is the number of detected photons due to background noise which follows a PDF noted q b (t), and T is the repetition period of the laser. If the fluorescence lifetime of the molecule is much smaller than the repetition period, only the first term of the sum in Eq. (S11) is significant.

In general, SPADs have negligible readout noise and the dark count rate contributes to the background noise. Thus, we can model the distribution of photons detected for each data point by a Poisson distribution of expectation f i . The PDF associated with the observation of X events on a given data point is then expressed by

p i (X; θ) = f X i X! e -fi . (S12)
The set of parameters that must be estimated from the data is θ = (N, Γ), while we estimate q irf (t), q b (t) and N b with independent measurements. Then, the information matrix can be calculated from Eq. (S8). b. Cramér-Rao bound After having experimentally measured the value of the parameters involved in the model, we can compute the Cramér-Rao bound on the standard error σ Γ on the decay rate estimates. The Cramér-Rao inequality can be expressed as

σ Γ Γ ≥ 1 √ N × F (T, N b , q irf , q b , n) , ( S13
)
where n is the number of data points and F is calculated by numerically inverting the information matrix [14].

  FIG. 1: Super-resolved LDOS mapping of a silver nanowire. (a) The EM-CCD camera acquires 31 frames per second, with an exposure time of 30 ms per frame and a field of view of tens of microns. A single fluorescent molecule is detected on the third frame of the sequence shown here. An image cropped around the molecule is shown here for the sake of simplicity. (b) The position of the molecule is estimated by fitting a two-dimensional Gaussian function to the measured point spread function (PSF). The PSF covers an area of ∼ 9 pixels (pixel size = 160 nm). (c) At the same time, the SPAD detects a fluorescence burst from this molecule. (d) Short laser pulses were used to excite this molecule. For each photon detected during the fluorescence burst, the time difference between excitation and emission can be determined with picosecond precision. (e) These photons are used to construct a decay histogram. The convolution of the IRF and a decreasing mono-exponential function is then fitted to this histogram in order to estimate the fluorescence decay rate of the molecule. The IRF is shown in orange. (f) Reconstructed decay rate map. Each dot represents the position of a detected molecule, and its diameter is fixed at 15 nm, which is the typical full width at half maximum (FWHM) of the probability density function followed by the position estimates. If several molecules are detected within the same area, we show their average decay rate on the map.

FIG. 2 :

 2 FIG. 2: (a) Number of detected molecules as a function of the distance to the nanowire axis. (b) Time-averaged intensity of the excitation field in the vicinity of the silver nanowire calculated from the results of a FDTD simulation. On these figures, dashed lines represent the estimated position of the nanowire edges. (c) Distribution of decay rate versus distance to the wire axis. The highest decay rate that can be measured, limited by the IRF of the setup, is 10 ns -1 . (d) Decay rate enhancement as a function of the distance to the nanowire for the three orientations of the dipole moment. The inset shows a cross-section of the system numerically studied.

FIG. 3 :

 3 FIG. 3: (a) Top: Cramér-Rao bound on the standard error on the position estimates as a function of the number of fluorescence photons detected by the camera. Bottom: Distribution of the number of fluorescence photons detected by the camera from the single molecules. (b) Top: Cramér-Rao bound on the standard error on the decay rate estimates as a function of the number of fluorescence photons detected by the SPAD. Bottom: Distribution of the number of fluorescence photons detected by the SPAD from the single molecules. On these figures, dashed lines represent the threshold condition N > 150 used for data analysis.

  FIG. S1: Optical setup. The excitation laser (λ = 640 nm), together with the photo-activation laser (λ = 405 nm) and the laser used for sample stabilisation (λ = 568 nm), illuminate the sample via a high numerical aperture oil objective (NA=1.4). A lens (f=300 mm) is located on the excitation path to ensure wide-field illumination. Fluorescence from the sample is filtered with a dichroic mirror (DM) and passes through a tube lens (TL). A 50:50 beamsplitter (BS) splits the light towards an EM-CCD camera and a SPAD. The sample under study contains photo-activated single molecules in the near-field of a silver nanowire (NW).

  FIG. S2: (a) Correlation coefficient calculated from binary representations of the SPAD events and the camera detections. A dashed line represents the estimated time offset between the two channels. (b) Measured response of the SPAD while scanning a fluorescent bead in the sample plane.

  FIG. S3: (a) Density and (b) intensity maps reconstructed from the 14,546 molecules detected by the EM-CCD camera.

  FIG. S4: Camera images and signal measured by the SPAD for a molecule far from the nanowire [sub-figures (a) and (b)], for a molecule close to the nanowire with Γ/Γ0 ∼7 [sub-figures (c) and (d)], and for a molecule close to the nanowire with Γ/Γ0 >15 [sub-figures (e) and (f)]. The associated decay histograms are shown in sub-figure (g).

  D. Bouchet, J. Scholler, G. Blanquer, Y. De Wilde, I. Izeddin, and V.Krachmalnicoff

	Institut Langevin, ESPCI Paris, CNRS, PSL University, 1 rue Jussieu, 75005 Paris, France
	1. EXPERIMENTAL SETUP
	A. Sample preparation

Acknowledgments

The authors thank S. Bidault for helping in sample preparation, A. C. Boccara for sharing his insights about the manuscript and I. Rech, A. Gulinatti and A. Giudice for providing the PMD-R detector.

Funding information

This work was supported by LABEX WIFI (Laboratory of Excellence ANR-10-LABX-24) within the French Program Investments for the Future under reference ANR-10-IDEX-0001-02 PSL*, by the Programme Emergences 2015 of the City of Paris, and by ANR-17-CE09-0006 SimpleLife.