Daniele De Bernardis

Mathieu Jeannin

Jean-Michel Manceau

Raffaele Colombelli

Alessandro Tredicucci

Iacopo Carusotto

Magnetic-field-induced cavity protection for intersubband polaritons

We analyse the effect of a strong perpendicular magnetic field on an intersubband transition in a disordered doped quantum well strongly coupled to an optical cavity. The magnetic field changes the lineshape of the intersubband optical transition due to the interface roughness of the quantum well from a Lorentzian to a Gaussian one. In this regime, a novel form of magnetic-field-induced cavity protection sets in, which strongly reduces the polariton linewidth to the cavity contribution only. Implications of our results for fundamental studies of nonlinear polariton dynamics and for technological applications to polariton lasers are finally highlighted.

Strongly confined quasi two-dimensional (2D) electronic systems have a long research history [START_REF] Ando | Electronic properties of two-dimensional systems[END_REF], and still remain under the researchers' spotlight due to their large importance for both fundamental science and technological developments. In recent years, 2D electronic systems supporting intersubband (ISB) transitions polarized perpendicular to the quantum well plane have been combined with good-quality electromagnetic cavities giving rise to strong light-matter coupling effects and opening the field of intersubband polaritons [START_REF] Dini | Microcavity Polariton Splitting of Intersubband Transitions[END_REF].

In particular, by adjusting the electronic density it is possible to vary the dipole moment of the ISB on a quite broad range, bringing these polaritonic devices from the usual strong-coupling regime to the more exotic ultrastrong coupling (USC) regime. This regime was experimentally achieved with ISB and other various solid-state platforms [START_REF] Dietze | Ultrastrong coupling of intersubband plasmons and terahertz metamaterials[END_REF][START_REF] Jouy | Transition from strong to ultrastrong coupling regime in mid-infrared metaldielectric-metal cavities[END_REF][START_REF] Askenazi | Ultrastrong light-matter coupling for designer Reststrahlen band[END_REF][START_REF] Todorov | Ultrastrong Light-Matter Coupling Regime with Polariton Dots[END_REF] and is of great interest for its strongly modified ground-state properties [START_REF] Ciuti | Quantum vacuum properties of the intersubband cavity polariton field[END_REF][START_REF] Frisk Kockum | Ultrastrong coupling between light and matter[END_REF][START_REF] Forn-Díaz | Ultrastrong coupling regimes of light-matter interaction[END_REF]. More in general, the large flexibility of ISB polaritons holds great promise for the study of a wide range of cavity quantum electrodynamics (cQED) effects from both a fundamental [START_REF] Ciuti | Quantum vacuum properties of the intersubband cavity polariton field[END_REF][START_REF] Knorr | Intersubband polariton-polariton scattering in a dispersive microcavity[END_REF][START_REF] Nguyen-Thê | Effective polariton-polariton interactions of cavity-embedded two-dimensional electron gases[END_REF][START_REF] Halbhuber | Non-adiabatic stripping of a cavity field from electrons in the deep-strong coupling regime[END_REF] and an applied perspective [START_REF] Colombelli | Quantum cascade intersubband polariton light emitters[END_REF][START_REF] Sapienza | Electrically Injected Cavity Polaritons[END_REF][START_REF] Liberato | Stimulated Scattering and Lasing of Intersubband Cavity Polaritons[END_REF]. Despite the great achievements reached using ISB polaritonic devices in both the mid-infrared and THz regime, a serious obstacle hindering a full development of this research program is caused by the relatively large linewidth of the ISB transition [START_REF] Colombelli | Perspectives for Intersubband Polariton Lasers[END_REF], mostly due to interface roughness of the semiconductor quantum well nanostructures [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF].

In this article we theoretically investigate a method to dramatically reduce the ISB polariton linewidth by suppressing the contribution of interface roughness and bringing it down to the cavity contribution only. The idea is to combine a strong light-matter coupling to a single electromagnetic cavity mode with a strong static magnetic field, so as to enter a regime of cavity protection [START_REF] Houdré | Vacuum-field Rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system[END_REF]. The strong perpendicular magnetic field quenches the in-plane kinetic energy of the electrons and strongly localises them in the disorder potential due to the interface roughness of the quantum well. As a consequence, the line shape of the ISB transition changes from a standard Lorentzian shape to an almost Gaussian one. Due to the fast decay of the Gaussian tails, the strong coupling to the cavity mode then effectively suppresses the inhomogenous broadening of the ISB transition, leaving a polariton linewidth which is then only limited by the cavity component. Differently from previous works on cavity protection in multi-quantum-well devices [START_REF] Manceau | Immunity of intersubband polaritons to inhomogeneous broadening[END_REF] and in various ensembles of emitters in cavity-QED devices [START_REF] Kurucz | Spectroscopic properties of inhomogeneously broadened spin ensembles in a cavity[END_REF][START_REF] Diniz | Strongly coupling a cavity to inhomogeneous ensembles of emitters: Potential for long-lived solid-state quantum memories[END_REF][START_REF] Putz | Protecting a spin ensemble against decoherence in the strong-coupling regime of cavity QED[END_REF][START_REF] Zhong | Interfacing broadband photonic qubits to on-chip cavityprotected rare-earth ensembles[END_REF][START_REF] Breeze | Room-temperature cavity quantum electrodynamics with strongly coupled Dicke states[END_REF][START_REF] Baghdad | Spectral Engineering of Cavity-Protected Polaritons in an Atomic Ensemble with Controlled Disorder[END_REF], the magnetic-field-induced cavity protection mechanism predicted in our work is based on the in-plane dynamics of electrons in a single quantum well, which completely changes nature under the effect of the strong magnetic field. Our calculations anticipate a dramatic enhancement of the polariton quality factor, which opens new perspectives to experimental research in this field and to technological applications.

The article is organised as follows. In section I we introduce the model used to describe ISB polaritons in the presence of interface roughness in the quantum well and calculate the cavity transmittivity in the strong coupling regime. In section II we characterise the impact of disorder on the ISB linewidth in the absence of a surrounding cavity and we extend the known results to the case of a strong perpendicular magnetic field. In section III we translate the concept of cavity protection to our context of ISB polaritons in strong magnetic field and we show its implications: the main regimes are characterized and experimental implementations are discussed. In section IV we summarize the main results and conclusions of this work. Additional information on the derivations are given in the Appendices.

I. MODEL

In this section we introduce the model to describe a planar quantum well (QW) coupled to the electromagnetic field of an optical cavity, as represented in Fig. 1 a) Sketch of the system. A semiconductor twodimensional quantum well is inserted in a metallic cavity. The electric field of the cavity mode oscillates along the zdirection, perpendicular to the quantum well plane. The electronic intersubband transitions in the quantum well give rise to an electronic polarization directed along the same direction of the cavity electric field, realising a strong light-matter coupling regime between the cavity and the quantum well. b-c) Schematic view of the electronic subbands in the disordered quantum well in absence/presence of the external magnetic field perpendicular to the quantum well plane. The lowest optically-active excitations are transitions between the lowest two subbands due to the electrons below the Fermi energy EF . Differently from the clean case, these transitions are not only fully vertical and each electron can jump to different eigenstates of the upper subband.

by the fabrication process. The inclusion of this disorder forces us to go beyond the standard theory of intersubband transitions based on the bosonization of the collective ISB electronic excitations [START_REF] Todorov | Intersubband polaritons in the electrical dipole gauge[END_REF] and include the underlying fermionic degrees of freedom. However in the small disorder, small excitation limit, it is still possible to derive a sufficiently simple and manageable theory of light-matter interactions with ISB transitions. This theory captures well the effect of disorder and can be immediately extended with the inclusion of a static external magnetic field orthogonal to the QW.

A. Electronic states in a disordered quantum well in the presence of a static magnetic field

We consider here an electronic system strongly confined in the z direction and free in the (x, y)-plane, in such a way to realise a two dimensional quantum well (QW) with well defined energy subbands [START_REF] Liu | Intersubband Transitions in Quantum Wells: Physics and Device Applications II , Semiconductors and semimetals[END_REF]. An homogeneous magnetic field with amplitude B is pointing in the z-direction, perpendicular to the two dimensional QW. Assuming the symmetric gauge for the magnetic vector potential, A( r) = B/2(-y, x, 0), the single parti-cle Hamiltonian for each electron reads

H e = p -e A 2 2m * + p 2 z 2m * + U (z, r), (1)
where p = (p x , p y ), r = (x, y) are respectively the electron's in-plane momentum and position, m * is the electron effective mass, e its charge and U (z, r) is the single particle confining potential. Under the assumption of a strong confinement on the z-direction we can employ a Born-Oppenheimer-like approximation between the different directions, and write the single particle eigenfunctions as

ψ(z, r) ≈ ζ n (z, r)ϕ n k ( r), (2)
where the subband wavefunctions ζ n (z, r) are assumed to depend only parametrically on the in-plane position r and for each r solve the eigenproblem

p 2 z 2m * + U (z, r) ζ n (z, r) = ε n ( r)ζ n (z, r) (3)
along the z direction. Here, n = 1, 2 . . . are integer numbers and ε n ( r) indicate the Born-Oppenheimer potential energy interfaces.

For each subband n, the in-plane eigenfunctions ϕ n k ( r) are thus obtained by diagonalising the in-plane Hamiltonian

H n = p -e A 2 2m * + ε n ( r) (4)
and are labeled generically by the k quantum number (or array of numbers). In the case of a clean QW with no magnetic field, this index reduces to the in-plane momentum.

We assume now that the z-confinement is characterised by a single length scale L z = L z ( r) which depends on the in-plane position r, and that the potential energy interface depends from the in-plane position only through this characteristic length ε n ( r) = ε n (L z ( r)). Moreover we assume that L z ( r) is given by a fixed length over the whole plane plus a small fluctuating part L z ( r) = L qw + δL( r) which describes the roughness of the QW interface. Expanding the potential energy interface at the first order in δL we arrive to the final expression of the in-plane electronic Hamiltonian of each subband as H n ≈ H n + ε qw n , where the in-plane electronic Hamiltonian is given by

H n = p -e A 2 2m * + δU n ( r) , (5)
ε qw n = ε n (L qw )
is the energy of the n subband of the clean QW, and δU n ( r) = ∂ L ε n (L z )•δL( r) is the positiondependent energy shift due to the disorder.

The effect of the external magnetic field B on the energy distribution of the in-plane eigenstates is sketched in Fig. 1(b,c): when the cyclotron frequency ω B = eB/m * is larger than the typical energy width of each potential energy interface δU n , the electronic dispersion of each subband changes from the usual parabola to discrete Landau levels.

