{'fraction_non_alphanumeric': 0.02396099373113536, 'fraction_numerical': 0.009890875319247736, 'mean_word_length': 4.848859315589354, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theory of coherent optical nonlinearities of intersubband transitions in semiconductor quantum wells

We theoretically study the coherent nonlinear response of electrons confined in semiconductor quantum wells under the effect of an electromagnetic radiation close to resonance with an intersubband transition. Our approach is based on the time-dependent Schr\"odinger-Poisson equation stemming from a Hartree description of Coulomb-interacting electrons. This equation is solved by standard numerical tools and the results are interpreted in terms of approximated analytical formulas. For growing intensity, we observe a red-shift of the effective resonance frequency due to the reduction of the electric dipole moment and the corresponding suppression of the depolarization shift. The competition between coherent nonlinearities and incoherent saturation effects is discussed. The strength of the resulting optical nonlinearity is estimated across different frequency ranges from Mid-IR to THz with an eye to on-going experiments on intersubband polariton condensation and speculative quantum optical applications such as single-photon emission.


I. INTRODUCTION
Intersubband (ISB) transitions in semiconductor quantum wells (QW) play a crucial role in a number of optoelectronic devices across a wide range of wavelengths, from the visible down to the IR and the THz ranges 1 . A most celebrated example is the quantum cascade laser, which is one of the most widespread semiconductorbased sources of coherent radiation for the Mid-IR and THz ranges of the electromagnetic spectrum 2 . Nonlinear optical effects in these wavelength regions are also attracting a great interest, in particular for what concerns the realization of passively mode-locked pulsed laser sources 3,4 , optical combs 5,6 as well as switching, modulation and harmonic generation when combined with meta-surfaces [7][8][9] .
Many among these developments are based on incoherent optical nonlinearities which result from a saturation mechanism due to the shelving of electrons into optically dark states and the consequent reduction of the effective oscillator strength of the transition [10][11][12][13] . Since they stem from an incoherent dynamics and have a relatively slow time-scale determined by the decay rate of the dark excitations, these nonlinearities can hardly result into coherent processes. On the other hand, coherent nonlinearities associated to ultrafast processes like Rabi oscillations of a two-level transition has been widely studied under strong pulsed illuminations [14][15][16] . The goal of this work is to contribute building a general theoretical picture of the nonlinear response of ISB transitions encompassing the two regimes.
The simplest theoretical description of the nonlinear optical response of electrons in quantum wells is based on master equations that only include a few discrete electronic states 9,16,17 . While this is accurate in the low electronic density limit, at higher densities Coulomb interactions start playing a significant role deforming and mixing the single-electron states. At the level of linear response to weak beams, they are responsible for the depolarization shift of the ISB transition [18][19][20][21] . For growing light intensities, many intriguing phenomena have been anticipated 17 , in particular the saturation of the optical transition was predicted to give to a corresponding reduction of the depolarization shift and, thus, a sizable frequency shift of the transition 22 . Effects of this kind have been experimentally investigated under quasi-cw illuminations 12 and theoretically analyzed for the case of pulsed illuminations 23 .
While existing approaches are sufficient to obtain accurate predictions for the effective nonlinear optical response in many specific configurations, an objective of this work is to develop a theory that is able to quantitatively capture the coherent dynamics of electrons in generic configurations including the subtle interplay between nonlinear effects and Coulomb interactions. Taking inspiration from earlier works 22,23 , we make use of a mean-field description of the electron dynamics based on a Hartree approximation of Coulomb interactions. This approach is based on a time-dependent Schrödinger-Poisson equation for the quantum electronic wavefunction in a potential that self-consistently includes the Coulomb interaction as a non-local nonlinear interaction term. Here, as a key advance, we numerically solve the time-dependent Schrödinger-Poisson equation to obtain predictions for the nonlinear optical response of the electronic system with no further approximation for all levels of excitations and electronic density.
The structure of the paper is the following. In Sec. II we introduce the physical system under investigation and we present our theoretical framework based on the Schrödinger-Poisson equation. This formalism is first applied in Sec. III to the calculation of the linear dynamics under weak illumination, recovering, among other, the depolarization shift of the ISB transition due to Coulomb interactions. We then move on in Sec.IV to study the op-tical response of a doped quantum well hosting a single intersubband transition. To show the predictive power of our technique, we carry out a microscopic investigation of the suppression of the depolarization shift under strong illuminations. Crucial further outcomes of this study are simple scaling laws connecting the nonlinear response in different wavelength regions from the Mid-IR to the THz as well as a critical discussion of the relation between the coherent nonlinearities and competing incoherent nonlinearities due to the shelving of electrons in dark states. In Sec. V, we investigate the consequences of these nonlinearities for a stack of doped QWs embedded in an optical cavity and operating in the strong lightmatter coupling regime. Besides the direct red-shift of the ISB transition due to the suppressed depolarization shift, in the microcavity case also the reduction of its oscillator strength leads to a sizable shift of the polariton branches via the reduction of the polariton Rabi splitting. These two effects are in competition on the lower polariton branch, but add up to reinforce the effective nonlinearity experienced by the upper polariton. These predictions are of great importance in the context of intersubband polariton physics with potential application to ISB polariton condensation and lasing 24,25 . As a final, more speculative topic, in Sec. VI we discuss the promise of these nonlinearities in view of translating to the Mid-and Far-IR domains the strong quantum nonlinearities in of polariton blockade originally predicted for interband exciton-polaritons in the near-infrared and visible domains [26][27][28] . Conclusions and perspectives are finally drawn in Sec. VII.

II. THE MODEL
We consider electrons in quantum well systems that are translationally invariant along the xy plane perpendicular to the growth axis z. All electrons are assumed to be initially located in the lowest subband of the QW and share the same wavefunction ψ(z) along the growth axis z. As usual, the total antisymmetry condition of the (fermionic) electronic wavefunction is ensured by the electronic motion along the xy plane, all electronic states of the lowest subband being filled up to a given Fermi energy. Future work will deal with the multi-subband case where the electron density is so large that the Fermi energy exceeds the ISB transition energy as well as with QW with parabolic energy potential where many subbands are populated 20,21,29 . Since such transitions are inherently independent of the thermal distribution of electrons, they can operate up to room temperature even at THz frequency, offering interesting applicative perspectives for this frequency range.
At the level of the Hartree description of Coulomb interactions, we neglect quantum correlations among electrons and we assume that each electron is subject to the electrostatic potential generated by the average electron density. This picture can be formalized in a time-dependent non-local nonlinear Schrödinger-Poisson equation, for the single-particle electron wavefunction φ(z), normalized to dz |ψ(z)| 2 = 1.
Here, m * is the effective mass of electrons in the solidstate material, V QW (z) is the effective confining potential of the QW, e is the electron charge, and is the background dielectric constant. Furthermore, σ imp (z) is the z-dependent three-dimensional density of the dopant impurities and the three-dimensional electronic density is given by σ el |ψ(z)| 2 in terms of the overall surface electron density σ el . The global neutrality is ensured by dz σ imp (z) = σ el . The non-locality of the nonlinear interaction term stems from the long-range nature of the Coulomb potential between electrons which, in the present planar geometry, is proportional to the distance |z − z |. In the limit of a vanishing electronic density, the nonlinear interaction term vanishes and the evolution of the single-particle wavefunction ψ reduces to a linear Schrödinger equation. Finally, the last term accounts for the time-dependent electric field E(t) of the electromagnetic wave that excites the electrons.

III. STATIC PROPERTIES AND LINEAR DYNAMICS
The first step of our calculation consists of determining the static distribution of charges in the ground state of the system. Because of the Coulomb interaction, the ground state wavefunction is in fact distorted from the lowest eigenstate of the QW potential V QW (z). In order to numerically determine the ground state configuration, we have used an imaginary-time evolution to determine the self-consistent minimum energy eigenstate of the SP equation (1). This imaginary-time evolution is a standard numerical technique used to determine the ground state of interacting many-body systems at the level of mean-field approximation, for instance dilute Bose-Einstein condensates of ultracold atoms 30 .
Examples of such calculations are illustrated in the figures. The convergence of the energy to the ground state energy is illustrated in the upper panel of Fig. 1, while snapshots of the electronic wavefunction approaching the ground state are shown as blue lines in the insets of this panel. The black lines display the total static potential felt by electrons, equal to the sum of the QW potential V QW (z) and the Coulomb interaction potential.
Examples of the final wavefunction are shown in the insets of Fig.2(a-b) for the same QW thickness and total electron density σ el = 3 · 10 12 cm −2 but different spatial energy during the imaginary-time SP evolution towards the ground state electronic wavefunction starting from a random initial wavefunction. The two curves refer to a substantial electron density (dashed) and a vanishing electronic density (solid). The insets show the spatial profile of the total potential including the QW and the Coulomb interaction potentials (black lines) and the squared modulus of the electronic wavevfunction (blue lines). Solid/dashed lines refer again to different values of the electronic density. The two insets refer to different times of the imaginary-time evolution, as indicated by the arrows. Lower panel (b): Time-dependence of the electric field amplitude of the incident pulse (red, left axis) and of the induced electronic polarization (blue, right axis). The insets show the effective potential (black) and the squared modulus of the electronic wavefunction (blue) at two different times indicated by the arrows. Solid/dashed lines refer again to different values of the electronic density.
distributions of the dopant impurities, either symmetrically located outside the well [top panel (a)] or inside in the well [bottom panel (b)]. As before, the solid black lines display the total static potential felt by electrons, while the blue lines show the squared modulus |ψ i (z)| 2 of the different eigenstates in the total potential, corre-sponding to the different subbands. The baseline of each wavefunction is vertically shifted by the energy of each electronic eigenstate. Throughout all this work, we will always assume that the electron population is initially concentrated in the lowest subband.
Even though the energy differences between the singleparticle energy levels shown in the insets of Fig.2 give a qualitative indication of the transition energies, an accurate determination of the collective mode frequencies requires the use of the full SP equation (1) since the electronic motion unavoidably leads to a time-dependent modification of the Coulomb potential. In the literature on intersubband transition, this frequency shift goes under the name of depolarization shift [18][19][20][21] .
The crucial difference between the single-particle states and the collective excitation frequencies is well known from the theory of the dilute Bose gas 31,32 . Here, the time-evolution is described at the mean-field level by a nonlinear equation very similar to (1) except for the spatially local form of the nonlinear interaction term, the so-called Gross-Pitaevskii equation (GPE). In a spatially homogeneous geometry, the single-particle states have the quadratic dispersion of the kinetic energy, while the Bogoliubov dispersion of the excitation modes obtained from the linearization of the GPE starts at low wavevector with a sonic dispersion dominated by the interaction effects and only later recovers the quadratic form. In our spatially confined geometry, an exact solution of the nonlocal Bogoliubov equations obtained by linearizing the SP equation (1) is not straightforwardly done with analytical techniques and requires numerical methods. For this reason, in this work we adopt an equivalent, fully numeric strategy that has the key advantage of being directly extended to the nonlinear regime.
The idea is to perturb the system with a suitably chosen weak electromagnetic pulse E(t). The carrier frequency of the pulse is chosen to be in the vicinity of the transition frequency of interest. The pulse duration is chosen to be long enough in time not to excite multiple excitation modes, but also short enough to give a sufficiently wide bandwidth that easily covers the desired transition. As a rule of thumb, Gaussian pulses with a duration in the several 10 fs range are typically a good choice to efficiently excite ISB transitions in the Mid-IR without having to fine-tune the carrier frequency.
The overall pulse strength is chosen to be weak enough to be in well in the linear regime and avoid (for the moment) all nonlinear effects. We then record the temporal evolution of the system at later times in response to the perturbation. As a key observable, we calculate the timedependence of the average electronic polarization An example of the temporal profile of the applied electric field pulse is shown as a red line in the bottom panel of Fig.1 together with the resulting time-dependence of z(t) (blue curve). From this latter, the transition fre-quency ω res is straightforwardly obtained by means of a sinusoidal fit of the form or, alternatively, by looking for the dominant peak in the Fourier transform of z(t). Both these strategies give equivalent results and, most importantly, are completely independent of the details of the specific initial excitation pulse used in the calculation. Note that the fact that in our calculation the oscillations persist for indefinitely long times is a consequence of having neglected all those electronic decoherence processes 33 that would lead to a decay of the dipole moment at a characteristic rate γ ISB typically on the order of a fraction of ps. Inclusion of decoherence processes into our Schrödinger-Poisson formalism will be the subject of future work. The result of a few such calculations for realistic QW parameters are shown by the left-most points of the different curves of Fig.2. These indicate the value of the oscillation frequency for different electron densities in the limit of a weak perturbation. For small electron densities (triangles), the Coulomb interactions have a negligible effect and the oscillation frequency recovers the energy separation ∆ o 21 = E 2 − E 1 of the two lowest eigenstates of the Schrödinger equation.
The significant blue-shift of the oscillation frequency that is observed for growing electron densities (rhombus, pentagons and circles) is a result of Coulomb interactions. Here, the importance of the depolarization shift 18-21 is easily appreciated comparing the numerically obtained oscillation frequency with the energy difference ∆ 12 between the lowest energy levels in the total potential shown in the insets for the highest electron density value, namely ∆ 21 93 meV for impurities located outside the well (top) and ∆ 21 111 meV for impurities located inside the well (bottom). In addition to the depolarization shift, note also the importance of static Coulomb interactions. This is visible by comparing the single-particle energy difference ∆ 21 with its value ∆ o 21 110 meV in the vanishing electron density limit: the difference is specially marked for the former configuration.

IV. NONLINEAR DYNAMICS
In the previous Section we have introduced our theoretical framework and we have given a first confirmation of its efficiency by numerically recovering the well-known depolarization shift of the oscillation frequency in the weak excitation regime. In the present Section we move on to the study of nonlinear effects. Among the many nonlinear effects that ISB may feature in different configurations, in this work we focus on the shifts of the oscillation frequency that appear under stronger excitations and on the related drop of the oscillator strength.
The numerical results will be interpreted in terms of simple analytical relations and scaling laws that allow to FIG. 2: Oscillation frequency versus deposited energy per electron for growing values of the electron densities, from a very small σ el ∼ 0 cm −2 value (red triangles), to higher ones σ el = 1 · 10 12 cm −2 (green diamonds), 2 · 10 12 cm −2 (yellow pentagons) and 3 · 10 12 cm −2 (blue dots). The two panels refer to different location of the dopant impurities outside (top) and inside (bottom) the QW. Solid lines are guides to the eyes based on a linear regression of the numerical data points. Within each panel, the inset shows (black line) the total potential experienced by the electrons and resulting from the sum of the QW potential and the Coulomb interaction energy and (blue lines) the squared modulus of the bound single-particle eigenstates in the total potential for the highest value of the electronic density considered in the main panels. Each wavefunction is shifted upwards by an amount corresponding to its eigenenergy. In all calculations of this work, we take values = 12.9 for the background dielectric constant and m * = 0.067m el for the effective mass of electrons. For the chosen QW parameters, the bare ISB transition is in the Mid-IR range at an energy around 110 meV. The vertical dashed lines indicate the position of ∆E = E0/2. quickly bridge different regimes and wavelength ranges. The experimental consequences of our results will be finally discussed in comparison with competing incoherent effects.

A. Numerical results
Nonlinear effects start being visible for larger values of the applied electric field amplitude. As a crucial step to unravel the different mechanisms underlying the nonlinear effects, we need to isolate the intrinsic nonlinear dependence of the oscillation frequency on the excitation level from the nontrivial dynamics of the excitation process during the electric field pulse. As it was pointed out in 23 , the amount of energy that gets effectively delivered to the electrons may in fact be strongly affected by nonlinear effects pushing the transition on-or off-resonance during the pulse itself.
To suppress these effects, we classify our numerically extracted value of the oscillation frequency as a function of the additional energy (per electron) ∆E that is deposited in the electronic motion at the end of the excitation sequence. As done in the linear regime, the nonlinear oscillation frequency ω res is extracted from the sinusoidal fit of the time-dependent dipole z(t) at late times or, equivalently, from its Fourier transform. The excitation energy is estimated from the total (kinetic plus potential plus interaction energy) of the electrons after the pulse has ended. Once again, provided the excitation pulse spectrum is narrow enough to only excite the desired ISB transition, we have verified that the extracted value of ω res does not depend on the details of the excitation pulse but only on the energy that gets deposited in the QW as expected on physical grounds. Fig.2 shows the oscillation frequency ω res as a function of the deposited energy ∆E for two different spatial distributions of the dopant impurities, namely outside (top) or inside (bottom) the well: as expected from the approximated calculations in Ref. 22, the oscillation frequency ω res shows an approximately linear red-shift for growing ∆E. As we have already seen for the depolarization shift, the slope of this dependence changes for different locations of the donors. For a given geometry, the different lines corresponding to different values of the electron density σ el are found to approximately cross at a single point on the ω res = E 0 = ∆ o 21 horizontal line. Together with the linear dependence of the depolarization shift on σ el , this points towards an approximated analytical form for the resonance frequency, where the adimensional a 0,1 parameters only depend on the well geometry and not on the electron density. It is interesting to note how the depolarization shift term in (4) has a similar functional form as the standard expression where ω pl = [4πe 2 σ el /(m * L QW )] 1/2 is the plasma frequency of electrons in the QW and a 0 is related to the constant geometrical coefficient 2 /(m * L 2 QW E 0 ) of order one.
For the configurations of Fig.2, the a 0,1 parameters are approximated by a 0 ∼ 0.35 and a 1 ∼ 1.5 for the top Similar plots for the amplitudez of the electron oscillations as a function of the excitation energy are shown in Fig.3 for the same well geometries and electron densities as in the previous figure. Physically,z corresponds to the effective electric dipole moment of the transition, a quantity closely related to its oscillator strength. The general features of this plot are easily understood in the case of a negligible electron density (red triangles) where one recovers the usual behaviour of two-level atoms. The squared dipole momentz 2 grows linearly for small excitation energies ∆E, then saturates to a maximum value for ∆E = E 0 /2 = 55 meV where the population is equally distributed between the two levels, and finally decreases again in the population-inversion regime.
Except for an overall reduction of the dipole moment due to the modification of the electronic wavefunctions by Coulomb interactions, a qualitatively very similar behaviour is found for higher electron densities σ el , with a linear growth ofz 2 at small ∆E and a saturation at larger ∆E. For all geometries and all values of the electron density considered, our simulations confirm the physical expectation that the transition is saturated when approximately half of the electronic population is in the excited states of the well.
This behaviour can be summarized in an approximated analytical form where the adimensional b 0,1 and s 1 coefficients only depend on the well geometry. In all considered cases, from the data displayed in the Fig.3 one can extract a value b 1 ∼ 1 for the coefficient in the nonlinear dipole moment reduction factor, which gives a factor 1/2 suppression of the squared dipole moment when ∆E/E 0 1/2 and electrons are equally (yet coherently) distributed between the two subbands. This is in stark contrast with the case of incoherent nonlinearities, where an equal population of the two states leads to a full quenching of the dipole moment 10-13 .

B. Scaling laws
The coefficients in the analytical forms (4-6) have an interesting interpretation in terms of scaling laws. For the purpose this discussion, let us restrict for simplicity to quantum wells of thickness L QW with infinitely high barriers. As compared to the more realistic configuration considered in the Figures, having infinite barriers only introduces quantitatively minor differences.
If lengths are measured in units of the well thickness ζ = z/L QW and, correspondingly, time in units of the inverse kinetic energy τ = t/(mL 2 QW ), the Schrödinger-Poisson equation can be recast in an adimensional form where the renormalized wavefunctionψ = L 1/2 QW ψ is normalized such that dζ |ψ(ζ)| 2 = 1, the renormalized impurity density isσ imp = σ imp L QW , and electric field isẼ = eEmL 3 QW / 2 . The functionṼ (ζ) describes the infinite barriers and is defined as V = 0 for |ζ| < 1/2 and +∞ otherwise. The relative strength of the Coulomb interactions is quantified by the coefficient η Coul = 2πe 2 mL 3 QW σ el / 2 on the second line.
The universality of this form provides simple scaling laws under which the coefficients in Eqs. (4)(5)(6) are invariant. This provides a straightforward way to extend our analytical results to all frequency ranges without the need of repeating the numerical calculation. For the sake of concreteness, let us consider an example of experimental interest: let us assume that the well width is changed by a factor α −1/2 so that the resonance frequency E 0 is multiplied by α. Given the typical experimental observation that the relative value γ ISB /E 0 of the linewidth of the ISB transition is approximately constant over the different wavelength regions, a useful benchmark is to keep the same value for the strength η Coul of the Coulomb interactions and thus of the relative value of the depolarization shift compared to the transition energy. To this purpose, one needs to vary the electron density by a factor α 3/2 . For instance, reducing the resonance frequency by a factor 10 from 110 meV (in the Mid-IR) to 11 meV (corresponding to ∼ 2.7 THz) requires a ∼ 3 times wider well; keeping the same η Coul then requires a ∼ 30 times lower two-dimensional electron density.

C. Competing incoherent nonlinearities
Our theoretical developments so far provide a simple, yet realistic model of those coherent optical nonlinearities that stem from the intrinsic nonlinearity of the electronic motion and the Coulomb interactions. As such, our results directly provide an estimate of the magnitude of coherent nonlinear processes. For instance, the strength of the parametric coupling underlying wave-mixing effects is directly related to the nonlinear frequency shift studied in the previous subsections.
It is however crucial to keep in mind that other nonlinear processes of incoherent nature are also typically at play for intersubband transitions in QWs. Electrons in QWs are in fact subject to different decoherence mechanisms that lead to a fast effective decay of the intersubband excitations into relatively long-lived dark electronic excitations 33 . This results in a sizable reduction in the density of active electrons participating to the electronic transition, which also leads to a quenching of the oscillator strength and, as pointed out in early works 22 , of the depolarization shift.
Leaving aside the shift of the transition frequency due to the static Coulomb interactions with the dark excitations (which is typically small in simple wells, but may become sizable in strongly asymmetric configurations 17 ), the coherent optical response can still be captured by our theory provided we identify at each time the electron density with the one of active electrons, σ el → σ act el and we estimate the evolution of this latter using a simple rate-equation model, where γ ISB is the decay rate of the coherent ISB excita-tions and γ d is the (typically much slower) decay rate of the dark excitations. Estimates for this latter are typically in the 10 ps range, much longer than the characteristic decay of the coherent ISB excitation on the order of a fraction of ps. At steady-state under a monochromatic excitation, the fraction of active electrons is reduced to Since in typical samples γ d γ ISB , the reduction in the density of active electrons can be important already at small excitation levels ∆E/E 0 1. This suggests that in quasi-CW illumination regimes the nonlinear shift of the resonances receives a dominant contribution from incoherent saturation effects. If one is interested in incoherent nonlinear processes such as bleaching and/or a frequency shift of the ISB transition, the long relaxation time γ d is a beneficial feature to reduce the required incident power 34 .
In spite of the presence of incoherent effects, the coherent nonlinearities that underlie wave-mixing effects remain however active and display an interestingly different scaling with σ el : given the form of the last, nonlinear term of (4), the excitation level ∆E/E 0 needed to obtain a given value of the coherent nonlinearity decreases as σ −1 el for growing electron density σ el . Via (9), this implies that the incoherent saturation effect can be reduced by increasing the electron density σ el . Furthermore, under the experimentally reasonable assumption that γ d does not change much when moving from the Mid-IR towards the THz range, this same equation suggests that the relative effect of the incoherent nonlinearities is reduced for longer wavelengths.
Finally, it is important to note that all these arguments hold a continuous-wave illumination of the sample by, e.g., a quantum cascade laser source 35 . A promising alternative to further suppress the incoherent effects is to use a pulsed excitation. In this regime, the experiment can be carried out on a fast enough time scale that the interesting coherent nonlinear dynamics occurs before a sizable amount of dark excitations is generated 36 .

V. OPTICAL NONLIENARITIES IN MICROCAVITIES
One of the most promising configuration to exploit the optical nonlinearities to generate useful optical processes is to embed the QWs within high-Q microcavity devices with a small enough size to enter the strong light-matter coupling regime where the elementary excitation modes have the mixed light-matter character of polaritons. This regime allows an efficient coupling of the electronic degrees of freedom to the optical fields and, at the same time, enhance the effect of the nonlinearities.
As a concrete benchmark of the strength of our ISB nonlinearities, we will consider the light intensity value that is needed to have a frequency shift comparable to the linewidth. This is the typical condition under which important nonlinear effects such as optical bistability 37 or optical parametric oscillation start occurring 25 . The discussion that follows will mostly concentrate on the latter effect, which is a promising strategy to achieve lasing and Bose-Einstein condensation effects in novel wavelength regions 24,38 . Since the amplitude of the parametric coupling between the pump and the signal/idler modes is quantitatively related to the frequency shift, it is natural to characterize the parametric oscillation threshold in terms of the ratio between the frequency shift and the decay rate.

A. Nonlinearity from the depolarization shift
Let us restrict for simplicity to the experimentally relevant case of a single-sided cavity containing N w wells at critical coupling with incident resonant radiation 39,40 ., Here, we can estimate the excitation density in each well for a given incident power P inc as where γ pol is the polariton decay rate. For simplicity, in what follows the |u X | 2 Hopfield coefficient quantifying the matter component of the polariton is assumed to be 1/2. Keeping in mind that the frequency shift of the polariton mode is |u X | 2 times the one of the matter excitation 37 , we can combine this formula with the analytical formula (4) for the ISB frequency shift to obtain an explicit expression for the power-dependent frequency-shift of the polariton, From this formula, we can estimate that a power around 1 MW/cm 2 is required to obtain a red-shift of the polariton mode comparable to the polariton linewidth γ pol = 5 meV in a cavity containing N w = 10 wells. Quite interestingly, note that this formula does not involve the electron density. This is of course valid as long as one remains in the strong coupling regime 34 . In contrast, as we have pointed out in the previous Subsection, the incoherent saturation gets suppressed for growing electron density at a fixed ∆E/E 0 . Based on the scaling laws discussed above, and plugging in the typical experimental observation that the quantity γ pol /E 0 (that is, the Q factor) is typically constant across the different frequency windows, one obtains that that the required power to achieve a red-shift comparable to γ pol scales as α 7/2 and thus quickly decreases as one moves to longer wavelengths.
As a concrete example, reducing the resonance frequency by a factor 10 from 110 meV to 11 meV (corresponding to ∼ 2.7 THz) reduces the P inc by a remarkable factor ∼ 3000 towards the few 100 W/cm 2 range. Further reductions could come from a reduction of the number N w of wells (keeping a constant a fixed overlap factor), an improvement of the cavity Q-factor, or a clever design of the cavity so to spatially concentrate the light intensity in subwavelength volumes 41,42 .

B. Nonlinearity from the saturation of polariton splitting
When dealing with microcavity configurations, it is important to remind that an additional frequency shift of the polariton modes arises from the nonlinear saturation of the dipole moment which induces a corresponding reduction of the polariton Rabi splitting 37 .
Within our theory, this effect is captured by the nonlinear dependence of the dipole moment in (6), which gives a corresponding variation of the oscillator strength, f f 0 (1−b 1 ∆E/E 0 ). In terms of the polariton splitting Ω R , this results into where the linear-regime Rabi frequency is given by Here, η is an adimensional parameter of geometric origin, typically of order one, while L cav is the thickness of the cavity.
It is interesting to quantitatively compare the magnitude of the nonlinear shift due to the ISB frequency shift (4) to the one coming from this reduction of Ω R . To this purpose, we can consider the ratio Plugging in the specific parameters for the Mid-IR QW considered above with an electron density σ el = 3 · 10 12 cm −2 , a QW density N w /L cav = 0.02 nm −1 and η = 1/2, one obtains Ω o R 16 meV. Inserting this value into (14), one finds that the ratio of the two nonlinearities is of order unity. Interestingly, the second factor on the RHS of (14) is constant under our usual scaling while, for a given overlap factor N w /L cav , the ratio Ω o R /E 0 displays a slow variation as α −1/4 . As a result, one can not expect major changes in this ratio when moving from Mid-IR towards the THz range. This confirms that our arguments on the scaling of the required incident power with operation wavelength remain valid when we also include this saturation nonlinearity.
On the other hand, a very different behaviour is expected when comparing the two polariton branches. For the upper polariton, both mechanisms result in fact in a red-shift and cooperate to reinforce the nonlinear effect. For the lower polariton, instead, they push in opposite directions and, depending on the actual value of the Hopfield coefficients of the two polaritons, they may cancel out, suppressing the final value of the effective nonlinearity experienced by the lower polariton branch. These arguments suggest that the upper polariton branch is more favourable for nonlinear optics experiments. This conclusion is to be contrasted with the exciton-polariton case where the frequency shift of the exciton under the effect of the repulsive binary interactions is in the blue direction, making the lower polariton more favourable to nonlinear experiments 37 .

VI. EFFECTIVE QUANTUM HAMILTONIAN
Even though the theory presented in this work is based on the excitation of the electronic system by classical light, our results are a good starting point to build a phenomenological quantum theory of optical nonlinearities of electrons in QWs. Such a development is of utmost importance if one is to extend quantum optics concepts, tools and applications originally developed for visible or near-IR light 43 to devices operating in longer wavelength ranges.
Indicating withΨ X (r) the (approximately bosonic) field operator describing the bright ISB excitation mode of the electrons in the QW 20,44 ,we can write a model Hamiltonian in the form: where ω lin = E 0 + a 0 e 2 L QW / is the linear oscillation frequency and the binary interaction energy ω nl = −a 1 e 2 L QW / accounts for the red-shift of the resonance. Note that this effective interaction term has an opposite sign compared to the exciton-polariton case 37 . In analogy to this exciton-polariton case, a positive interaction coefficient between ISB polaritons g > 0 was erroneously assumed by two of us for the calculations in 25 , with no major prejudice for the overall conclusions of that work. The term describing the coupling to the applied electric field E(r, t) has the physical meaning of an effective transition dipole. At linear regime its value is while at higher densities displays a saturation behaviour of saturation density n sat = σ el /b 1 . As usual, the oper-atorΨ † X (r)Ψ X (r) indicates the in-plane density of ISB excitation quanta. Within the microscopic framework of Ref. 45, our quantum Hamiltonian is applicable in the "bosonic" regime of a relatively large number of electrons and relatively weak excitation.
As a simplest example of application of this model Hamiltonian, it is interesting to estimate the strength of the single-excitation nonlinearity in a subwavelength resonator of lateral area S cav where the electromagnetic field is confined in all three-dimensions [46][47][48][49] . As a figure of merit, we will consider the frequency shift ∆ 1 of the polariton resonance when a single quantum of excitation is injected into the device. When ∆ 1 exceeds the linewidth γ pol , the presence of a single quantum of excitation is able to push the oscillation frequency away from resonance with the incident light and, in this way, prevent the injection of a second quantum of energy into the device. This phenomenon goes under the name of photon/polariton blockade and is experimentally visible as strong non-classical features in the transmitted and reflected light such as antibunching 26,37,50 .
In order to estimate ∆ 1 , we first note that the energy of a single quantum of excitation will distribute among the N el = σ el S cav electrons present in the resonator, giving an excitation density ∆E/E 0 = 1/N el . As a result, For the Mid-IR configuration considered in this work, the first fraction on the RHS is of order 0.4, so blockade ∆ 1 /γ pol ≥ 1 requires the number of electrons N el to be a sizable factor below the Q factor of the cavity, indicated here by the E 0 /( γ pol ) factor. Assuming that electrons are uniformly distributed in the cavity area with a given two-dimensional density, this imposes an upper bound on the cavity area. As a quantitative benchmark, for the electron density value σ el = 3 · 10 12 cm −2 used so far, each electron effectively occupies a region of ∼ 30 nm 2 . Using subwavelength mid-infrared nano-antennas and comparable electron densities, it was possible to achieve a lateral confinement of the field strong enough to observe strong light-matter coupling with a few 10 3 electrons confined in an area of a characteristic size of 100 nm 49 .
Based on our scaling arguments for the different coefficients in (17), it is immediate to see that the criterion based on the Q factor and the number N el of electrons directly extends to longer wavelength regimes and leads to comparable if not more promising predictions. The larger dipole moment of the transition reduces in fact the required electronic density and thus weakens the constraint on the maximum physical size of the patch cavity to observe blockade. As a result, ultra-strong light-matter coupling has been observed using a cyclotron transition of less than 100 electrons coupled to a sub-THz nanogap hybrid LC microcavities 47 and using a 3 THz ISB transition of a few 1000 electrons coupled to a LC resonator 48 .
Comparing these values with the Q factors in the 20 range that are presently available and considering the perspectives of further improvement sketched in the original works, these results are extremely promising in view of reaching polariton blockade in the Mid-IR and THz domains in the next future. Given the present state of technology, it is likely that a main experimental challenge will concern our ability of single photon detection for such long-wavelength radiation.

VII. CONCLUSIONS AND PERSPECTIVES
To summarize, in this work we have developed a general theory of the coherent optical nonlinearities associated to intersubband transitions in semiconductor quantum wells including the quantum mechanical motion of electrons and the Coulomb interactions. As most relevant observable quantities, simple expressions for the intensity-dependence of the oscillation frequency and the dipole moment of the intersubband transition are derived. Interesting scaling laws in the operation wavelength are highlighted and crucial differences from competing processes such as incoherent saturation effects are pointed out in the different regimes of excitation. The consequences of these optical nonlinearities on intersubband polaritons in microcavity geometries are investigated and quantitative estimates across different ranges of wavelengths from the Mid-IR to the THz are put forward. These predictions appear promising in view of the observation of novel phenomena such as parametric gain and Bose-Einstein condensation of intersubband polaritons. As a more speculative direction, we have finally explored the potential of intersubband polaritons as a platform for exploting blockade effects to generate antibunched light in longer wavelength ranges where quantum optics is still much less developed.
The theoretical framework discussed in this work will be of great use in future work to design structures maximizing the strength of the nonlinear response for different processes such as intensity-dependent frequency shifts, harmonic generation, and coherent wave-mixing processes. On the longer run, our results will be a useful starting point to build a fully quantum optical theory to guide experiments aiming at extending, e.g., photon blockade phenomena and single-photon emission to novel wavelength regimes.