In our description we completely neglect the effect of electron-electron interactions. These are known to have a small effect in the ISB optical linewidth [START_REF] Luin | Interplay between disorder and intersubband collective excitations in the two-dimensional electron gas[END_REF], which is mainly limited by the interface roughness [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF][START_REF] Nikonov | Collective Intersubband Excitations in Quantum Wells: Coulomb Interaction versus Subband Dispersion[END_REF][START_REF] Ullrich | Theory of the Linewidth of Intersubband Plasmons in Quantum Wells[END_REF]. For this reason, we chose to avoid further complications and focus on the single particle dynamics only, so to highlight in the most transparent way the interplay between roughness disorder and magnetic field. Including the effect of Coulomb interactions in our description of the in-plane electron dynamics goes beyond the scope of this paper and is left for a future work.

In spite of this, it must be kept in mind that the effect of electron-electron interactions is still present in the presence of magnetic field [START_REF] Kempa | Intersubband Transport in Quantum Wells in Strong Magnetic Fields Mediated by Single-and Two-Electron Scattering[END_REF][START_REF] Kempa | Electron-electron scattering in strong magnetic fields in quantum well systems[END_REF] and, in particular, the effect of Coulomb interactions onto the electron motion along the z-axis is completely unaffected by the magnetic field along z. As such, it keeps giving an important contribution to the nonlinearity of the ISB response to strong electromagnetic fields [START_REF] Cominotti | Theory of coherent optical nonlinearities of intersubband transitions in semiconductor quantum wells[END_REF]. This implies that our results are directly applicable to the on-going quest of developing ISB polaritonic devices with strong optical non-linearities and polariton lasers, where a reduced polariton linewidth can be a game-changing improvement [START_REF] Knorr | Intersubband polariton-polariton scattering in a dispersive microcavity[END_REF].

B. Intersubband transitions and cavity-polaritons

After having characterized the single-particle eigenstates of electrons, we consider here that the two dimensional QW is placed inside a cavity enclosed within a pair of metallic plates parallel to the QW (x, y)-plane, as represented in Fig. 1. We assume that the electronic ISB transitions mostly couple to the single TM 0 cavity mode, which, in the long-wavelength approximation, can be taken to be homogeneous along the plane, with frequency ω c and zero-field amplitude E 0 = ω c /(2 0 V ) directed in the z-direction and determined by the cavity volume V .

The light-matter Hamiltonian is then given by the dipole-gauge Hamiltonian [START_REF] Todorov | Intersubband polaritons in the electrical dipole gauge[END_REF],

H int ≈ dzd 2 r P ( r, z) • E cav ( r, z) , ( 6
)
where P is the electronic polarization of the QW and E cav ≈ u z E 0 (a + a † ) is the electric field of the cavity. Assuming we are below the ultra-strong light-matter coupling regime, we can for simplicity discard the ∼ P 2term [START_REF] Bernardis | Cavity quantum electrodynamics in the nonperturbative regime[END_REF][START_REF] Rokaj | Light-matter interaction in the long-wavelength limit: no ground-state without dipole self-energy[END_REF].

Introducing the creation/annihilation spinless electronic operators Ψ( r, z), satisfying the Fermionic anticommutation relation {Ψ( r, z) † , Ψ( r , z )} = δ(zz )δ (2) ( r -r ) we can rewrite the electronic polarization as P ( r, z) ≈ ezΨ † ( r, z)Ψ( r, z) u z [START_REF] Todorov | Intersubband polaritons in the electrical dipole gauge[END_REF], where e is the electron charge. The second quantized light-matter coupled Hamiltonian reads

H tot = ω c a † a + dz d 2 r Ψ † ( r, z)H e Ψ( r, z)+ + e dz d 2 r zΨ † ( r, z)Ψ( r, z) • E 0 (a + a † ). (7)
Restricting ourselves to QWs with a fixed number of electrons N e , we define the plasma frequency of the electronic transitions as

ω P = e 2 0 m * N e V . (8)
The lowest QW intersubband frequency transition and its z-direction oscillator strength are defined as

ω qw = ε qw 2 -ε qw 1 f qw = 2m * ω qw z 2 21 , (9)
in terms of the dipole matrix element z 2

n n = | ζ n |z|ζ n | 2 between
the n, n single particle eigenfunctions defined in Eq. (2). To avoid complications stemming from the implicit parametric dependence of the matrix elements z n n ∼ L qw +δL( r) on the in-plane position r, we consider an average value z 2 n n over the (x, y)-plane, neglecting its dependence on the fluctuating QW thickness. This assumption is well justified within the weak disorder assumption and corrections are of higher order in δL( r).

We then re-express the electron field operator in terms of the single particle eigenfunctions Ψ( r, z) = n k ζ n (z, r)ϕ n k ( r)c n k , where the creation operators c n k labelled by the subband n and in-plane k indices satisfy Fermionic commutation rules {c n k , c † n k } = δ n,n δ k,k . Within a rotating wave approximation, we neglect the counter-rotating terms in the light-matter interaction and we restrict our attention to the lowest intersubband transition between n = 1 and n = 2. In this way, we obtain the cavity-plasma Hamiltonian

H tot ≈ ω c a † a + n=1,2 k ε qw n + ω n (k) c † n k c n k + Ω R 2   a • k,k Λ k k c † 2 k c 1 k + h.c.   , (10)
where the strength of the light-matter coupling is quantified by the Rabi frequency

Ω R = ω P f qw ω c ω qw . (11)
Here, ω n (k) are the eigenfrequencies of H n as defined in Eq. ( 5) and the k, k transition matrix element is given by the in-plane wavefunction overlap between states in the two different subbands

Λ k k = ϕ 2 k |ϕ 1 k √ N e . ( 12
)
Note that the use of the dipole gauge allows us to truncate the Hilbert space to the two lowest subbands without introducing potentially dangerous spurious effects [START_REF] Bernardis | Breakdown of gauge invariance in ultrastrongcoupling cavity QED[END_REF].

C. Cavity transmission and quantum well optical density

Assuming that the cavity is resonant with the lowest QW intersubband transition ω c ∼ ω qw , we can use the cavity-plasma Hamiltonian in Eq. ( 10) to derive the input-output equations for the system [START_REF] Gardiner | Quantum Noise[END_REF] and then the cavity optical transmission of a weak probe of frequency ω at normal incidence (see App. A for further details on the derivation ),

T c (ω) = - γ c /2 ω -ω c + i γc 2 + Ω 2 R 4 χ qw (ω) , ( 13
)
where γ c is the cavity loss rate and microscopic details of the QW are summarized by the optical response

χ qw (ω) = - k ≤k F ,k |Λ k,k | 2 ω -ω qw -ω k k + iκ . ( 14
)
Here

ω k k + ω qw = ω 2 (k) -ω 1 (k ) + ω qw is the frequency of the (1, k ) → (2, k
) transition, κ is a small phenomenological energy loss rate introduced to regularize the response, and the sum over initial states k is restricted to the occupied states below the Fermi level k F of the lowest subband.

Using the Sokhatsky identity lim κ→0 Im [1/(ω + iκ)] = -πδ(ω) we can introduce the so-called QW optical density

ρ qw (ω) = lim κ→0 1 π Im [χ qw (ω)] = k ≤k F ,k |Λ k,k | 2 δ(ω -ω qw -ω k k ). ( 15
)
This quantity is normalised to one dω ρ qw (ω) = 1, as can be verified by considering that 1 = k |ϕ n k ϕ n k | and using again the Sokhatsky identity, and is equivalent to the density of states typically used in the literature [START_REF] Ando | Lineshape of inter-subband optical transitions in space charge layers[END_REF], which is derived from the electronic Green's function, as it is shown in App. B. The QW optical response can be rewritten only in terms of the optical spectral density as

χ qw (ω) = -dω ρ qw (ω ) ω -ω + iκ , (16)
making this quantity the central object of our investigation.

II. LINEWIDTH BROADENING DUE TO ROUGHNESS DISORDER

To understand the behaviour of the cavity trasmission T c (ω) in the different regime, we first need a clear un- derstanding of the QW optical response χ qw (ω), which represents the optical susceptibility of the QW to an external probe with incident frequency ω. This reduces to study how the optical density ρ qw (ω) is modified by the disorder and the external magnetic field.

In a clean sample, or when the disorder is exactly the same in the two subbands the transition matrix is Λ k k = δ k k and all transition frequencies coincide, so to have ω k k = 0. From Eq. ( 15), it then follows that the QW optical density is a delta function centered at the QW frequency ρ qw (ω) ∼ δ(ω -ω qw ).

In this section we will examine how this picture is destroyed by the disorder and, successively, what is the interplay between disorder and magnetic field. To be concrete we keep the discussion in the simple, paradigmatic, example of the QW with infinite well confinement. However, in order to make our results general and independent from the precise shape of the confining potential in the z-direction, we will introduce a set of natural units to express all the quantities in adimensional form. In this way all the results of this section applies as well to any other type of z-confinement after providing the appropriate rescaling of variables.

A. Simple model for the interface roughness disorder

When the system is instead affected by the interface roughness giving a different disorder potential δU n for each subband, several k → k transitions contribute to the QW optical density ρ qw (ω) and, in the large-system limit, this becomes a smooth continuous distribution, broadened around ω qw , with a Lorentzian linewidth Γ qw . The quality factor of the QW ISB transition is then given by Q qw = ω qw /Γ qw .

Specifically, we focus here on the case of an infinite well (or square box) confinement in the z-direction. In this way we have that

ε n (L) = 2 (πn) 2 2mL 2 . ( 17
)
We also introduce a reference length scale that depends only from the QW effective electron mass and central frequency

l qw = m * ω qw . (18)
In the case of an infinite well with length L qw , where the central frequency is given by the lowest transition l qw = 2/3L qw /π ≈ L qw /4. From now on we will give all the lengths in units of this length scale.

For simplicity, we specialise on the case of a Gaussiandistributed interface roughness, for which δL( r) = 0, and δL( r)δL( r

) = ξ 2 0 /(2π) exp -| r -r | 2 /ξ 2 c .
Here the overline bar indicates the disorder average and ξ 0 is the roughness amplitude, while ξ c is its correlation length. We can then define a dimensionless parameter which controls the disorder amplitude also in the general case of an arbitrary confinement and an arbitrary disorder roughness

η dis = |∂ L ε 1 (L qw )|ξ 0 ω qw = 2 3 3/2 ξ 0 πl qw ≈ 1 6 ξ 0 l qw , (19)
where the last equality holds for the specific case of the infinite square box potential (see App. C for major details). The disorder potential for each subband n reads

δU n ( r)/ω qw = n 2 η dis ∆( r), (20)
where ∆( r) = δL( r)/ξ 0 . In Fig. 2(a) an example of the disorder potential is reported together with its Fourier transform in Fig. 2(b). Here the parameter's choice is inspired by the typical values found in GaAs/AlGaAs systems, see for instance [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF].

It is worth to highlight that the use of energy, length units and disorder amplitude in terms of ω qw , l qw and η dis (as detailed in the Appendix C) is important to keep our description as general as possible and to make it independent from the specific shape of the confining QW potential in the direction of the ISB transition (zdirection). For instance, the use of these units makes the application to harmonic z-confinement case straightforward. The only difference is in how the results are re-parametrised in the disorder ξ 0 and correlation ξ c lengths. Specifically, we will have that η dis = 2ξ 0 /l qw and δU n ( r)/ω qw = n η dis ∆( r). Notice the more favourable scaling of the harmonic confinement ∼ n.

B. Intersubband linewidth at B = 0

In the absence of magnetic field, B = 0, the QW optical density is approximatively given by a Lorentzian distribution [START_REF] Ando | Lineshape of inter-subband optical transitions in space charge layers[END_REF][START_REF] Ndebeka-Bandou | Relevance of intraand inter-subband scattering on the absorption in heterostructures[END_REF], Fermi energy EF / ωqw ≈ 0.5. In both panels, the red dashed line is the Unuma-linewidth given by Eq. [START_REF] Zhong | Interfacing broadband photonic qubits to on-chip cavityprotected rare-earth ensembles[END_REF]. To numerically implement Eq.( 15) we used a finitewidth delta function, with small linewidth γ δ /ωqw ≈ 0.002, which is then subtracted from the numerical data.

ρ qw (ω) ≈ 1 2π Γ qw (ω -ω qw ) 2 + Γ 2 qw /4 , (21)
where Γ qw is its full width at half maximum (FWHM). This behaviour is confirmed looking at Fig. 3(a), but we also notice some additional features. The long tails of the distribution decay somehow faster than a proper Lorentzian; moreover the optical density ρ qw (ω) is not completely symmetric and its asymmetry can be increased or reduced by changing the Fermi energy E F , see Fig. 3(b). In particular reducing the Fermi energy reduces the extension of the low energy tail (on the left side of the peak). All these effects can be understood considering that the left tail is due to the electrons that jump to an energy level of the second subband which has lower energy than their initial level, so to have ω k k < 0. Assuming the Fermi energy larger than the disorder energy scale, we have that the lowest value of the in-plane electronic transition is approximately given by the Fermi energy ω qw + ω k k min ≈ ω qw -E F / . On contrary the value of the Fermi energy has no influence on the tail on the right side of the peak, since the electrons have no limitation in jumping toward higher energies.

As we can see from Fig. 4, the FWHM Γ qw scales quadratically both in the disorder strength and in the correlation length [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF][START_REF] Ando | Lineshape of inter-subband optical transitions in space charge layers[END_REF][START_REF] Grange | Atomic-Scale Insights into Semiconductor Heterostructures: From Experimental Three-Dimensional Analysis of the Interface to a Generalized Theory of Interfacial Roughness Scattering[END_REF]

Γ qw /ω qw ∼ η 2 dis , ξ 2 c /l 2 qw , (22)
as long as the Fermi length is much longer than the disorder correlation length ξ c k F 1. Following the semi-analytical approach developed in [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF] by Unuma et al., it is possible to derive an analytic expression for the ISB linewidth which holds in the zero energy limit of the electronic scattering against the disorder

Γ qw, U ω qw ≈ 9 2 ξ 2 c l 2 qw η 2 dis . (23)
As we can see from Fig. 4 this formula, that we call Unuma-linewidth, fits very well the ISB linewidth extracted from the numerical simulations in the regime of small ξ c k F 1.

C. Intersubband linewidth at B = 0

When the magnetic field is turned on, B = 0, the situation changes drastically. In particular when the cyclotron frequency ω B = eB/m * exceeds the energy scale of disorder ω dis = ω qw η dis ξ c /l qw we can no longer think about the electrons in terms of free particles diffusing in a disordered landscape, but instead we need to switch to a description in terms of Landau levels. Having this in mind we can define a reference magnetic field that approximatively sets the border of this transition

B 0 = η dis m * ω qw e ξ c l qw = η dis ξ c l qw ω qw 2µ * B ( 24
)
Here, µ * B = e /(2m * ) is the effective Bohr magneton relative to the effective mass m * of the electron. In the rest of the article we will give all magnetic fields expressed in this rescaled unit. In the last part we will give a more precise analysis on the quantitative conditions to achieve the regime of cavity protection.

For a completely filled Landau level, the Lorentzian optical density is broken into a series of equispaced Gaussian peaks, as we can see from Fig. 5(a). The central peak represents all transitions from a given Landau level 1 below the Fermi energy in the 1-th subband, to the same respective Landau level in the 2nd subband, 1 = 2 . The side peaks represent instead transitions to other Landau levels of the upper subband, 2 = 1 . For partially filled bands, the peaks corresponding to transitions starting from the uppermost populated Landau level become asymmetric, yet with no qualitative consequence on the overall conclusions of our study. In Fig. 5(b) we plot the optical density with fixed magnetic field when the Fermi energy falls in the middle of the zeroth/first/second broadened Landau level, E F = 0, ω B , 2ω B . We see that the asymmetry is more prominent when the Fermi energy is in the lowest Landau level, and then is progressively washed out when E F is in the higher Landau levels.

The spectral distance between neighbouring Gaussian is approximately given by the cyclotron frequency, ∼ ω B , as a consequence of the transitions between neighbouring Landau levels. Focusing on the properties of the central Gaussian, centered on ω qw , we call its linewidth Γ B . Contrary to the non magnetic case, discussed in the previous section, here we expect that the linewidth of each of these Gaussian lobes scales linearly with the disorder strength Γ B /ω qw ∼ η dis . This can be understood from the fact that the disorder plays the role of a small perturbation on a degenerate system, bringing corrections at the linear order instead of the usual second order.

As in the non-magnetic case, the disorder correlation length ξ c plays an important role in determining the width of each Gaussian, but in contrast to the nonmagnetic case, here the system has an intrinsic length scale, given by the magnetic length

l B = eB . ( 25
)
We then expect that the FWHM of each Gaussian depends non-trivially on ξ c /l B . Our aim is now to have more quantitative insights on the dependence between the width of the Gaussian lobes and the system parameters. In order to do so we start ). The red/green dashed lines are given by the Hikami linewidth ΓB, H/ωqw while the black solid line is the Unuma linewidth Γqw, U/ωqw. Disorder parameters: ξc/lqw = 0.5. b) Gaussian linewidth ΓB (purple dots) extracted from numerical simulations of the optical density as a function of the magnetic field B. The dashed red line is given by the Hikami linewidth ΓB, H/ωqw, while the dashed-dotted line is given by the small disorder expansion of the Hikami linewidth Γ

B, H . Disorder parameters as in Fig. 5: η dis = 2/15 ≈ 0.13 and ξc/lqw = 0.5. c) Same plot as in b) but with different parameters in order to highlight the saturation of the linewidth at high magnetic field. The red dashed line is given by the Hikami-linewidth ΓB, H/ωqw. The purple dotted line is the saturation value of the Hikami linewidth given by Γ ∞ B, H /ωqw ≈ 0.094. Disorder parameters: η dis = 1/30 ≈ 0.03 and ξc/lqw = 4. In all plots the Fermi energy is EF = 1.5 ωB, in such a way that the central Gaussian is determined by only the two lowest Landau levels. To numerically implement Eq.( 15) we used a finite-width delta function, with small linewidth γ δ /ωqw ≈ 0.003, which is subtracted from the numerical data.

from the observation that not all the transitions are important in this regime of strong magnetic field. Indeed, considering the central Gaussian lobe of the optical density, we realise that only the intra-Landau-level transitions with the same k n are actually relevant. This is justified by the fact that the disorder landscape is the same in the two subbands, and differs only in its amplitude. Specifically, in the 2-subband it is 4 times larger than in the 1-subband. The system is then strongly localised but each state in the ground subband can overlap only with the state in the upper subband localised in the same region. So, in Eq. ( 15), we can approximate Λ ( 1k1 ),( 2,k2) ≈ δ 1 2 δ k1k2 . The energy difference of each transition from a state ( 1 , k 1 ) in the 1-subband to the corresponding state ( 2 , k 2 ) in the 2-subband is then given by

ω qw + ω ( 2k2)( 1 k1) ≈ ω qw + 3(ω 1k1 -ω B 1 )
, where ω 1k1 is the k 1 -th eigenenergy of the 1 -th Landau level of the 1-subband.

The considerations above are particularly important when we restrict each band to only include the lowest Landau level (LLL). This is well justified when the Fermi energy is smaller than the cyclotron frequency, E F < ω B , and it is a crucial assumption in order to carry on the calculation analytically. However, we will see that the result extracted in the LLL remains a very good approximation also in the more general case. Calling the linewidth of the LLL Γ LLL we arrive to conclude that the linewidth of the optical density is ap-proximately three times the linewidth of the 1-subband LLL, Γ B ≈ 3 Γ LLL . In [START_REF] Hikami | Anderson localization in a spatially correlated random potential under a strong magnetic field[END_REF] we can find an exact expression for Γ LLL extracted from the density of states of a non-interacting 2-dimensional electron gas in presence of a Gaussian-correlated random potential in the limit of large correlation length ξ c /l B 1 (a very similar result can be derived also for the case of Gaussianly distributed short-range scatters, see [START_REF] Ando | Lineshape of inter-subband optical transitions in space charge layers[END_REF][START_REF] Ando | Theory of Hall Effect in a Two-Dimensional Electron System[END_REF][START_REF] Ando | Electron Localization in a Two-Dimensional System in Strong Magnetic Fields. I. Case of Short-Range Scatterers[END_REF]). A brief summary of the calculation is reported in the Appendix D. Since the resulting density of state is Gaussian, the FWHM reads

Γ LLL ω qw = 2 log(2) π ξ c /l B ξ 2 c /l 2 B + 2 η dis . ( 26
)
The FWHM of the central Gaussian peak of the ISB optical density is thus given by

Γ B, H ω qw ≈ 6 log(2) π ξ c /l B ξ 2 c /l 2 B + 2 η dis . (27)
This formula, that we call Hikami-linewidth, fits extremely well the optical density's FWHM extracted numerically as we can see from Fig. 6, even in the case in which higher Landau level contribute to the optical density.

As expected the width of the central Gaussian scales linearly with the disorder strength ∼ η dis , which is is perfectly captured from Eq. ( 27), see Fig. 6(a). In this plot we also included the corresponding values from the Unuma linewdith Γ qw, U from Eq. ( 23), which gives an estimate of the ISB linewidth without the magnetic field.

It is worth noticing that at this stage the effect of the magnetic field is actually to broaden the ISB transitions and so to worsen the quality factor of our bare QW. In the next section we will see that an opposite result occurs when the QW is embedded in a cavity. For small values of the magnetic field we also found a linear scaling for the width in ∼ ξ c /l B , which means a square-root scaling in the magnetic field intensity ∼ √ B, as can be seen from Fig. 6(b). Quite surprisingly we realise that the Hikami formula gives an accurate estimation of the linewidth even in this regime of moderate correlation length, where ξ c l B . Expanding it at lowest order in ξ c /l B we recover the square-root behaviour in the magnetic field strength. Using Eq. ( 23) we can re-express the small disorder expansion of the Hikami linewidth in a nice and compact expression

Γ B, H ω qw ≈ Γ (1) B, H ω qw = 2 log(2) π 2µ * B B ω qw Γ qw U ω qw . (28)
In Fig. 6(b) we see that both Eqs. ( 27)-( 28) fits very well the numerical data in our regime of interest. It is worth noticing once again the broadening effect of the magnetic field, by comparing the Hikami linewidth to the Unuma linewidth Γ

B, H /Γ qw, U = 2 log(2)/π ω B /ω qw Q qw . Even if the typical magnetic field that we are using through this work is such that ω B /ω qw ∼ 0.1 -0.5, the ISB quality factor is always taken to be Q qw 10, in order to fulfil the condition of small disorder, and so Γ B, H /Γ qw, U 1. In Fig. 6(c) instead the behaviour at larger values of ξ c /l B is reported, and we see how the linewidth saturates for large magnetic field, approaching the asymptotic value

Γ ∞ B H ω qw ≈ 6 log(2) π η dis . (29)

D. Discussion

As we have seen in the previous discussion, there is a large difference in the optical properties of the ISB QW between the non-magnetic and the strong magnetic case. In particular the magnetic field turns the shape of the optical density from a Lorentzian to a Gaussian (a central Gaussian separated from a series of smaller Gaussian side lobes), and sensibly broadens the ISB transition. These differences can be traced back to the localisation properties of the electrons in the disordered potential due to the interface roughness.

Without the magnetic field the electronic states are partially localised due to the interplay of multiple scattering and interference processes on the complex landscape of the disorder potential. The Lorentzian shape of the optical density can be then understood in terms of the finite effective lifetime that the scattering on disorder gives to the electronic coherence.

When the magnetic field is turned on the kinetic energy is completely quenched and we switch to a new regime of strong spatial localisation in Landau levels, where the localisation length is set by the magnetic length l B . In Fig. 7 we can see an example how eigenfunctions at similar eigenenergies localize in a very different way in the two cases. An alternative, possibly more intuitive understanding of this physics can be obtained within a semiclassical picture: in the strong magnetic field regime with l B ξ c , the electrons follow the semi-classical guiding center trajectories along the equipotential lines of the disorder of each subband [START_REF] Huckestein | Scaling theory of the integer quantum Hall effect[END_REF]. Their energy levels are approximately given by ω k ∼ ω B + δU ( r k ), where for each value of the quantum number k labelling the states, one has to pick a different representative realspace position r k along the orbit. Since according to [START_REF] Kurucz | Spectroscopic properties of inhomogeneously broadened spin ensembles in a cavity[END_REF] the disorder potential felt by the two subbands are proportional to each other, the electronic transitions are only between states with the same localization pattern (as described in the last section) and their frequencies are effectively sampling the values of the disorder potential, which follows a Gaussian distribution. This gives a basic intuitive explanation why the resulting optical density is Gaussian-distributed with a linewidth set by the width of the disorder potential. Despite the much larger linewidth, it is the much faster decay of the tails of the Gaussian distribution compared to the ones of a Lorentzian distribution which will be at the heart of our developments in the next sections.

III. INTERSUBBAND POLARITONS AND CAVITY PROTECTION

In this section we show how the dramatic change in the shape of the QW optical density ρ qw (ω) due to the magnetic field that was displayed in Fig. 5 can be used to sensibly improve the properties of the ISB cavity polaritons. It is in fact a well-known fact that a sufficiently strong coupling of a collection of emitters with an overall Gaussian optical density to a single-mode cavity may give rise to polaritonic peaks whose linewidth is only limited by the cavity losses [START_REF] Houdré | Vacuum-field Rabi splitting in the presence of inhomogeneous broadening: Resolution of a homogeneous linewidth in an inhomogeneously broadened system[END_REF]. On contrary, the linewidth of the polariton peaks resulting from an emitter with a Lorentzian-shaped optical density is set by the average of the cavity and emitter linewidths. The physical mechanism underlying these two different behaviours goes under the name of cavity protection [START_REF] Kurucz | Spectroscopic properties of inhomogeneously broadened spin ensembles in a cavity[END_REF][START_REF] Diniz | Strongly coupling a cavity to inhomogeneous ensembles of emitters: Potential for long-lived solid-state quantum memories[END_REF], and is a general feature of polaritonic systems independently of their material realization [START_REF] Manceau | Immunity of intersubband polaritons to inhomogeneous broadening[END_REF][START_REF] Putz | Protecting a spin ensemble against decoherence in the strong-coupling regime of cavity QED[END_REF][START_REF] Zhong | Interfacing broadband photonic qubits to on-chip cavityprotected rare-earth ensembles[END_REF][START_REF] Breeze | Room-temperature cavity quantum electrodynamics with strongly coupled Dicke states[END_REF][START_REF] Baghdad | Spectral Engineering of Cavity-Protected Polaritons in an Atomic Ensemble with Controlled Disorder[END_REF].

A. Magnetic-field-induced cavity protection

We consider here the case in which the magnetic field is strong and exceeds the disorder strength, B > B 0 . Since in a clean sample the polariton frequencies are given by ω ± = ω c ±Ω R /2, when the Rabi frequency becomes comparable with the cyclotron frequency Ω R ∼ ω B , the polariton modes sit in between a pair of Gaussian peaks of the optical density. Here, the optical density has a much smaller value than in the standard non-magnetic case,

ρ qw (ω ± )| B =0 ρ qw (ω ± )| B=0 .
If we think of the system in terms of the clean polariton eigenstates coupled to another continuum of states given by the disorder, we can apply the standard Weisskopf-Wigner theory, expecting that the polaritonic linewidth Γ ± is proportional to the optical density calculated at the polariton frequency, Γ ± ∼ ρ qw (ω ± ). If the optical density drops in correspondence of the polariton frequencies a strong linewidth-narrowing effect is expected.

This behaviour is indeed clearly visible in the numerical results shown in Fig. 8, where we plot the cavity transmission T c (ω) defined in Eq. ( 13) as a function of the transmitted frequency ω and the Rabi frequency Ω R . It is clear that, for a sufficiently large value of B/B 0 and for a Rabi frequency comparable to the cyclotron frequency Ω R ≈ ω B , the polaritonic transmission linewidth is several order of magnitude smaller than the B = 0 case and is only limited by the cavity linewidth.

A further increase of Ω R makes the polaritonic transmission to broaden again as in the weak Rabi frequency case. This happens specifically when the polariton frequencies correspond to one of the side Gaussians of the optical density, that is for Ω R ∼ 2n•ω B , with n = 1, 2 . . . . This behaviour is well captured following the theory developed by Diniz et al. [START_REF] Diniz | Strongly coupling a cavity to inhomogeneous ensembles of emitters: Potential for long-lived solid-state quantum memories[END_REF], where the polaritonic linewidth can be estimated using the following formula

Γ ± ≈ 1 2 γ c + π 2 Ω 2 R ρ qw ω c ± Ω R 2 , ( 30
)
where γ c is the cavity linewidth. This formula holds in the strong coupling regime, when Ω R γ c , Γ qw (see App. E for more details). Having Ω R > γ c , Γ qw means that the polaritonic linewidth Γ ± is determined by the shape of the tails of the QW optical density. If ρ qw has a Lorentzian shape, as in the regular non-magnetic case, we have that π/2 Ω 2 R ρ qw (ω c ± Ω R /2) ≈ Γ qw and the resulting polaritonic linewidth is the average between the cavity and ISB linewidths, Γ ± ≈ (γ c + Γ qw )/2. On contrary, if the QW optical density has a Gaussian shape (or, more generally, decays faster than 1/ω 2 ), we have that π/2 Ω 2 R ρ qw (ω c ± Ω R /2) ≈ 0 and the resulting polaritonic linewidth is only given by the contribution due to the cavity linewidth Γ ± ≈ γ c /2.

In Fig. 9(a) we plot the result of [START_REF] Ullrich | Theory of the Linewidth of Intersubband Plasmons in Quantum Wells[END_REF] in the two cases of B = 0 and B ≈ 4.5 B 0 . When the vacuum Rabi coupling Ω R is comparable to the QW linewidth Γ qw the polaritonic peaks have similar linewidth in both the non-magnetic and magnetic cases, Fig. 9(b) left panel. Instead, when the strong coupling regime is fully reached the mechanism behind the cavity protection sets in: the contribution of the ISB transition to the polariton linewidth drops to zero, as shown in Fig 9(a), and the linewidth is only given by half the cavity linewidth. In Fig. 9(b) central panel, we illustrate the cavity protection effect for a cavity quality factor Q c ≈ 150 close to the typical experimental values. However, as pointed out at the beginning of the section, this is true only when the polariton frequency is located in a gap between two Gaussian peaks, i.e. when Ω R ∼ ω B : when the polariton frequency reaches the next Gaussian peak, at Ω R = 2ω B the linewidth is again broadened to its non protected value, see the right panel of Fig. 9(b). In Fig. 9(c) we plot the quality factors of the lower and upper polaritons defined as

Q ± = ω ± Γ ± , (31)
as a function of the Rabi frequency, in the presence and in the absence of the magnetic field. As in the other plots of Fig 9, we set ω c = ω qw , we take the cavity quality factor as Q c = ω c /γ c = 150 and, when it is on, the magnetic field is set to B/B 0 = 4.5. We extract numerically the polaritonic frequencies ω ± and the linewidths Γ ± as, respectively, the peak frequencies of the cavity transmission T c (ω), and the FWHM of each peak. When the magnetic field is on, the polaritonic quality factor increases sensibly (the light and dark solid lines in Fig. 9(c)), reaching its limiting value

Q ± max = 2 ω c ± Ω R /2 γ c , (32)
marked by the light and dark red dashed lines. For instance, when Ω R /ω qw ≈ 0.3 we have that the polaritonic quality factors increase by a factor between 3 and 4. Of course, the improvement factor would be far more dramatic for higher values of the cavity quality factors Q c as shown in Fig. 8. Such improvements are presently the subject of intense research [START_REF] Abujetas | Spectral and temporal evidence of robust photonic bound states in the continuum on terahertz metasurfaces[END_REF][START_REF] Pisani | Highly resolved ultra-strong coupling between graphene plasmons and intersubband polaritons[END_REF][START_REF] Chalimah | Midinfrared Dispersion Relations in InP-Based Photonic Crystal Slabs Revealed by Fourier-Transform Angle-Resolved Reflection Spectroscopy[END_REF][START_REF] Van Hoof | Unveiling the Symmetry Protection of Bound States in the Continuum with Terahertz Near-Field Imaging[END_REF].

As a final remark, we have to notice that even in the absence of magnetic field both lower and upper polaritonic quality factors Q ± might display a slight increase as a function of the Rabi frequency Ω R (light and dark blue solid lines in Fig. 9(c)). This is a side effect due to the non-perfectly Lorentzian and asymmetric shape of the QW optical density ρ qw (ω). However this is a rather marginal effect if compared to the dramatic linewidth suppression that is induced by the magnetic field.

B. Experimental implementations

In this section we discuss the actual experimental feasibility of the magnetically induced cavity protection for ISB polaritons. In doing this, we will focus on the minimal values of the magnetic field and of the Rabi frequency that are needed to observe the effet. We restrict to devices operating in the THz or mid-infrared regime, mainly focusing on the parameters reported in [START_REF] Colombelli | Perspectives for Intersubband Polariton Lasers[END_REF]. A typical size for a QW in this regime is L z ∼ 8 -40 nm. Using the GaAs electron's effective mass m * ≈ 0.067 m e , where m e is the free electron mass, we have that the energy range for the fundamental ISB transition is given by Here, Bcp is given in units of Tesla, ωqw is given in meV, and ξc in nm. The disorder amplitude is kept constant to ξ0 = 0.3 nm. The effective Bohr magneton corresponding to the effective mass of electrons in GaAs is used, µ * B ≈ 0.86 meV / T. The purple star indicates the position of the setup inspired by [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF] described in the text.

ω qw ∼ 10 -250 meV. (33)
The typical roughness fluctuation scale is around ξ 0 ∼ 0.1 -1 nm and its typical correlation length is ξ c ∼ 1 -10 nm. The typical quality factor of the bare ISB transitions in these devices is in the range Q qw ∼ 10-50.

In order to estimate the typical magnetic field and Rabi frequency needed to quench the linewidth, there are two conditions that we need to satisfy:

1. We need to break the Lorentzian optical density in well separated Gaussians. This is achieved when

ω B > Γ B . (34)
2. We need to have a Rabi frequency that is large enough to overcome the linewidth of the central Gaussian of the optical density and, thus, reach the strong-coupling regime. This condition is maximally fulfilled when the polariton frequency falls in the middle of a gap between two Landau levels. This occurs when

Ω R = ω B . (35)
Because of the very fast decay of the Gaussian tails of the magnetic optical density, the first condition can be considered fulfilled already when ω B = 2 Γ B , from which we define the cavity-protection magnetic field value

B cp = Γ B µ * B . (36)
Using the Hikami-linewidth Eq. ( 27), we can derive an cp as predicted by Eq.( 38) as a function of the intersubband frequency ωqw and the QW quality factor Qqw (in log-scale). Here, B

cp is given in units of Tesla and ωqw in meV. The red hexagon and green circle mark the position of, respectively, typical THz [START_REF] Todorov | Ultrastrong Light-Matter Coupling Regime with Polariton Dots[END_REF][START_REF] Deimert | Realization of Harmonic Oscillator Arrays with Graded Semiconductor Quantum Wells[END_REF][START_REF] Raab | Ultrafast two-dimensional field spectroscopy of terahertz intersubband saturable absorbers[END_REF] and MIR [START_REF] Knorr | Intersubband polariton-polariton scattering in a dispersive microcavity[END_REF] devices.

analytical expression for this quantity,

B cp = eξ 2 c 1 + 144 log(2) π η dis eξ 2 c ω qw 2µ * B 2 1/2 - eξ 2 c . (37)
Using some typical values such as ξ 0 = 0.75 nm, ξ c = 5 nm and L z = 8 nm [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF], and considering the typical effective Bohr magneton µ * B ≈ 0.86 meV / T for the GaAs electron's effective mass, we obtain ω qw ≈ 262 meV, and B cp ≈ 33.4 T, which looks like a quite extreme value for the magnetic field. Already for a slightly smaller value for the QW height fluctuations, for instance ξ 0 = 0.3 nm, we however obtain a much smaller and accessible result B cp ≈ 7.6 T. A complete plot of B cp in unit of Tesla, as a function of the ISB frequency ω qw and the disorder correlation length ξ c in the typical range for THz-MIR ISB polaritons is displayed in Fig. 10 for a fixed value of the QW length fluctuations ξ 0 = 0.3 nm.

Since the values of ξ 0 and ξ c are hardly accessible to direct measurement, it is useful to derive an expression for B cp as a function of quantities that are directly accessible to experiments, like the central ISB frequency ω qw and the bare QW quality factor Q qw . This can be achieved by expanding Eq. [START_REF] Gardiner | Quantum Noise[END_REF] to lowest order in the small disorder limit, and combining it with the Unuma-linewidth in Eq. ( 23)

B (1) cp ≈ 16 log(2) π 1 Q qw ω qw 2µ * B . ( 38
)
This expression is plotted in Fig. 11, where we can clearly identify the regime of interest for THz-MIR devices, in the range of magnetic field between B ∼ 1 -10 T. Note that for the realistic parameters of disorder considered here, the actual values of B

cp and B 0 are quite close.

As a final point, we need to address the second condition for cavity protection, regarding the Rabi frequency Ω R . Combining Eq. [START_REF] Rokaj | Light-matter interaction in the long-wavelength limit: no ground-state without dipole self-energy[END_REF] with Eq. ( 38) we obtain

Ω R = 2µ * B B (1) cp = 16 log(2) π 1 Q qw ω qw . (39)
For the range of quality factors estimated above, we obtain a range for the minimal Rabi frequency necessary to implement the magnetic-induced cavity protection in the order of

Ω R ω qw ∼ 0.05 -0.5. ( 40
)
From this simple estimation it appears that the required Rabi frequency is not far from the ISB transition frequency, Ω R ∼ ω qw , pushing the system from the strongcoupling towards the ultra-strong coupling (USC) regime.

Given the complexities of the theoretical description of the USC regime [START_REF] Frisk Kockum | Ultrastrong coupling between light and matter[END_REF][START_REF] Forn-Díaz | Ultrastrong coupling regimes of light-matter interaction[END_REF] we have chosen to restrict here our attention to the strong coupling physics and we leave a specific investigation of the USC features to a future work.

As a final point, it is important to comment on the role of other processes that, in addition to the scattering on disorder potential, may contribute to the linewidth of the ISB [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF]. While all decoherence channels stemming from static external potentials are tamed by the cavity protection mechanism discussed here for interface roughness, e.g. alloy disorder and ionized impurity, a special attention must be paid to phonon scattering processes. Longitudinal acoustic (LA) phonons typically have low frequencies, much lower than the intrinsic broadening of polaritons: as such, they can be considered as quasi-static and the virtually elastic LA phonon scattering thus fall in the same category of static potentials.

Longitudinal optical (LO) phonons have instead a much higher energy E LO on the order of a few ten meV in typical materials and their contribution to the linewidth is small but sizable in many cases, giving a lower bound to the achievable linewidth on the order of 1meV for transitions in the MIR according to [START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF]. Still there exist regimes, e.g. THz ISB transitions with ω 21 < E LO and a low electron density such that E F < E LO for which LO phonon emission is kinematically not allowed and the lower bound disappears.

IV. CONCLUSIONS

In this work, we have proposed and characterized a strategy to dramatically improve the quality factor of intersubband polaritons in semiconductor-based devices.

By applying a strong magnetic field perpendicular to the quantum well plane, the dominant Lorentzian broadening of the intersubband transition due to scattering of electrons onto interface roughness disorder turns into a Gaussian one due to the strong electron localization in the disorder potential. When strongly-coupled to the cavity mode, a cavity protection mechanism sets in which removes the Gaussian linewidth, leaving only the cavity contribution to the polariton linewidth.

In combination with higher-Q cavity configurations, our proposal has the potential to lead to polariton devices with unprecedented performances. Based on available experimental evidence, we anticipate that a narrower polariton linewidth will be a game-changing step in view of technological applications of intersubband polaritons, including nonlinear polariton devices [START_REF] Cominotti | Theory of coherent optical nonlinearities of intersubband transitions in semiconductor quantum wells[END_REF] and polariton lasing [START_REF] Knorr | Intersubband polariton-polariton scattering in a dispersive microcavity[END_REF][START_REF] Colombelli | Perspectives for Intersubband Polariton Lasers[END_REF]. On a longer perspective, our proposal highlights intersubband polariton physics in the presence of magnetic fields as a novel arena where to study the interplay between light-matter interaction with quantum Hall physics [START_REF] Smolka | Cavity quantum electrodynamics with many-body states of a two-dimensional electron gas[END_REF][START_REF] Ravets | Polaron Polaritons in the Integer and Fractional Quantum Hall Regimes[END_REF][START_REF] Knüppel | Nonlinear optics in the fractional quantum Hall regime[END_REF][START_REF] Appugliese | Breakdown of topological protection by cavity vacuum fields in the integer quantum Hall effect[END_REF][START_REF] Ciuti | Cavity-mediated electron hopping in disordered quantum Hall systems[END_REF][START_REF] Rubio | A new Hall for quantum protection[END_REF][START_REF] Rokaj | Quantum Electrodynamical Bloch Theory with Homogeneous Magnetic Fields[END_REF].

the ISB transitions. We have

i∂ t a = ω c -i γ c 2 a + Ω R 2 λ Λ λ s λ -+ γ 2 α in (t) + b c , (A3) i∂ t s λ -= (ω qw + ω λ )s λ --Ω R Λ λ s λ z a+ + i2κs λ z s λ --2s λ z b λ , (A4) i∂ t s λ z = Ω R 2 Λ * λ s λ + a -Λ λ s λ -a † + -i2κs λ z -iκ + b λ s λ + -b † λ s λ -. (A5)
Here γ c is the cavity loss rate, κ is the ISB transition relaxation rate, b c and b λ are respectively the cavity and ISB transition quantum noise operators. α in (t) = α 0 in e -iωt is the coherent input field , assumed to be monochromatic at a given frequency ω. The input α in , reflected α r , transmitted α t and cavity field a are related by the following input-output formulas

α r = α in + γ c /2a, α t = γ c /2a. (A6)
We now consider the mean-field and weak driving regime. We then replace the operators a, s λ -, s λ z in Eqs. (A3)-(A4)-(A5) with their respective expectation values (for simplicity in the notation we still keep the same notation for the operators and their mean value). The quantum noise operators drop out since that at temperature T = 0 they have zero mean value. Because of the weak drive we can also approximate s λ z ≈ -1/2. In the rotating frame at the input frequency ω, the mean-field Langevin equations take the form of a driven dissipative system of coupled harmonic oscillators

i∂ t a = ω c -ω -i γ c 2 a+ Ω R 2 λ Λ λ s λ -+ γ 2 α in , (A7) i∂ t s λ -= (ω qw + ω λ -ω)s λ -+ Ω R 2 Λ λ a -iκs λ -. (A8)
Solving these equations for the steady state, and defining the optical transmission through the cavity as

T c (ω) = α t α in , (A9)
we immediately arrive to Eq. ( 13).

roughness δL is a random variable, in general with nontrivial spatial correlations. We can thus write

δL( r)δL(0) = ξ 2 0 C( r), (C3)
where ξ 0 quantifies the magnitude of the fluctuations of the QW thickness and C( r) is an arbitrary correlation function normalised to have C(0) = 1/(2π). We can then introduce the adimensional random variable ∆( r) = δL( r)/ξ 0 .

We then rescale the position and the momenta as

r -→ r • l qw p -→ p • /l qw (C4)
in terms of an equivalent harmonic-oscillator length within the QW

l qw = m * ω qw 1/2 . ( C5
)
Defining a = (-y, x)/2, and considering the usual cyclotron frequency as ω B = eB/m * , we have

H n = ω qw ( p -η B a) 2 2 + g n η dis ∆( r) , (C6)
where we introduced the two adimensional energy magnitudes

η B = ω B ω qw η dis = ∂ L ε 1 (L qw )ξ 0 ω qw , (C7)
and the numerical coefficients

g n = ∂ L ε n (L qw ) ∂ L ε 1 (L qw ) . ( C8
)
In the specific case of an infinite box potential g n = n 2 , while in the harmonic oscillator case g n = n -1/2, with n = 1, 2 . . . . In our numerical calculations, the adimensional disorder is generated by sampling ∆( r) from a Gaussian distribution with unit variance at every position r. We then take its Fourier transform and we impose a Gaussian cut-off

∆( k) -→ ξ c √ 2 e -ξ 2 c k 2 /8 ∆( k). (C9)
After transforming back to real-space, we are left with the desired correlator

∆( r)∆( r ) = e -| r-r | 2 /ξ 2 c 2π . ( C10
)
We then proceed to diagonalize the rescaled Hamiltonian Eq. (C6) on a given basis. When the magnetic field is not present B = 0 we consider a finite box of lengths L x , L y whose basis wave functions are

ϕ nx,ny ( r) = 2 L x L y sin( πn x L x x) sin( πn y L y y), (C11)
with n x , n y = 1, 2 . . .. Same results are obtained considering a periodic system and considering a planewave eigenbasis.

In the presence of a finite magnetic field B = 0, we first need to fix a gauge for the vector potential. For instance, in the symmetric gauge where a = (-y, x)/2. In order to simulate an infinite system with open boundary conditions we use the basis set

ϕ k ( r) = 1 √ 2πl B ! k! ξ k-e -|ξ| 2 /2 L k-(|ξ| 2 ), (C12)
where ξ = (x + iy)/( √ 2l B ) and L α (x) are the generalised Laguerre polynomials. Since these states are concentric circles, centered in the origin, we include as many states as possible to fill the spatial extension of the system, cutting the basis right before touching the border of the numerical space grid. If the total extension of all basis states L basis is such that L basis l B , ξ c , averaging over many realisations is then approximately similar to consider a very large system, where every realisation is a smaller patch of the whole system.

In order to compute the matrix elements on a given basis we consider a two dimensional spatial grid N x ×N y , where we typically use a number of grid points in the range N x = N y = 100 -150, with a grid step in the range ∆r/L z = 0.1 -0.5, dependently from the chosen values of ξ c and l B . All final results are typically averaged over N dis ∼ 20 -100 realisations.

Notice that we do not include the edge states in our calculation. This is motivated by the quantitative smallness of their contribution, which is completely negligible in an extended system. Indeed in a very large system the edge modes represents a much smaller fraction of the whole system, and, since the dipole moment along z is the same for every electron, localised or non-localised, their contribution is negligible with respect to the total number of bulk states. In this section we briefly review the calculation of the disordered density of states in the lowest Landau level (LLL) contained in [START_REF] Hikami | Anderson localization in a spatially correlated random potential under a strong magnetic field[END_REF]. For simplicity we work in adimensional units where ω qw = 1.

The calculation starts from the identity in Eq. (B6). What we need to calculate is the Green's function for the disordered Schrödinger equation, projected in the LLL and averaged over the disorder. To do so we introduce a quantum field theory representation of the Schrödinger Green's function in terms of complex scalar field path integral [START_REF] Zee | Quantum Field Theory in a Nutshell[END_REF] G(x, y, ω) = x| 1 ω -H |y = 1 i

DφDφ * e i d 2 x [φ(x ) * (ω-H)φ(x )] φ(x)φ * (y) DφDφ * e i d 2 x [φ(x ) * (ω-H)φ(x )] .

(D1)

We then use the relation between path integral and functional determinant to transform the Bosonic path integral in the denominator into a Grassmannian path integral following the so-called super-symmetric approach [START_REF] Brézin | Density of states in the presence of a strong magnetic field and random impurities[END_REF] 1 DφDφ * e i d 2 x [φ(x ) * (ω-H)φ(x )] = det

(D3)

We now take the average over the disorder. Since we are considering a Gaussian disorder, we can safely apply the second cumulant expansion, exp -i d 2 x δU (x )f (x ) = = exp -1/2 d 2 x 1 d 2 x 2 δU (x 1 )δU (x 2 )f (x 1 )f (x 2 ) .

(D4)

where • • • indicates the average over many disorder realisations. We then introduce the disorder correlator as δU (x 1 )δU (x 2 ) = C(x 1 -x 2 ). In the case of infinite range correlations (ξ c → ∞) we have C(x 1 -x 2 ) = C 0 , and the averaged Green's function reads To get the density of states, we need to set x = y, integrate over the whole area A of the system and finally normalise over the total electron number. Since

k |ϕ 0 k (x)| 2 = 1/(2πl 2
B ) we get a factor A/(2πl 2 B ) from the spatial integration. In a finite system the total number of states in the LLL is N LLL ∼ A/(2πl 2 B ) [START_REF] Huckestein | Scaling theory of the integer quantum Hall effect[END_REF]. Completing the square, performing the Gaussian integral in Eq. (D6) and using Eq. (B6) we finally find the density of states for the LLL in the infinite correlation-length approximation

ρ ∞ LLL (ω) ≈ 1 η dis e -πω 2 η 2 dis , (D9)
where we considered that C 0 = η 2 dis /(2π). From this distribution we extract the Hikami linewidth in the infinite correlation-length case Γ ∞ B H defined in Eq. [START_REF] Nikonov | Collective Intersubband Excitations in Quantum Wells: Coulomb Interaction versus Subband Dispersion[END_REF]. In order to include information about the finite correlation length, we expand the correlator in series as C(x) ≈ C 0 (1 -x 2 /ξ 2 c ). In Eq. (D5) we need to add to the argument of the exponential in the last line the following term

C 0 2 d 2 x 1 d 2 x 2 |x 1 -x 2 | 2 ξ 2 c
(φ * (x 1 )φ(x 1 ) + η(x 1 )η(x 1 )) × × (φ * (x 2 )φ(x 2 ) + η(x 2 )η(x 2 ))

≈ 2l 2 B ξ 2 c k (a * k a k + b * k b k ). 2 . ( D10
)
We see that this term contributes only by a shift of C 0 which is replaced by C 0 -2C 0 l 2 B /ξ 2 c . We then re-sum all the terms similar to this one from the higher order correlator expansion considering that C 0 -2C 0 l 2 B /ξ 2 c ≈ C 0 /(1 + 2l 2 B /ξ 2 c ). In this way, we arrive to the LLL density of states in the case of a finite correlated disorder potential .

ρ LLL (ω) ≈ 1 √ 2πγ LLL e -ω 2

(D12)

The FWHM reported in Eq. ( 26) is then equal to Γ LLL = 2 2 log(2) γ LLL .

  FIG.1.a) Sketch of the system. A semiconductor twodimensional quantum well is inserted in a metallic cavity. The electric field of the cavity mode oscillates along the zdirection, perpendicular to the quantum well plane. The electronic intersubband transitions in the quantum well give rise to an electronic polarization directed along the same direction of the cavity electric field, realising a strong light-matter coupling regime between the cavity and the quantum well. b-c) Schematic view of the electronic subbands in the disordered quantum well in absence/presence of the external magnetic field perpendicular to the quantum well plane. The lowest optically-active excitations are transitions between the lowest two subbands due to the electrons below the Fermi energy EF . Differently from the clean case, these transitions are not only fully vertical and each electron can jump to different eigenstates of the upper subband.

FIG. 2 2 .

 22 FIG. 2. a) Spatial profile of a single realisation of the disorder potential δUn( r) for the n = 1 subband. b) Fourier transform of the same disorder potential δ Ũn( k) as a function of the wavenumber k. Disorder parameters: η dis ≈ 0.06, ξc/lqw = √ 2. The numerical methods and the integration parameters are discussed in the Appendix C.

3 FIG. 3 . 2 . 16 FIG. 4 .

 332164 FIG. 3. a) Optical density ρqw(ω) (blue solid line) and Lorentzian fit (black dashed line). b) Optical density ρqw(ω) calculated for increasing values of the Fermi energy EF . Disorder parameters: η dis ≈ 0.06, ξc/lqw = √ 2. In a) the Fermi energy is fixed to EF / ωqw ≈ 0.6, while in b) a range of values is used as indicated in the legend.

FIG. 5

 5 FIG.5. a) Plots of the optical density ρqw(ω) for increasing values of the magnetic field B/B0 = 1, 4.5, 7 (red solid lines). Blue lines show the same quantity in the absence of magnetic field. Disorder parameters: η dis = 2/15 ≈ 0.13, ξc/lqw = 0.5. Fermi energy EF / ωqw = 0.5. For this choice, the quality factor of the intersubband transition in the absence of magnetic field is numerically found to be Qqw ≈ 40. b) Plots of the optical density ρqw(ω) with fixed magnetic field B/B0 = 2.5 for increasing values of the Fermi energy EF = 0, ωB, 2ωB. Here the lowest Landau level is set to be at zero energy, such that the Fermi energy falls in the middle of the disorder-broadened 1 = 0, 1, 2 Landau level. The vertical dashed line highlights the asymmetry of the central Gaussian. Disorder parameters: η dis = 0.08, ξc/lqw = 1.

15 FIG. 6

 156 FIG.6. a) Gaussian linewidth (FWHM) ΓB of the central peak of the optical density as a function of the disorder strength η dis for a constant value of the magnetic field B. The red squares/green dots are the numerical data corresponding to increasing values of B, that is B/ B0 = 0.2, B/ B0 = 0.6. The disorderless reference magnetic field is defined as B0 = ωqw/(2µ * B ). The red/green dashed lines are given by the Hikami linewidth ΓB, H/ωqw while the black solid line is the Unuma linewidth Γqw, U/ωqw. Disorder parameters: ξc/lqw = 0.5. b) Gaussian linewidth ΓB (purple dots) extracted from numerical simulations of the optical density as a function of the magnetic field B. The dashed red line is given by the Hikami linewidth ΓB, H/ωqw, while the dashed-dotted line is given by the small disorder expansion of the Hikami linewidth Γ

4 FIG. 7 .

 47 FIG. 7. Illustrative examples of the real-space wavefunction of a n = 1 subband eigenstate in the presence of disorder in the B = 0 (left) and B = 0 (right) cases.

2 FIG. 8 .

 28 FIG. 8. Logarithmic colorplot of the cavity transmission as a function of the Rabi coupling ΩR and the external incident frequency ω for various values of the magnetic field B in the resonant regime ωc = ωqw. Disorder parameters as in Fig. 5: η dis = 2/15 ≈ 0.13, ξc/lqw = 0.5. Fermi energy EF / ωqw = 0.5.The ISB quality factor (without cavity) in the three panels is numerically found to be Qqw ≈ 40, 12, 10, from left to right. An extreme value of the cavity quality factor Qc = 10 5 is taken for illustrative purposes, more realistic values will be considered in the next figures. In physical units, for a QW of thickness Lz = 26 nm and m * ≈ 0.067 me, these adimensional parameters correspond to ωqw ≈ 25 meV and a reference magnetic field B0 ≈ 1 T. In the central panel we thus have B ≈ 4.4 T while in the right panel B ≈ 6.9 T.

FIG. 9

 9 FIG.9. a) Plot of the electronic contribution to the polariton linewidth estimated using Eq. (30) as a function of Rabi coupling ΩR, left (right) parts corresponding to the upper (lower) polariton. The blue and red lines refer to the B = 0 and B = 4.5 B0 cases, respectively. b) Three examples of the cavity transmission Tc(ω) given by Eq. (13) as a function of the incident frequency ω in the presence/absence of magnetic field (blue and red lines, respectively) for increasing values of the Rabi coupling ΩR (from left to right) as indicated by he dashed, dotted and dashed-dotted lines in panel a). c) The quality factor of the lower (light) and upper (dark) polaritons Q± = ω±/Γ± as a function of the Rabi frequency ΩR. Blue and red lines refer to the cases in the absence/presence of magnetic field B/B0 = 0, 4.5. The black dotted and black dashed-dotted lines indicate the usual upper/lower polariton quality factors obtained averaging between cavity and ISB linewidths, Q ±, std = 2(ωc±ΩR/2)/(γc+Γqw). The light/dark red dashed lines indicate the upper bound to the polariton quality factor set by the cavity losses as defined in Eq.[START_REF] Kempa | Electron-electron scattering in strong magnetic fields in quantum well systems[END_REF]. Same system parameters as in Fig.8except for the realistic value Qc = 150 of the cavity quality factor.

FIG. 10 .

 10 FIG.10. Log-scale contour-plot of the cavity-protection magnetic field Bcp as defined in[START_REF] Gardiner | Quantum Noise[END_REF] as a function of the intersubband frequency ωqw and the disorder correlation length ξc.Here, Bcp is given in units of Tesla, ωqw is given in meV, and ξc in nm. The disorder amplitude is kept constant to ξ0 = 0.3 nm. The effective Bohr magneton corresponding to the effective mass of electrons in GaAs is used, µ * B ≈ 0.86 meV / T. The purple star indicates the position of the setup inspired by[START_REF] Unuma | Intersubband absorption linewidth in GaAs quantum wells due to scattering by interface roughness, phonons, alloy disorder, and impurities[END_REF] described in the text.

10 FIG. 11 .

 1011 FIG. 11. Log-scale contour-plot of the cavity-protection magnetic field B (1)

  Appendix D: Lowest Landau level disordered density of state

1 i

 1 [ω -H] = = DηDηe i d 2 x [η(x )(ω-H)η(x )] . (D2)Here η(x) is an anticommuting Grassmann field. Considering our Hamiltonian composed by kinetic energy T and disorder potential δU termsx| 1 ω -H |y = DφDφ * DηDη φ(x)φ * (y)× × exp i d 2 x φ * ω -T φ + η ω -T η × × exp -i d 2 x δU (x ) (φ * φ + η η) .

2 d 2 x 2 .

 222 DφDφ * DηDη φ(x)φ * (y)×× exp i d 2 x φ * ω -T φ + η ω -T η × × exp -C 0 (φ * (x )φ(x ) + η(x )η(x )) the fields over the LLL considering φ(x) = k ϕ 0 k (x)a k and η(x) = k ϕ 0, k b k , where ϕ ,kis the wave function of the k-th state in the -th Landau level. We also insert in the integral an identity operator using a delta function and we shift ω → ω -ω B /2. The Green's function then readsx| 1 ω -H |y = 1 i +∞ -∞ dσ e -C 0 2 σ 2 +iωσ D(x, y, σ) (D6) with D(x, y, σ) = k ϕ 0 k (x)ϕ * 0 k (y)I k (σ),(D7)andI k (σ) = Π k d 2 a k d 2 b k a * k a k δ σk (a * k a k + b * k b k ) = lim ε→0 dt 2π Π k d 2 a k d 2 b k a * k a k e itσ-i(t-iε) k (a * k a k +b * k b k ) the step function Θ(σ).

ACKNOWLEDGMENTS

We are grateful to Giacomo Scalari, Ivan Amelio, Alberto Nardin and Jacopo Nespolo for fruitful discussions. We acknowledge financial support from the European Union FET-Open grant MIR-BOSE (737017).

Appendix A: Input-output theory of the cavity plasma Hamiltonian in the weak excitation regime Consider the cavity-plasma Hamiltonian defined in Eq. [START_REF] Knorr | Intersubband polariton-polariton scattering in a dispersive microcavity[END_REF]. Since we are looking only at its low energy excitations we can truncate the electronic Hilbert space keeping only the unperturbed Fermi sea state |F S and its lowest single-electron-hole-pair excitations between the first two subbands, given by

where k ≤ k F , with k F the Fermi momentum. All the Fermionic operators in Eq. ( 10) are then replaced by Pauli matrices, each of them representing a transition between the Fermi sea and one electron-hole state ( = 1)

Here every couple λ = (k, k ) identify a single two-level system and s λ z , s λ -, s λ + are the usual spin-1/2 operators. The range of the vectorial index λ is limited to the semirectangular area (-∞,

We then derive the usual quantum Langevin equations [START_REF] Gardiner | Quantum Noise[END_REF] assuming a two-sided cavity and a generic bath for Appendix B: Electronic response and subband Green's function

The electronic optical response defined in Eq. ( 14) can be rewritten within a Green's function formalism. We start from the initial formula

Using the Sokhatsky identity

we can write

where the Green's function is defined in the Lehman representation as

is the Fermi distribution. The QW optical spectral density is then given by

where

is the local spectral density of a single n-subband.

Appendix C: Technical details on the numerical calculation

Here we consider the Hamiltonian as defined in Eq. ( 5), including the presence of the homogeneous magnetic field via the symmetric-gauge vector potential A( r) = B/2(-y, x, 0),

Here, the disorder potential is

in terms of the energy ε n (L qw ) of the electronic state trapped in the QW in the z-direction. The QW interface

Here we summarise the results contained in [START_REF] Diniz | Strongly coupling a cavity to inhomogeneous ensembles of emitters: Potential for long-lived solid-state quantum memories[END_REF] regarding the cavity protection effect.

The poles of the cavity transmission Eq. ( 13) give the system's resonant frequencies and their linewidth, through the equation T -1 c (ω) = 0. In a perfectly clean system or when the disorder potential is identical for the two subbands, we have Λ k k = δ kk . We thus obtain the standard cavity transmission modified by the coupling to the QW

where

is the resonant response of the QW. Assuming the lightmatter resonance condition ω c ∼ ω qw in the strong coupling regime Ω 2 R /(γ c κ) 1, we have that the cavity transmission has two resonant peaks, corresponding with the polaritonic frequencies ω ± ≈ ω c ± Ω R /2. The linewidth of these two peaks is the same and is given by the average between the cavity and the QW linewidths Γ ≈ (γ c + κ)/2.

In a situation of small disorder we still expect to have the two polariton peaks around the frequencies ω ± but with a modified linewidth. In order to derive this new linewidth we split the QW optical response in real and imaginary part χ(ω) = χ R (ω) + iχ I (ω), considering that χ I = πρ(ω). We also work in a rotating frame, such that ω -ω c -→ ω, with ω c = ω qw . The poles of the cavity transmission are then given by

Assuming the strong coupling regime, where Ω R is much larger than the linewidth of ρ(ω)

and we approximate ρ(ω) ≈ ρ(ω ± ) , (E5) From these relations we immediately obtain the approximated pole equation, valid only in the proximity of

Solving it we find

where the polaritonic linewidth is given by


N Lesparre
email: lesparre@unistra.frackerer@unistra.fr

P Ackerer

J F Girard
email: jf.girard@unistra.fryounes@unistra.fr

A Younes

T Mara

Global sensitivity analysis of Magnetic Resonance Sounding simulated responses to the hydraulic parameters of a catchment hydrogeological model

Keywords: MRS, hydrological model, sensitivity analysis, catchment scale

INTRODUCTION

Magnetic resonance sounding (MRS) is a geophysical method directly sensitive to the underground water content. Beyond the estimate of the water content state at a given time, the method has been successfully applied to constrain the parameters of hydrological models. At the scale of one MRS acquisition setup (a few tens of meters), [START_REF] Costabel | Noninvasive estimation of water retention parameters by observing the capillary fringe with magnetic resonance sounding[END_REF] explored how MRS data can help estimating the water retention parameters. At a similar scale, [START_REF] Herckenrath | Coupled hydrogeophysical inversion using time-lapse magnetic resonance sounding and time-lapse gravity data for hydraulic aquifer testing: Will it work in practice?[END_REF] studied the relevance of coupled MRS and gravity measurements repeated over time to retrieve the aquifer characteristics during a pumping test. [START_REF] Boucher | Constraining groundwater modeling with magnetic resonance soundings[END_REF] used MRS measurements to constrain the hydraulic conductivity at saturation of a hydrological model applied at a catchment scale. Similarly, the fitting of MRS data covering a catchment with estimates from a hydrological model highlighted the necessity to divide the aquifer in two vertical layers with a deeper less porous medium [START_REF] Lesparre | Magnetic resonance sounding measurements as posterior information to condition hydrological model parameters: Application to a hard-rock headwater catchment[END_REF]. More recently, the acquisition of timelapse MRS measurements helped modelling the water flow in an unsaturated medium [START_REF] Legchenko | Time-lapse magnetic resonance sounding measurements for numerical modeling of water flow in variably saturated media[END_REF].

Applications of time-lapse MRS acquisitions are progressing and a recent study showed its interest to monitor water infiltration during a rainy event [START_REF] Mazzilli | Surface Nuclear Magnetic Resonance Monitoring Reveals Karst Unsaturated Zone Recharge Dynamics during a Rain Event[END_REF]. In this paper, we explore how such time-lapse MRS simulated measurements are sensitive to the parameters of a hydrological model in order to discriminate which parameters MRS could help constraining in practice. For this purpose, we perform a global sensitivity analysis (GSA) of MRS data to the hydraulic parameters of a hydrogeological model applied on the Strengbach catchment (Fig. 1). We applied the measured meteorological forcing during a six-month period and analyse the variability of the MRS estimate at two stations presenting contrasting hydrological behaviours. We perform the same analysis with water table depth estimates to illustrate the complementarity between MRS and piezometric head data.

METHODS

The hydrological model applied covers the whole Strengbach catchment of 0.8 km², located in the Vosges mountains (North-East of France). The catchment climate is temperate of oceanic mountainous type with a mean temperature of 6°C. 90% of the catchment is covered by forest with 80% of spruces and 20% of beeches. The underlying bedrock is of granitic type and the soils are coarsely grained, sandy and rich in gravels and pebbles.

SUMMARY

We apply a global sensitivity analysis of magnetic resonance sounding (MRS) data to the hydraulic parameters of a hydrological model at a catchment scale. The methodology used provides a quantitative estimate of the proportion to which the explored parameters influence the simulated MRS response. The sensitivity analysis is performed through a variance-based approach that provides indicators to infer how the MRS sensitivity to the different parameters might vary in space and time. We also apply the sensitivity analysis to water table depth estimates. We show that MRS measurements complement fruitfully classical water table depth acquisitions. Indeed, MRS is sensitive to the model water content at saturation while this parameter does not influence significantly the water table depth estimates.

MRS 2021, 26-28 October, Strasbourg, France. https://mrs2021.sciencesconf.org We use the Normally Integrated Hydrological Model (NIHM) to simulate the water flow in the underground and at the surface through a coupling between the two compartments [START_REF] Weill | A lowdimensional subsurface model for saturated and unsaturated flow processes: ability to address heterogeneity[END_REF][START_REF] Weill | A lowdimensional subsurface model for saturated and unsaturated flow processes: ability to address heterogeneity[END_REF][START_REF] Jeannot | A low-dimensional integrated subsurface hydrological model coupled with 2-D overland flow: application to a restored fluvial hydrosystem (Upper Rhine River-France)[END_REF]. We run 2500 synthetic models with randomly selected hydrological model parameter values. We focus the analysis on, the water content at saturation s  , the hydraulic conductivity at saturation s K and the  parameter of the Van Genuchten equation. The aquifer is divided in two vertical layers corresponding to the soil and the saprolite, the last one being less porous. The soil and the saprolite thicknesses are defined through a statistical analysis of seismic profiles covering the catchment (Fig. 1). Geostatistical models are generated to reproduce the measured seismic velocity at the profile location and their statistical variability at the catchment scale. The soil and saprolite bottom interfaces are then determined by fixing a seismic velocity threshold of 700 m/s and 2000 m/s, respectively. The soil and saprolite thicknesses vary then normally in each element NIHM. s  and  vary uniformly, while s K follows a uniform distribution in the logarithmic domain. The variation ranges of those three parameters are summarized in Table 1.

 (m -1 ) 0.1 -2 0.1 -2
The models run from the 1 st of June 2012 to the 31 st of May 2013, the first 6 months are used to "warm-up" the simulations and the GSA is applied from the 1 st of December 2012. The studied period shows several distinct rainy events (Fig. 2), so we can examine the influence of different water content states on the MRS and water level depth variability. The measured precipitation and the estimated evapotranspiration are provided as a model input, such data come from the OHGE meteorological station (Fig. 1). The temperature is also given to determine the nature of the precipitation (rain or snow) and in case of snow, the infiltration is delayed depending on the temperature and on the previous height of snow accumulated.

NIHM provides as outputs the flow at the outlet, the water head at each MRS and piezometric stations and daily maps of the underground water content used to estimate the MRS signal as described in [START_REF] Lesparre | Magnetic resonance sounding measurements as posterior information to condition hydrological model parameters: Application to a hard-rock headwater catchment[END_REF]. The global sensitivity of each data type to the parameters explored is then analysed through variance-based sensitivity indices [START_REF] Mara | Application of global sensitivity analysis of model output to building thermal simulations[END_REF]. Such indices distinguish the contribution of each parameter to the signal total variance [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] We compute the Sobol' indices through Polynomial Chaos Expansion (PCE; [START_REF] Wiener | The homogeneous chaos[END_REF]. Such representations are indeed particularly suited to perform variance-based GSA [START_REF] Younes | Sensitivity and identifiability of hydraulic and geophysical parameters from streaming potential signals in unsaturated porous media[END_REF]. The Sobol indices are directly computed from the PCE coefficients, without running the surrogate models, at a daily rate at each measurement point and for each data type.

[ [ ( ) | ]] [ ( ] i i V E y X x S V y X  . ( 1

RESULTS

The contribution of the explored parameters to MRS signal variance as function of time is represented for station 7 and 22 Fig. 3. For sake of concision we limit here the analysis to those stations since they present a distinct hydrogeological behaviour. Station 7 is situated in a steep region that favours water drainage. Conversely, station 22 is located on a relatively flat region more suitable for water storage. We note that the variance of the MRS signal is much stronger at station 22 than at the 7 th . We observe that at both stations the MRS signal is mainly sensitive to the soil s  and the saprolite . Both stations are also sensitive to the saprolite s  , station 7 is more significantly sensitive to the soil  value, while station 22 is sensitive, at times, to the saprolite  . Station 22 is also much more sensitive to the soil thickness than station 7. We also observe that a large part of the signal variance is related to contributions concerning several parameters in interaction since the total variance (blue lines) is well above the respective contributions of each explored parameter. The application of the same GSA to water table depth underlines the complementarity between MRS and piezometric data. In particular, the water table depth variance is not significantly impacted by the s  parameter (Fig. 4). Under both stations, the water table depth variance is mostly sensitive to the saprolite s K , and to the  parameter of both the soil and saprolite. We note that the variance at station 7 is very smooth. At station 22, the variance of the water table depth presents much more irregulrities corresponding to the different precipitation events. We note that the impact of saprolite s K on the water table depth signal increases with the spring arrival.

As for the MRS signal, a part of the water table depth total variance is related to interactions between several parameters. analysis of these results invites to deepen the study to distinguish which coupled parameters impact both the MRS and piezometric signal. A study of the whole measurement stations variance would also help understanding the respective role of the different parameters to the MRS and water table depth signals. The GSA performed here encourages time-lapse MRS measurements at different places of a catchment to confirm our findings by including aspects such as the measurements noise that could vary in space and time.

Figure 1 .

 1 Figure 1. Distribution of MRS measurement stations and seismic profiles acquired on the Strengbach catchment.

Figure 2 .

 2 Figure 2. Rain minus the evapotranspiration measured on the Strengbach catchment used as the meteorological forcing inserted in NIHM (top). Flowrate measured at the catchment outlet (below).

  26-28 October, Strasbourg, France. https://mrs2021.sciencesconf.org However, we note that this sensitivity varies with time, notably concerning the s K parameter. At station 7, the contribution of s K is higher in winter during the drainage periods that follow high flow events. At station 22, s K shows also a higher contribution in a period of drainage, but in spring when precipitations are less intense and more scattered across time.

Figure 3 .

 3 Figure 3. Variance of the MRS signal as function of time estimated at station 7 (top) and 22 (bottom). The blue line represents the total signal variance while the coloured areas indicate the respective contribution of the explored parameters.

Figure 4 .

 4 Figure 4. Variance of the water table depth as function of time estimated at station 7 (top) and 22 (bottom). The blue line represents the total signal variance while the coloured areas indicate the respective contribution of the explored parameters.

Table 1 :

 1 Range of variations of the explored parameters.

ACKNOWLEDGMENTS

Meteorological and flow rate data collection was funded by the Observatoire Hydro-Géochimique de l'Environnement (OHGE), which is financially supported by CNRS/INSU France and the University of Strasbourg. OHGE is part of the OZCAR research infrastructure (http://www.ozcar-ri.org). The authors would like to acknowledge the ANR HYDROCRIZSTO-15-CE01-0010-02 project. Computing time was provided by the HPC-UdS.