C Barrera

V Castro

F Escudero

J J Cruz

I Verdugo

J Yon

A Fuentes
email: andres.fuentes@usm.cl

Characterization of anisole sooting properties in a laminar coflow diffusion lame

  L'archive ouverte pluridisciplinaire

Introduction

Nowadays, combustion processes of conventional fuels have become a major environmental concern due to the emission of greenhouse gases and particulate matter generated on high temperature processes, which not only contribute to global warming but also put the health of the human population at risk [START_REF] Kumfer | The role of soot in the health effects of inhaled airborne particles[END_REF]. In an attempt to reduce these emissions, oxygenated fuels such as alcohols, ethers, and non-edible oils have been suggested to be used as blends with biodiesel [START_REF] Verma | A review on performance, combustion and emissions utilizing alternative fuels[END_REF]. It is expected that as the presence of oxygenated fuels is increased, the combustion efficiency will improve and emissions will be reduced [START_REF] Kalligeros | An investigation of using biodiesel/marine diesel blends on the performance of a stationary diesel engine[END_REF]. Thus, these alternative fuels have been proposed as partial or total replacement of conventional fuels currently used. Among oxygenated alternative fuels, one that has gained increased attention is anisole (C 7 H 8 O). Anisole is an oxygenated aromatic fuel which is obtained from lignin processing. It has the potential to be used as biofuel in diesel engines through its blending with conventional fuels like n-heptane [START_REF] Xiong | Sooting tendency of n-heptane/anisole blends in laminar counterflow flames with oxygen enrichment[END_REF], which is known by its low toxicity, tendency to suppress soot formation, and contribution to reduce hydrocarbon (HC) and CO emissions [START_REF] Rajesh | Effect of anisole addition to waste cooking oil methyl ester on combustion, emission and performance characteristics of a DI diesel engine without any modifications[END_REF]. However, there are limited studies that deal with soot formation characteristics in anisole flames. For example, Zhang et. al [START_REF] Zhang | Chemical effects of anisole and toluene addition to n-heptane on PAH characteristics in laminar premixed flames by LIF measurement and kinetic model[END_REF] studied the effects of anisole and toluene addition to n-heptane. Through combination of experimental measurements and a kinetic model, they found that toluene and anisole addition resulted in notable differences on the PAH growth processes. Wang et. al [START_REF] Wang | Study of soot formations in co-flow laminar diffusion flames of n-heptane and oxygenated aromatic biofuels from atmospheric condition to 2.3 bar[END_REF] studied the characteristics of soot formation on anisole and benzyl alcohol blends with n-heptane at atmospheric and elevated pressure, observing an enhancement on soot formation processes with increasing oxygenated additives. In addition, they found a power-law increasing dependence of soot volume fraction with pressure, which was larger with the addition of anisole than pure n-heptane.

Besides paying attention to the fuel used, it is also important to take into account the oxidizing environment in which the fuel is burned. The more oxygen on the coflow oxidizer, the less nitrogen is available, reducing the energy losses by sensible heat. In addition, increased oxygen will result in higher flame temperatures, enhancing soot formation processes, but also enhancing oxidation processes both for availability of oxygen and increased temperature. Several studies have analyzed the effect of the oxygen molar fraction on the coflow oxidizer (oxygen index, OI) of a diffusion flame in terms of soot production, temperature, and thermal radiation from the flame [START_REF] Escudero | Unified behavior of soot production and radiative heat transfer in ethylene, propane and butane axisymmetric laminar diffusion flames at different oxygen indices[END_REF][START_REF] Cortés | Effect of Fuels and Oxygen Indices on the Morphology of Soot Generated in Laminar Coflow Diffusion Flames[END_REF][START_REF] Cepeda | Influence of the oxygen index on acoustically forced laminar ethylene non-premixed flames[END_REF]. However, only limited studies have focused on the effect of the oxygen index on sooting propensity of anisole [START_REF] Xiong | Sooting tendency of n-heptane/anisole blends in laminar counterflow flames with oxygen enrichment[END_REF].

In this study, we focus on characterizing experimentally the sooting propensity and maturity of an anisole laminar coflow diffusion flame, under different oxygen-enriched environments. Specifically, the flame characterization is based on non-intrusive optical diagnostics of multiwavelength line-of-sight attenuation measurements, allowing to compute soot volume fraction ( 𝑓 𝑠 ), temperature (𝑇 𝑠 ), local radiative power per unit volume (∇ • 𝑞 𝑅 ), and a degree of maturity of soot particles (through a maturity coefficient 𝛽). The purpose of this study is twofold. First, to use the retrieved flame properties under different oxygen indices to better understand the soot formation and oxidation processes of anisole, and second, to construct a comprehensive database of two-dimensional fields of these properties for validation of detailed soot production models for complex liquid fuels, such as anisole.

Experimental methodology

Experimental apparatus and procedures

Figure 1 shows a schematic of the experimental setup used for the multiwavelength line-of-sight attenuation (MW-LOSA) measurements. A modified version of the Gülder burner (6) [START_REF] Patiño | Soot primary particle sizing in a n-heptane doped methane/air laminar coflow diffusion flame by planar two-color TiRe-LII and TEM image analysis[END_REF] was used to produce the laminar coflow diffusion flames. To avoid fuel condensation, the fuel tube and outer surface are heated with silicon heater bands and insulated by a Ertalon black cover. The temperature of the heater was controlled by a temperature controller Omega CN7533 1/32 [START_REF] Rajesh | Effect of anisole addition to waste cooking oil methyl ester on combustion, emission and performance characteristics of a DI diesel engine without any modifications[END_REF]. To generate a stable coflow laminar flame, anisole is mixed with a carrier gas, and injected into the inner tube of the burner. The carrier gas is methane, which is expected to not greatly affect soot formation processes of the anisole fuel [START_REF] Zhu | Experimental study on the effect of hydrogen addition on methane/ethylene diffusion flame soot formation based on light extinction measurement[END_REF]. The liquid anisole was stored (1) and pressurized (2) using N 2 , and then delivered to a vapour delivery module (VDM, Bronkhorst SW-100 (3)) at a pressure of 3 bar. The VDM is used to vaporize the fuel at 175ºC, and a Coriolis mass flow controller (Bronkhorst F-201CV-500 [START_REF] Cortés | Effect of Fuels and Oxygen Indices on the Morphology of Soot Generated in Laminar Coflow Diffusion Flames[END_REF]) is used to precisely control the fuel mass flow. The vaporized fuel was delivered to the burner through a 175ºC heated hose [START_REF] Xiong | Sooting tendency of n-heptane/anisole blends in laminar counterflow flames with oxygen enrichment[END_REF]. Based on a previous study [START_REF] Patiño | Soot primary particle sizing in a n-heptane doped methane/air laminar coflow diffusion flame by planar two-color TiRe-LII and TEM image analysis[END_REF], the volumetric flow rate of methane is set at 257.5 ml/min, whereas the mass flow rate of anisole is set to 0.25 g/h, resulting on a constant molar fraction of anisole of 𝑋 𝑣𝑎 𝑝 = 4.7 × 10 -2 . These conditions result on a laminar flame with a height slightly below the smoke point (SP) for all conditions. An oxidizer coflow composed of OI% of oxygen and (1 -OI)% of nitrogen was injected into the coannular section of the burner, flowing at a constant volumetric flow rate of 140 slpm and controlled by MFC Brooks SLA5800 mass flow controller [START_REF] Escudero | Unified behavior of soot production and radiative heat transfer in ethylene, propane and butane axisymmetric laminar diffusion flames at different oxygen indices[END_REF]. Variations on the OI were accomplished by controlling the ambient air flow (assumed to be composed of 21% oxygen and 79% nitrogen) and doping it with pure oxygen, according to Table 1. A fast photodiode [START_REF] Cepeda | Influence of the oxygen index on acoustically forced laminar ethylene non-premixed flames[END_REF] coupled to a DAQ [START_REF] Patiño | Soot primary particle sizing in a n-heptane doped methane/air laminar coflow diffusion flame by planar two-color TiRe-LII and TEM image analysis[END_REF] was focused on the flame tip to check the flame stability through a fast Fourier transform (FFT) analysis. Briefly, MW-LOSA consists on measuring the line-of-sight integrated flame transmissivity (𝜏 𝜆 ) and radiative intensity (𝑃 𝜆 ) at different wavelengths and then, through deconvolution processes and data treatment, retrieve the local soot properties of interest. To accomplish these measurements, a modulated back-light source [START_REF] Zhu | Experimental study on the effect of hydrogen addition on methane/ethylene diffusion flame soot formation based on light extinction measurement[END_REF] was used to generate the projected images, following the procedure described by Yon et al. [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF]. A LED source [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF] coupled to an integrating sphere ( 14) is used to establish a diffusive light source, which is collimated by an achromatic lens of 70 mm diameter and focal length f = 150 mm [START_REF] Eaves | CoFlame: A refined and validated numerical algorithm for modeling sooting laminar coflow diffusion flames[END_REF]. This beam, after traversing the flame, passes through another achromatic lens f = 400 mm which focuses the transmitted light on the sensor of an EMCCD camera [START_REF]Non-grey gas radiative transfer analyses using the statistical narrow-band model[END_REF] coupled to a filter wheel [START_REF] Thurgood | The TN Quadrature Set for the Discrete Ordinates Method[END_REF]. Flame images with and without back-light (from the LED) were captured at each wavelength, setting the camera acquisition rate at twice the frequency of the back-light [START_REF] Mohammad | The coreshell internal nanostructure of soot-A criterion to model soot maturity[END_REF]. During the procedure, 250 images were acquired by the EMCCD camera for the different OI conditions, with 4 different narrow-band detection wavelengths (𝜆 = 500, 532, 660 and 810±5 nm), using different exposure times for each wavelength (0.125, 0.07, 0.013 and 0.013 s, respectively) in order to enhance signal-to-noise ratio. Finally, the EMCCD was calibrated by absolute radiance using a broadband NIST radiometric source (Labsphere USLR-V) The EMCCD camera response was recorded for each filter used, obtaining calibration factors (𝛾 = 15.5 × 10 3 , 31.2 × 10 3 , 14.6 × 10 3 , 35.5 × 10 3 W/m 3 /sr) from visible to infrared, respectively.

Soot temperature

The soot temperature 𝑇 𝑠 distribution was estimated using the 𝜏 𝜆 and 𝑃 𝜆 measurements from the MW-LOSA campaign. The wavelength 𝜆 = 810 nm was selected to reduce the scattering effects, which may not be negligible for wavelengths lower than 800 nm [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF]. To retrieve 𝑇 𝑠 , we need information about the local absorption coefficient (𝑘 𝑎,𝜆 ) and flame emission 𝐽 𝜆 = 𝐼 𝑏𝑏 𝜆 𝑘 𝑎,𝜆 , where 𝐼 𝑏𝑏 𝜆 is the black-body intensity. An onion-peeling deconvolution process with Tikhonov regularization was used to retrieve 𝑘 𝑎,𝜆 and 𝐽 𝜆 at each height (𝑧) and radial position (𝑟) of the flame. With these properties and using Planck's law, soot temperature can be determined as:

𝑇 𝑠 = - 𝑘 𝑏 𝜆 ℎ𝑐 ln 𝜆 5 𝐽 𝜆 2ℎ𝑐 2 𝛾𝑘 𝑎,𝜆 -1 (1)
where 𝑘 𝑏 is the Boltzmann constant, ℎ the Planck constant, 𝑐 the speed of light in vacuum, and 𝛾 = 35.5 × 10 3 W/m 3 /sr is the calibration constant determined as described previously. Details of the procedure can be found elsewhere [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF].

Soot maturity

Soot maturity determination is based on the different spectral dependence of the absorption function 𝐸 (𝑚) (𝑚 is the wavelength-dependent refractive index), according to the presence of organic and graphitic components on a soot particle [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF][START_REF] Bescond | Soot optical properties determined by analyzing extinction spectra in the visible near-UV: Toward an optical speciation according to constituents and structure[END_REF]. From the 𝐽 𝜆 and 𝑇 𝑠 fields, the experimental ratios of 𝐸 (𝑚, 𝜆 𝑖 ) were obtained at each pair of detection wavelengths (𝜆 𝑖 , 𝜆 𝑗 ):

𝐸 (𝑚, 𝜆 𝑖 ) 𝐸 (𝑚, 𝜆 𝑗 ) = 𝐽 𝜆 𝑖 𝜆 𝑖 𝐼 𝑏𝑏 𝜆 𝑗 (𝑇 𝑠 ) 𝐽 𝜆 𝑗 𝜆 𝑗 𝐼 𝑏𝑏 𝜆 𝑖 (𝑇 𝑠 ) . (2)
The spectral dependence of the ratio of 𝐸 (𝑚) can be expressed as a power-law function of the ratio of wavelengths, controlled by an exponent 𝛽, termed maturity coefficient [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF]:

𝐸 (𝑚, 𝜆) 𝐸 (𝑚, 𝜆 𝑟𝑒 𝑓 ) = 𝜆 𝑟𝑒 𝑓 𝜆 𝛽 , (3)
where 𝜆 𝑟𝑒 𝑓 is a reference wavelength and 𝐸 (𝑚, 𝜆 𝑟𝑒 𝑓 ) the corresponding value for the absorption function. Note that 𝛽 = 𝜒 -1, where 𝜒 is the Ångström absorption exponent. 𝛽 can be obtained as the best fit parameter which minimizes deviations of the experimental ratios of Eq. ( 2) and the theoretical ones from Eq. (3). From the relation described by the theoretical ratio (3), 𝛽 = 0 means a non-spectral dependency (mature soot) whereas 𝛽 > 0 indicates the presence of certain amount of organic compounds on the soot particle. Through the mixing law employed by Bescond et al. [START_REF] Bescond | Soot optical properties determined by analyzing extinction spectra in the visible near-UV: Toward an optical speciation according to constituents and structure[END_REF], a relationship between 𝐸 (𝑚, 𝜆), 𝛽 and the amount of organic compounds of soot can be obtained, considering absence of amorphous soot, i.e., only organic and graphitic soot. Finally, 𝐸 (𝑚, 𝜆) can be obtained from this relationship and Eq. ( 3).

Soot volume fraction

Assuming that soot particles follow the Rayleigh-Debye-Gans approximation and the absorption cross section is proportional to their volume, the soot volume fraction ( 𝑓 𝑠 ) can be obtained as 𝑓 𝑠 = 𝜆𝑘 𝑎,𝜆 /6𝜋𝐸 (𝑚, 𝜆). Soot volume fraction was obtained from the measurements of 𝑘 𝑎,𝜆 and considering the spatial distribution of soot absorption function 𝐸 (𝑚, 𝜆) determined at 𝜆 = 810 nm, as detailed by Yon et al. [START_REF] Yon | Revealing soot maturity based on multi-wavelength absorption/emission measurements in laminar axisymmetric coflow ethylene diffusion flames[END_REF].

Soot radiation

Using the experimental fields of 𝑓 𝑠 and 𝑇 𝑠 , we can solve the radiative transfer equation (RTE) in cylindrical coordinates to estimate the local flame radiation (∇ • 𝑞 𝑟 ), and total radiation emitted ( 𝑄 𝑅 ) by soot particles within the flame, as it is performed in numerical codes [START_REF] Eaves | CoFlame: A refined and validated numerical algorithm for modeling sooting laminar coflow diffusion flames[END_REF]. To solve the RTE, the discrete ordinates method (DOM) is used to discretize the whole solid angle (𝜔), using a 𝑇 3 quadrature [START_REF] Thurgood | The TN Quadrature Set for the Discrete Ordinates Method[END_REF][START_REF]Non-grey gas radiative transfer analyses using the statistical narrow-band model[END_REF]. The solution of the RTE provide the directional radiative intensities 𝐼 𝜆 . Then, ∇ • 𝑞 𝑟 is computed by integrating the contributions for all wavelengths and directions as:

∇ • 𝑞 𝑅 = ∫ ∞ 0 4𝜋𝑘 𝑎,𝜆 𝐼 𝑏𝑏 𝜆 -𝑘 𝑎,𝜆 ∫ 4𝜋 𝐼 𝜆 𝑑𝜔 𝑑𝜆 (4)
The total radiation emitted by the anisole flame is computed integrating Eq. 4 over the complete flame volume:

𝑄 𝑟 = 2𝜋 ∬ ∞ 0 ∇ • 𝑞 𝑟 𝑟𝑑𝑟𝑑𝑧. ( 5
)
Error propagation A Montecarlo error-propagation analysis was performed to quantify the uncertainties of each soot property. The measured fields of 𝜏 𝜆 and 𝑃 𝜆 were perturbed by adding a Gaussian noise of zero mean and standard deviation derived from experimental measurements using the 250 images. Then the whole procedure for retrieving each property was performed for 75 generated samples. The statistical analysis of these provide values for the mean and standard deviation of each soot property, which are the error bars provided on the subsequent plots. Further information can be found on the Supplementary Material.

Results and discussions

Fields distribution

Figure 2 shows the experimentally retrieved fields for soot temperature, volume fraction, radiation and maturity coefficient (from top to bottom) for increasing oxygen enrichment conditions (from left to right). The overall features of the flame when increasing the oxygen content are similar to those encountered on typical hydrocarbons [START_REF] Escudero | Unified behavior of soot production and radiative heat transfer in ethylene, propane and butane axisymmetric laminar diffusion flames at different oxygen indices[END_REF]. Figure 2(a) shows the temperature fields computed using Eq. 1. Overall, the flame temperature increases as the OI is increased, specially near the flame edges. The mean temperature is reported on the top of each field, showing monotonic increase from 1638 K to 1830 K as the OI is incremented from 21% to 35%. The increase on the soot temperature results on an increase of soot volume fraction, as shown on Fig. 2(b). The total soot loading, defined as the volume integral of 𝑓 𝑠 , Γ = 2𝜋 ∬ ∞ 0 𝑓 𝑠 𝑟𝑑𝑟𝑑𝑧, is indicated at the top of each field. The 𝑓 𝑠 distribution is spread through the flame zone at OI = 21%, meaning that both the center and wing zones of the flame present similar soot loadings. As the OI is increased, 𝑓 𝑠 is higher on the whole, and 𝑓 𝑠 on the wings is promoted more than 𝑓 𝑠 on the center zone of the flame. This suggests that the hydrogen-abstraction carbonaddition (HACA) mechanism, that controls the soot loading on the flame wings, is enhanced in a more pronounced fashion as compared to the PAH condensation growth mechanism, that dominates the flame center. These effects on 𝑇 𝑠 and 𝑓 𝑠 also impact the luminous flame height. There is a monotonic decrease of the luminous flame height from approximately 53 mm for OI = 21% to 26 mm for OI = 35%. As 𝑇 𝑠 increases, soot formation processes are enhanced, resulting in higher 𝑓 𝑠 , but at the same time the oxidation processes are also enhanced due to both 𝑇 𝑠 and higher O 2 content, resulting on a shorter flame [START_REF] Escudero | Unified behavior of soot production and radiative heat transfer in ethylene, propane and butane axisymmetric laminar diffusion flames at different oxygen indices[END_REF][START_REF] Cortés | Effect of Fuels and Oxygen Indices on the Morphology of Soot Generated in Laminar Coflow Diffusion Flames[END_REF]. It is interesting to note that, although the local soot content is increased for higher OIs, the total soot loading Γ monotonically decreases, indicating that the decrease of flame volume dominates over the increase on local 𝑓 𝑠 . Fig. 2(c) presents the divergence of the radiative heat flux, which combines the effects of 𝑇 𝑠 and 𝑓 𝑠 . 𝑄 𝑅 is indicated at the top. As the radiative calculations are performed from the infrared (IR) to the far-infrared (FIR) region of the spectrum, the wavelength dependent absorption coefficient is assumed to be a power-law function of the ratio of wavelengths, controlled by maturity coefficient 𝛽, in the same way as for the absorption function. The reference wavelength and absorption coefficient are taken at 𝜆 = 810 nm, being the most IR measurement available. These assumptions may play a role on the absolute value of ∇ • 𝑞 𝑅 , but the qualitative tendencies should be well-captured. As 𝑓 𝑠 and 𝑇 𝑠 are increased, the local values of ∇ • 𝑞 𝑅 are also increased as more oxygen is added to the coflow. Now it is interesting to see that the total radiation emitted by the flame (analogous to the total soot loading) presents an increase until 31%, and then starts to decrease again, indicating that the increase of total radiation due to the fourth-power of 𝑇 𝑠 and linear with 𝑓 𝑠 through the absorption coefficient, dominate over the decrease in flame volume until 31%, where these effects start to compensate.

Finally, Fig. 2(d) shows the field of maturity coefficient 𝛽. In general, for all conditions there is some spectral dependence of the ratio of absorption function (𝛽 ≠ 0), indicating the presence of organic compounds within the soot particle. As the OI is increased, temperature is increased and 𝛽 is decreased, evidencing that higher temperatures promoted the graphitization of soot particles. This observation is in line with the view that 𝑇 𝑠 should play a predominant role on the maturation of soot particles [START_REF] Mohammad | The coreshell internal nanostructure of soot-A criterion to model soot maturity[END_REF][START_REF] Hope A Michelsen | A review of terminology used to describe soot formation and evolution under combustion and pyrolytic conditions[END_REF].

Local analysis

In this section, a more detailed analysis of the experimental fields is developed. First, it is appropriate to define locations of interest to be analyzed. In general, centerline values of soot properties, values at the location of maximum 𝑓 𝑠 (wings), radial profiles and integrated quantities are of interest for validation of numerical models. A quantity that is useful to analyze the global amount of soot and its production rates at each height above the burner (HAB) is the radially integrated soot volume fraction, defined as 𝐵 𝑠 = 2𝜋 ∫ ∞ 0 𝑓 𝑠 𝑟𝑑𝑟. Its derivative, 𝑑𝐵 𝑠 /𝑑𝑧 can be used as a measure of soot production rates, where 𝑑𝐵 𝑠 /𝑑𝑧 > 0 means that formation processes dominates over oxidation, and vice-versa. In addition, a mean formation/oxidation ( ᾱ 𝑓 𝑜𝑟/𝑜𝑥 ) rate is defined as the mean value of 𝑑𝐵 𝑠 /𝑑𝑧 before/after 𝐵 𝑠,𝑚𝑎𝑥 . The location of maximum 𝐵 𝑠 , or 𝑑𝐵 𝑠 /𝑑𝑧 = 0, establishes the HAB where formation and oxidation processes are balanced. This point will be taken as reference HAB for the subsequent analysis. Figure 3 shows the results for 𝐵 𝑠 for each HAB, as continuous lines at the left axis, and 𝑑𝐵 𝑠 /𝑑𝑧 as dashed lines at the right axis. Fig. 3(a) shows these values as a function of the HAB, and Fig. 3(b) shows the maximum 𝐵 𝑠 and mean formation and oxidation rates. As the OI increases, the plots of 𝐵 𝑠 are shifted to lower HABs and become narrower due to the decrease in flame height. All 𝐵 𝑠 curves finish on 𝐵 𝑠 = 0, confirming that all flames are non-smoking. The production rates become more intense, similar to other studies [START_REF] Escudero | Unified behavior of soot production and radiative heat transfer in ethylene, propane and butane axisymmetric laminar diffusion flames at different oxygen indices[END_REF][START_REF] Cepeda | Impact of water-vapor addition to oxidizer on the thermal radiation characteristics of non-premixed laminar coflow ethylene flames under oxygendeficient conditions[END_REF], the oxidation rates being more sensitive to the OI than the formation rates. However, the maximum values of 𝐵 𝑠 remain in a range between 177-215 ppm mm 2 , indicating that the localized increase in 𝑓 𝑠 in a narrow zone near the flame wings is compensated by the decrease on flame radius. Fig. 4 shows the plots of the soot properties as a function of the HAB for two zones of interest, namely, the position of maximum soot volume fraction at each height (or flame wing, first column) and the flame centerline (second column). As seen on Fig. 4(a), temperature increases as the HAB is increased. There is a sharp decrease in 𝑇 𝑠 for each curve at the 𝑓 𝑠,𝑚𝑎𝑥 location approximately at middle of the flame for all conditions, that simply denotes the point where the 𝑓 𝑠,𝑚𝑎𝑥 shifts from the wing to the flame centerline. In general, flame temperature is higher at the flame wings near the reaction zone, and at the flame tip. As the presence of oxygen in the flame increases, the maximum temperature gradually increases, reaching a peak temperature of 2234 K for an OI = 35%. Fig. 4(b) shows that 𝑓 𝑠 is promoted on the flame wings as the OI is increased, whereas 𝑓 𝑠 on the centerline remains approximately constant, between 4-6 ppm for all conditions. The HAB where 𝐵 𝑠 is maximum is almost coincident with the maximum 𝑓 𝑠 , except for an OI of 21%, where 𝑓 𝑠,𝑚𝑎𝑥 continues to increase until approximately an HAB of 40 mm. The increase in temperature and soot volume fraction due to the enriched OI conditions, results on an increase of ∇ • 𝑞 𝑅 , as shown in Fig. 4(c), specially at lower HABs on the flame wings where the production of soot is more promoted. The maximum values of ∇ • 𝑞 𝑟 almost coincide with the maximum 𝐵 𝑠 , near the zone where oxidation starts to dominate over formation processes. Finally, the maturity coefficient 𝛽 is lower at the location of maximum 𝑓 𝑠 than at the flame centerline, showing that more mature particles are found near the hightemperature reaction zone of the flame near the wings. Figure 5 shows the radial profiles of the flame properties at the height of maximum 𝐵 𝑠 . Temperatures increase with increasing presence of oxygen on the coflow, and from the flame center to the wing, in a similar fashion for the different OI conditions. The non-smoothness of soot volume fraction comes from the consideration of spatial variation of 𝐸 (𝑚), that in turn comes from the results of 𝛽. The global features already discussed are also present on this plot, where the soot volume fraction at the wings is affected in a more pronounced way than the centerline as the OI is increased, and this behavior is repeated for ∇ • 𝑞 𝑅 . Finally, the maturity coefficient 𝛽 follows a decreasing tendency from the center of the flame wing, the variation being more prominent as OI is increased. As the OI is increased, there is higher occurrence of elevated 𝑓 𝑠 , but the low value of its PDF denotes that these values are localized, as seen on the fields distributions. There is a clear shift to the right (higher temperature) for the PDFs of 𝑇 𝑠 as the OI is increased, and a similar behavior is observed for the local radiation. Finally, it is interesting to note that the PDF of the maturity coefficient seems to correlate to the 𝑇 𝑠 ; higher temperatures results on lower 𝛽. Figure 6(e) shows the plot 𝑇 𝑠 with their corresponding 𝛽, for each position (𝑟, 𝑧) within the flame. We can observe that as the temperature increases, there is an overall decrease of the maturity coefficient 𝛽, meaning that on the whole, particles are more graphitized at higher OIs. At an OI of 35%, the plot is more spread, meaning that for the same temperature, there is a wide range of maturity levels of the particles that are exposed to that temperature. This may be due to the effect of amorphous components that are not taken into consideration in this study, or the time-temperature histories of soot particles. Further studies with other fuels and coupled with numerical simulations are needed to decouple these chemical effects, and to provide more insights about the maturation and soot production processes on complex fuels such as anisole.

Conclusions

An experimental study consisting on multi-wavelength line-of-sight attenuation measurements was performed to provide further insights of soot production and maturation processes in anisole fueled flames at different oxygen-enriched conditions. Two dimensional field distributions of soot temperature, volume fraction (taking into account the spatial variation of soot absorption function), divergence of the radiative heat flux, and maturity coefficient, were computed from flame extinction and emission measurements and analyzed in a global and local fashion. In general, as the oxygen index (OI) is increased, soot formation and oxidation processes are enhanced, resulting on higher local soot loadings, and flames that are shorter and more luminous. These results are inline with observations on typical hydrocarbons at different OIs. In addition, the maturity coefficient 𝛽, was found to decrease on the whole as the OI is increased, which suggested that there is a correlation between high temperatures and lower 𝛽 (higher amount of graphitic contents). Through the analysis of the PDFs of the different properties, and visualizing the correlation between 𝛽 and temperature, it was found that higher temperatures produce more graphitic soot in the studied anisole flames. For the same temperature, there is a range of maturity levels that the particles can achieve, which may depend on the time-temperature history of the soot particle.

Figure 1 :

 1 Figure 1: Experimental setup used for multi-wavelength line-of-sight attenuation measurements in the anisole flame.

Figure 2 :

 2 Figure 2: (a) Spatially-resolved soot volume fraction, (b) temperature, (c) divergence of the radiative heat flux computed with DOM and (d) 𝛽 maturity index for different oxygen index conditions.

Figure 3 :

 3 Figure 3: Radially integrated soot volume fraction 𝐵 𝑠 (continuous lines, left axis) and mean formation and oxidation rates through the slope of the 𝐵 𝑠 curve (dashed lines, right axis) as a function of the HAB.

Figure 4 :

 4 Figure 4: (a) Soot volume fraction, (b) temperature, (c) divergence of the radiative heat flux and (d) maturity coefficient at the of the maximum 𝑓 𝑠 (left column) and the flame centerline (right column) as function of the HAB. As insets graphs, variation of each property with the OI at the HAB of maximum 𝐵 𝑠 is presented.

Figure 5 :

 5 Figure 5: Radial profiles of (a) 𝑓 𝑠 , (b) 𝑇 𝑠 , (c) ∇ • 𝑞 𝑟 and (d) 𝛽 on the 𝐵 𝑠,𝑚𝑎𝑥 HAB location

Figure 6 :

 6 Figure 6: (a)-(d) Probability density function generated from each property of the soot over the entire flame. (e) Plot of 𝛽 vs 𝑇 𝑠 for every point of the flame, with a linear fit of all the data and its equation.

Figure 6 (

 6 Figure6(a)-(d) shows the probability density function (PDF) of soot properties, representing the whole information of each flame. As the OI is increased, there is higher occurrence of elevated 𝑓 𝑠 , but the low value of its PDF denotes that these values are localized, as seen on the fields distributions. There is a clear shift to the right (higher temperature) for the PDFs of 𝑇 𝑠 as the OI is increased, and a similar behavior is observed for the local radiation. Finally, it is interesting to note that the PDF of the maturity coefficient seems to correlate to the 𝑇 𝑠 ; higher temperatures results on lower 𝛽. Figure6(e) shows the plot 𝑇 𝑠 with their corresponding 𝛽, for each position (𝑟, 𝑧) within the flame. We can observe that as the temperature increases, there is an overall decrease of the maturity coefficient 𝛽, meaning that on the whole, particles are more graphitized at higher OIs. At an OI of 35%, the plot is more spread, meaning that for the same temperature, there is a wide range of maturity levels of the particles that are exposed to that temperature. This may be due to the effect of amorphous components that are not taken into consideration in this study, or the time-temperature histories of soot particles. Further studies with other fuels and coupled with numerical simulations are needed to decouple these chemical

Figure S2 :

 S2 Figure S2: Relative errors for (a) soot volume fraction, (b) temperature, (c) divergence of the radiative heat flux, (d) 𝛽 maturity index and (e) absorption function at 810 nm, for different oxygen index conditions.

Table 1 :

 1 Experimental conditions for the oxidizer flow

	OI (%)	21	25	31	35
	𝑄 𝑎𝑖𝑟 (slpm) 140 132.91 122.28 115.19
	𝑄 𝑂 2 (slpm)	0	7.09	17.72 24.81

Acknowledgments

This work was partially funded by Chile's ANID projects FONDECYT/Regular XXXX, FONDE-CYT/Postdoctoral 3210498 and 3190860, new PIA/ANILLO ACTXXXX?.

Supplementary material Soot absorption function

Figure S1 shows the fields of absorption function obtained at 810 nm.

Uncertainty analysis

The fields of relative errors for 𝑇 𝑠 , 𝑓 𝑠 , ∇ • 𝑞 𝑅 , 𝛽, and 𝐸 (𝑚, 810), defined as the values of the standard deviation (denoted as 𝛿) divided by the mean value (denoted by an overbar), is shown in Fig. S2. These error were computed using 75 perturbed fields using the mean and standard deviation from the experimental measurements of 𝜏 𝜆 and 𝑃 𝜆 .


R Cominotti

H A M Leymann

J Nespolo

J.-M Manceau

M Jeannin

R Colombelli

I Carusotto

Theory of coherent optical nonlinearities of intersubband transitions in semiconductor quantum wells

Keywords:

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ödinger-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 Bose-Einstein condensation of intersubband polaritons and to the speculative exploration of quantum optical phenomena such as single-photon emission in the Mid-IR and THz windows.

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 in-teractions 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, the objective of this work is to develop a theory that is able to quantitatively capture the coherent dynamics of electrons in generic configurations and include 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 yields 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 for arbitrary levels of excitation strength 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 optical response under weak illumination, recovering, among other, the depolarization shift of the ISB transition due to Coulomb interactions at high electronic densities.

In Sec.IV we move on to study the nonlinear optical response to stronger excitations. In particular we point out a marked intensity-dependence of the ISB resonance and of the effective dipole moment. The numerical results are interpreted in terms of a nonlinear quenching of the oscillator strength and the consequent suppression of the depolarization shift. Analytical scaling laws connecting the nonlinear response in different wavelength regions from the Mid-IR to the THz are then extracted from the numerics. A critical discussion of the relation between the coherent nonlinearities and competing incoherent nonlinearities due to the shelving of electrons in dark states is finally provided.

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 light-matter coupling regime. Besides the direct red-shift due to the suppressed depolarization shift, in the microcavity case a sizable shift of the polariton branches also occurs from the reduction of the oscillator strength at high intensities. While these two effects are in competition on the lower polariton branch, they add up constructively to reinforce the upper polariton nonlinearity. 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 ISB nonlinearities in view of translating to the Midand Far-IR domains those polariton blockade effects that were originally predicted for interband exciton-polaritons in the near-infrared and visible range [26][27][28] , so to explore quantum nonlinear optics phenomena in a novel spectral window. Conclusions and perspectives are finally drawn in Sec. VII.

II. THE MODEL AND THE STATIC PROPERTIES

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. 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 extension of our formalism to the multi-subband case 20,21,29 where the electron density is so large that the Fermi energy exceeds the ISB transition energy and the coupled evolution of several wavefunctions corresponding to the different subbands has to be simultaneously determined.

At the level of the Hartree description of Coulomb interactions, we neglect quantum correlations among elec-trons 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 timedependent non-local nonlinear Schrödinger-Poisson (SP) equation,

iℏ ∂ψ(z, t) ∂t = - ℏ 2 2m * ∂ 2 ψ ∂z 2 + V QW (z)ψ(z) + 2π e 2 ϵ dz ′ |z -z ′ | σ imp (z ′ ) -σ el |ψ(z ′ )| 2 ψ(z) -e E(t) z ψ(z) (1)
for the single-particle electron wavefunction ψ(z), normalized to dz |ψ(z)| 2 = 1.

V QW (z) is the effective confining potential of the QW and, for simplicity, we have assumed a parabolic electronic band with a z-independent effective mass m * at all points. Generalization to non-parabolic and/or spacedependent mass case can be carried out with straightforward modifications. 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.

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) and can be obtained as the self-consistent minimum energy eigenstate ψ g (z) of the nonlinear SP equation,

E g ψ g (z) = - ℏ 2 2m * ∂ 2 ψ g ∂z 2 + V QW (z)ψ g (z)+ + 2π e 2 ϵ dz ′ |z -z ′ | σ imp (z ′ ) -σ el |ψ g (z ′ )| 2 ψ g (z) . (2)
In order to numerically solve this eigenvalue equation for ψ g (z) and the single-particle energy E g in specific configurations, an imaginary-time evolution technique has been adopted. This is a standard numerical technique used to determine the ground state of interacting manybody systems at the level of mean-field approximation, for instance dilute Bose-Einstein condensates of ultracold atoms 30 . A brief account of the principles of this numerical method can be found in Appendix A. Examples of such calculations are illustrated in Fig. 1. The convergence of the energy to the ground state energy is illustrated in Fig. 1(a) for two different values of the electron density (solid and dashed lines). Snapshots of the electronic wavefunction approaching the ground state are shown as blue lines in the insets of this panel. Here, the black lines display the total potential V tot (z) felt by electrons, equal to the sum of the QW potential V QW (z) and the Coulomb interaction potential.

Plots of the final wavefunction are shown in Fig. 1(b) for the same QW thickness and total electron density σ el = 3 • 10 12 cm -2 but different spatial distributions of the dopant impurities, either symmetrically located outside the well (modulation doping, left panel) or inside in the well (right panel), indicated by the red arrows. As before, the solid black lines display the total potential V tot (z) felt by electrons, while the blue lines show the squared modulus |ψ i (z)| 2 of the different eigenstates in the total potential, corresponding to the different subbands. The baseline of each wavefunction is vertically shifted by the energy of each electronic eigenstate. For the sake of simplicity, throughout all this work, we will always assume that the electron population is initially concentrated in the lowest subband. The energy differences between the single-particle energy levels shown in the Fig. 1(b) give a qualitative indication of the transition energies, but 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] .

An exact prediction for the collective oscillation frequencies in the linear regime of weak excitations could be obtained by linearizing the SP equation (1) around the ground state ψ g determined in the previous Section. This gives the linearized equation of motion

iℏ ∂ δψ(z, t) ∂t = - ℏ 2 2m * ∂ 2 δψ(z) ∂z 2 + V QW (z) δψ(z) + 2π e 2 ϵ dz ′ |z -z ′ | σ imp (z ′ ) -σ el |ψ g (z ′ )| 2 δψ(z)+ -σ el δψ * (z ′ ) ψ g (z ′ ) + ψ * g (z ′ ) δψ(z ′ ) ψ g (z) (3)
for the perturbation δψ(z) defined as

ψ(z) = e -iEgt/ℏ [ψ g (z) + δψ(z, t)].
Differently from the Bogoliubov equations for atomic gases with contact interactions [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF][START_REF] Castin | Coherent atomic matter waves[END_REF] , here the linearized equation features a nonlocal Coulomb interaction term. An analytical solution for the eigenmodes of (3) can be obtained only in suitable limits as discussed in Appendix B. In the general case, it typically requires numerical methods to diagonalize the linear problem ? . For this reason, in this work we adopt an equivalent, yet fully numeric strategy that has the key advantage of being directly extended to the nonlinear regime. As it is sketched in Fig. 2(a), the idea is to simulate the real-time evolution starting from the ground state wavefunction ψ g (z) obtained by the imaginary-time evolution discussed in the previous Section and look at the small oscillations induced by 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 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 time-dependence of the average electronic polarization

z(t) = dz z |ψ(z, t)| 2 . ( 4
)
An example of the temporal profile of the applied electric field pulse is shown as a red line in Fig. 2 together with the resulting time-dependence of z(t) (blue curve). From this latter, the transition frequency ω res and the oscillation amplitude z are straightforwardly obtained by means of a sinusoidal fit of the form ? z(t) ≃ z cos(ω res t + φ) .

(5)

Note that the fact that in our calculation the oscillations persist for indefinitely long times is a consequence of having neglected all the electronic decoherence processes [START_REF] Unuma | [END_REF] 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 in Fig. 2(b), where the oscillation frequency ω res in response to a very weak perturbation in the linear regime is plotted as a function of the electron density for dopant impurities located outside the well (modulation doping) or inside the well. For small electron densities, the Coulomb interactions have a negligible effect and the oscillation frequency recovers the energy separation E o 21 of the two lowest eigenstates of the linear Schrödinger equation in the bare potential V QW (z) of the well.

The significant blue-shift of the oscillation frequency that is observed for growing electron densities is a result of Coulomb interactions. Here, the importance of the depolarization shift [18][19][20][21] is easily appreciated comparing the numerically obtained oscillation frequency (solid and dashed lines) with the static Hartree prediction given by the energy difference E 21 between the two lowest energy single-particle orbitals in the total static potential V tot including interactions (square symbols) ? . Additional analytical insight on the physical origin of this frequency difference is presented in Appendix B.

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 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 done in the previous Section for the linear regime, the oscillation frequency ω res and the oscillation amplitude z are extracted from the sinusoidal fit of the time-dependent dipole z(t) at late times. 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 and z do 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.

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 Ref. 23, the amount of energy that gets effectively delivered to the electrons may in fact be strongly affected by nonlinear effects that push the transition frequency 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. Once the pulse has gone, energy is conserved, so ∆E can be estimated by numerically calculating the increase in the SP energy

E SP = dz ℏ 2 2m * ∂ψ ∂z 2 + V QW (z) |ψ(z)| 2 + + π e 2 ϵ dz dz ′ |z -z ′ | σ imp (z ′ ) -σ el |ψ(z ′ )| 2 |ψ(z)| 2 (6)
between the wavefunction ψ at a generic late-time t late and the one at the initial time t = 0, ∆E = E SP (t late ) -E SP (t = 0). An intuitive understanding of the physical meaning of the ∆E parameter can be obtained by relating it to experimentally accessible quantities. For sufficiently weak Coulomb interactions, we can neglect nonlinearities and approximately write ∆E = p 2 E o 21 , where p 1,2 are the fractional occupations of the ground and excited subband, p 1 + p 2 = 1. The point ∆E = E o 21 /2 then corresponds to the saturation condition p 2 = p 1 = 1/2. Whereas these arguments are exact in the low-density limit, they remain qualitatively valid also in the presence of sizable nonlinearities as long as Coulomb energy is a relatively small correction to the bare transition energy E o 21 . Additional remarks on the relation between ∆E and the incident intensity in a microcavity geometry will be given at the beginning of Sec.V 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 second term in (7) has a similar functional form as the standard expression for the depolarization shift

ℏω res E o 21 = E o 21 2 + (ℏω pl ) 2 E o 21 ≃ 1 + ℏ 2 ω 2 pl 2E o 21 2 (8)
where ω pl = [4πe 2 σ el /(m * ϵL QW )] 1/2 is the plasma frequency of QW electrons and the a 0 coefficient in ( 7) is related to the adimensional geometrical parameter ℏ 2 /(m * L 2 QW E o 21 ) of order one. For the configurations of Fig. 3, the a 0,1 parameters are approximated by a 0 ∼ 0.35 and a 1 ∼ 1.5 for the top panel and by a 0 ∼ 0.6 and a 1 ∼ 1.3 for the bottom one. As it is discussed in detail in Appendix B, the smaller value of the weak-excitation frequency shift parameter a 0 in the upper panel is mostly the effect of the static Coulomb interactions with the doping impurities, which are more significant when these latter are located outside the well and counterbalance the effect of the dynamical distortion of the electronic distribution induced by the excitation. On the other hand, the slope a 1 of the excitation-dependent shift has comparable values in the two cases as it mostly depends on depolarization shift effects. Quite interestingly, for the experimentally realistic values of the electron density considered in the figure, the nonlinear red shift can be comparable to the typical linewidths of ISB transitions (around a few meV's) already for deposited energies well below saturation of the two-level ISB transition, that is ∆E ≪ E o 21 /2 ≃ 55 meV (vertical dashed line in the Figure ). As we are going to see in the next Subsec.IV C, this statement may need revisiting if a sizable fraction of the electrons gets accumulated in the QW as incoherent dark excitations.

Oscillation amplitude

Similar plots for the amplitude z of the electron oscillations as a function of the excitation energy are shown in Fig. 4 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). In this case, one recovers the usual behaviour of two-level atoms. The squared dipole moment z2 grows linearly for small excitation energies ∆E, then saturates to a maximum value for ∆E = E o 21 /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 of z2 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 gets 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

z L QW 2 ≃ b 2 0 1 -b 1 ∆E E o 21 1 -s 1 e 2 σ el L QW ϵE o 21 ∆E E o 21 , (9)
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. 4 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 o 21 ≃ 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][11][12][13] .

B. Scaling laws

The coefficients in the analytical forms (7) and ( 9) have an interesting interpretation in terms of scaling laws. For the purpose of 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

i ∂ ψ ∂τ = - 1 2 ∂ 2 ψ ∂ζ 2 + Ṽ (ζ) ψ -Ẽ(τ ) ζ ψ+ + η Coul dζ |ζ -ζ ′ | [σ imp (ζ ′ ) -σ el | ψ(ζ ′ )| 2 ] , (10)
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.( 7) and ( 9) are invariant. This provides a straightforward way to extend our analytical results to all frequency ranges without the need of repeating the numerical calculation. Specifically, let us consider that the QW width is varied by a factor α -1/2 so that the resonance frequency E o 21 is multiplied by α. Under this change, the Coulomb interaction parameter η Coul (i.e. the relative value of the depolarization shift) stays constant if the electron density is varied 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. In the next Section, we will see how this scaling impacts the value of light intensity that is needed to observe nonlinear effects.

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 nonlinear frequency shift of the ISB resonance or, equivalently, the strength of the parametric coupling in wave-mixing processes.

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 [START_REF] Unuma | [END_REF] . This results in a sizable reduction in the density of active electrons participating to the electronic transition and thus in 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 σ act el in time using a simple rate-equation model,

σact el = γ d (σ el -σ act el ) -γ ISB σ act el ∆E E o 21 ( 11
)
where γ ISB is the decay rate of the coherent ISB excitations (of density approximately given by σ act el ∆E/E o 21 ≃ p 2 σ act el ) 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

σ act el σ el = 1 1 + ∆E E o 21 γ ISB γ d . ( 12
)
Since in typical samples γ d ≪ γ ISB , the reduction in the density of active electrons can be important already at small excitation levels ∆E/E o 21 ≪ 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 ( 7), the excitation level ∆E/E o 21 needed to obtain a given value of the (coherent) nonlinear shift decreases as σ -1 el for growing electron density σ el for fixed QW geometry. Via (12), 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 NONLINEARITIES IN MICROCAVITIES

One of the most promising configurations to exploit the optical nonlinearities discussed in the previous Section to observe useful optical processes is to embed the QWs within high-Q microcavity devices so to enter the so-called strong light-matter coupling regime. In this regime, the elementary excitation modes have the mixed light-matter character of polaritons, which allows for an efficient coupling of the electronic degrees of freedom to the optical fields and therefore enhances the effect of the nonlinearities. Building a complete theory of the nonlinear dynamics of such microcavity devices is a task that goes far beyond this work for which preliminary steps have been reported in Ref. 25,34,36. As such, the goal of this Section is to obtain a quantitative estimate of the actual strength of the ISB nonlinearities in a configuration that is most promising for applications.

As a specific benchmark quantity, we will consider the light intensity value that is needed to have a frequency shift of the polariton mode 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.

Even though our theory is fully general and can be applied to generic devices, for the sake of concreteness we keep in mind the specific example of double-metal microcavities, a most promising work-horse for studies of ISB polaritons in both the THz and mid-infrared frequency ranges 39 . In these devices, the QWs are sandwiched between two metallic layers. The back metallic layer is left unpatterned and acts as a perfectly reflecting plane, while the front one is periodically patterned to allow optical access from the far field. Owing to this single-sided geometry, all spectroscopic information can be obtained from reflectivity measurements since there is no transmitted beam. Furthermore, the efficiency of the nonlinear process can be optimized by independently tailoring the different decay channels 39,40 so to reach the so-called critical coupling regime with external radiation, where radiative and non-radiative losses are equal and on resonance all incident light is funneled into the cavity 41 .

Under this condition, the energy density stored in the cavity is simply related to the incident power P inc by ε st γ pol = P inc where γ pol is the polariton decay rate. In the strong-coupling regime, the energy of a polariton mode is shared by its light and matter component in proportion to the Hopfield coefficients, so the excitation density in each well and per electron is given by

σ el ∆E = |u X | 2 ε st N w = P inc |u X | 2 N w γ pol (13)
where N w is the number of QWs coupled to the cavity mode and |u X | 2 is Hopfield coefficient quantifying the matter component of the polariton.

A. Nonlinearity from the depolarization shift

Deep in the strong coupling, the frequency shift of the polariton mode is |u X | 2 times the one of the matter excitation 37 . This result can be combined with the analytical formula (7) for the ISB frequency shift to obtain an explicit expression for the power-dependent frequency-shift of the polariton,

∆(ℏω res ) = -a 1 e 2 L QW ϵ |u X | 4 N w γ pol E o 21 P inc . (14)
From this formula, assuming for simplicity |u X | 2 = 1/2 and a typical number N w = 10 of wells, we can estimate that an intensity around 1 MW/cm 2 is required to obtain a red-shift of the polariton mode comparable to the polariton linewidth ℏγ pol = 5 meV. 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 Sec.IV C, the incoherent saturation effect at a given value of the nonlinear frequency shift is smaller for a large electron density.

Based on the scaling laws discussed above, and plugging in the typical experimental observation that the quantity ℏγ pol /E o 21 (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 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 42,43 .

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 (9), which gives a corresponding variation of the oscillator strength,

f ≃ f 0 (1 -b 1 ∆E/E o 21 )
. In terms of the polariton splitting Ω R , this results into

Ω R = Ω o R + ∆Ω R ≃ Ω o R 1 -b 1 ∆E E o 21 . ( 15
)
where the linear-regime Rabi frequency is given by

Ω o R = 2πe 2 N w σ el η ϵ m * L cav 1/2 . ( 16
)
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 (7) to the one coming from this reduction of Ω R . To this purpose, we can consider the ratio

∆Ω R ∆(ℏω res ) = b 1 a 1 Ω R E o 21 ϵE o 21 e 2 σ el L QW . ( 17
)
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 (17), one finds that the ratio of the two nonlinearities is in the order of unity. Interestingly, the last factor on the RHS of ( 17) is constant under our usual scaling while, for a given overlap factor N w /L cav , the ratio Ω o R /E o 21 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 include this saturation nonlinearity.

Note that a very different behaviour is expected for the two polariton branches. For the upper polariton, both mechanisms give rise to 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, they may cancel out, suppressing the final value of the effective nonlinearity. These arguments suggest that in the ISB case 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 a more favourable choice for nonlinear optics 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 attempt 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 [START_REF] Walls | Quantum optics[END_REF] to devices operating in longer wavelength ranges. Building a complete theory of quantum nonlinearities is a task that goes way beyond this work, so we will restrict here to some semi-quantitative reasonings that offer an intuitive feeling of the strength of the effect.

Indicating with ΨX (r) the (approximately bosonic) field operator describing the bright ISB excitation mode of the electrons in the QW 20,[START_REF] Ciuti | [END_REF] ,we can write a model Hamiltonian in the form:

H = ℏω lin d 2 r Ψ † X (r) ΨX (r)+ + ℏω nl 2 d 2 r Ψ † X (r) Ψ † X (r) ΨX (r) ΨX (r)+ -e d 0 d 2 r E(r, t) Ψ † X (r) 1 - Ψ † X (r) ΨX (r) n sat + -h.c. ( 18
)
where ℏω lin = E o 21 + 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 .

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

d 0 = b 0 L QW 1 -s 1 e 2 σ el L QW ϵE o 21 √ σ el , (19)
while at higher densities displays a saturation behaviour of saturation density n sat = σ el /b 1 . As usual, the operator Ψ † X (r) ΨX (r) indicates the in-plane density of ISB excitation quanta. In the language of Ref. 46, our quantum Hamiltonian (18) refers to 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 [47][48][49][50] . 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,51 .

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 o 21 = 1/N el . As a result,

∆ 1 γ pol ≃ a 1 e 2 L QW σ el ϵ E o 21 E o 21 ℏγ pol 1 N el . (20)
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 o 21 /(ℏγ 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 sub-wavelength 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 in the mid-infrared with a few 10 3 electrons confined in an area of a characteristic linear size of 100 nm 50 .

Based on our scaling arguments for the different coefficients in (20), 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 for THz radiation. The larger dipole moment of the transition reduces in fact the required electronic density for strong coupling 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 48 and using a 3 THz ISB transition of a few 1000 electrons coupled to a LC resonator 49,52 .

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 hurdle along this path will consist of the development of efficient single photon detectors to measure quantum correlations 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 their 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. 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. Finally, as a more speculative direction, we have 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 with more complex potentials, so to maximize 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. with

L =    E o 21 + ∆ H + ∆ x ∆ x -∆ x -E o 21 -∆ H -∆ x    . (B3)
Here, E o 21 = E o 2 -E o 1 is the energy difference between bare electronic levels. The terms accounting for the Coulomb interactions are proportional to the dimensional η = 2πe 2 σ el /ϵ coefficient quantifying the effective strength of Coulomb interactions. In detail,

∆ H = V 21 -V 11 (B4)
is the static Hartree shift of the transition under the effect of the Coulomb interactions, expressed in terms of the static Coulomb shift of the ψ 1,2 states in the charge distribution determined by the ground state electrons and the impurities,

V 11 = -η dz dz ′ |z -z ′ | |ψ o 1 (z)| 2 -σ imp (z) |ψ o 1 (z ′ )| 2 (B5) V 21 = -η dz dz ′ |z -z ′ | |ψ o 1 (z)| 2 -σ imp (z) |ψ o 2 (z ′ )| 2 (B6)
and ∆ x accounts for the dynamical distortion of the electronic distribution induced by the excitation,

∆ x = -η dz dz ′ |z -z ′ | ψ o 1 (z) ψ o 2 (z) ψ o 1 (z ′ ) ψ o 2 (z ′ ) .
(B7) For low electron densities, all these quantities are small compared to the bare transition energy E o 21 , so the resonance energy ℏE res can be analytically calculated within perturbation theory as

ℏω res ≃ E o 21 + ∆ H + ∆ x - ∆ 2 x 2E o 21 . (B8)
These quantities can be easily computed in the idealized case of an infinite well of thickness L w for which the wavefunctions and the energies of the two lowest states and their energies have the following analytical forms,

ψ o 1 (z) = 2 L w cos πz L w , ℏE o 1 = ℏ 2 2m * π L 2 (B9) ψ o 2 (z) = 2 L w sin 2πz L w , E o 2 = ℏ 2 2m * 2π L 2 . (B10)
In particular, the value of the integrals in (B5-B7) can be straightforwardly estimated for different locations of the impurities as summarized in the upper part of Table I.

Inserting these values in the matrix (B3), one obtains the prediction for the shift of the ISB transition listed in the bottom part of Table I It is interesting to remind that this crucial difference between the single-particle energies and the collective excitation frequencies also occurs in dilute Bose gas with local interactions [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF][START_REF] Castin | Coherent atomic matter waves[END_REF] . For instance, in a spatially homogeneous geometry, the single-particle states have the quadratic dispersion of the kinetic energy, while the Bogoliubov dispersion of the collective 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.

FIG. 1 :

 1 FIG. 1: (a) Example of imaginary-time SP evolution of the energy starting from a simply guessed wavefunction (a Gaussian centered on the well) and converging towards the ground state. 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 wavefunction (blue lines) at different imaginary times. The top inset refers to an early imaginary time indicated by the arrow, well before the convergence of the algorithm. The bottom inset refers instead to a late imaginary time after convergence. Within each panel, the solid/dashed lines refer to different values of the electron density, respectively negligible (solid) and substantial (dashed). (b): Results of the calculation for two different locations of the impurities, as indicated by the red arrows. The black solid line shows the total potential experienced by the electrons. The blue lines show the squared modulus of several single-particle eigenstates for an electron density of 3 • 1012 cm -2 . Each wavefunction is shifted upwards by an amount corresponding to its eigenenergy. In all calculations within this work, we take values ϵ = 12.9 for the background dielectric constant and m * = 0.067m el for the effective mass of electrons.

  FIG. 2: (a): Example of the time-dependence of the electric field amplitude of the incident pulse (red) and of the induced electronic polarization (green). The insets show the total potential Vtot(z) (black) and the squared modulus of the electronic wavefunction (blue and orange) at two different times indicated by the stars in the main graph. (b) Numerically evaluated ISB transition frequency in the linear regime of weak excitations, as a function of the electron density (solid and dashed lines), compared with the prediction E21 of the static Hartree approximation (squares) discussed in the text and in the Appendix B. The bare transition frequency E o 21 is indicated by the horizontal black dotted line. The dopant impurities are located either outside the well (modulation doping, dashed line and open symbols) or inside the well (solid line and filled symbols.

FIG. 3 :o 21 -

 321 FIG.3: Oscillation frequency ℏωres as a function of the deposited energy per electron ∆E. Different curves refer to growing values of the electron density, 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 (a,b) panels refer to different location of the dopant impurities outside (a) and inside (b) 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 Vtot(z) in the ground state and (blue lines) the squared modulus of several 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. 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 the saturation point ∆E = E o 21 /2.

FIG. 4 :

 4 FIG. 4: Normalized squared amplitude of the oscillating electric dipole z as a function of the deposited energy ∆E for the same configurations shown in Fig.3. The vertical dashed lines indicate the position of the saturation point ∆E = E o 21 /2.

  The different columns refer to different locations of the impurities, namely outside the well [as in the left panel of Fig.1(b)], inside the well [as in the right panel of Fig.1(b)], at the center of the well. The energy shifts on the last two rows are evaluated for an electronic density σ = 3 • 10 12 cm -2 corresponding to the rightmost points in Fig.2(b).

in

  Fig.3(b). Note how the resonance shift strongly departs from the static Hartree prediction ∆ H based on the static Coulomb energy shift of single-particle orbitals.

TABLE I :

 I . In spite of the severe approximations done, this prediction is in quite good quantitative agreement with the full numerical calculation shown Table of the Coulomb interaction integrals (B5-B7) and of the energy shifts (B4) and (B8) for an infinite well of thickness Lw = 12.7 nm for which the transition energy ∆eg=109.6 meV is comparable to the case considered in the main text.

		outside ±10 nm	inside ±2 nm	center
	V11/η	0.79	-0.011	-0.06
	V21/η	0.73	-0.013	-0.02
	∆H /η	-0.06	-0.002	0.04
	∆x/η	0.11	0.11	0.11
	ℏωres -E21 9.5 meV	26.6 meV	37.7 meV
	∆H	-17.2 meV	-0.55 meV	10.3 meV

S. Barbieri, M. Ravaro, P. Gellie, G. Santarelli, C. Manquest, C. Sirtori, S. Khanna, E. Linfield, and G. Davies, Nat. Photon. 5, 306 (2011).

4 M. Bagheri, C. Frez, L. Sterczewski, I. Gruidin, M. Fradet, I. Vurgaftman, C. Canedy, W. Bewley, C. Merritt, C. S.

Acknowledgments

We acknowledge financial support from the European Union FET-Open grant MIR-BOSE (737017). The Trento team acknowledges support from the Provincia Autonoma di Trento. Continuous stimulating discussions with J. Hawecker, S. Dhillon, M. Knorr, C. Lange, R. Huber, A. Tredicucci, G. Scalari, J. Faist are warmly acknowledged.

Appendix A: Imaginary-time evolution

In this Appendix, we briefly review the main principles underlying the imaginary-time evolution method used to find the wavefunction ψ g (z) and the energy E g of the lowest-energy eigenstate of the time-independent Schrödinger-Poisson equation (2).

The partial differential equation encoding the imaginary-time evolution is obtained from the real-time evolution equation (1) by rotating the time variable t → -iβ in the complex plane,

The factor involving the norm of the wavefunction

is required at the denominator of the interaction term since the imaginary-time evolution (in contrast to the real-time one) does not conserve the norm. At long times β → ∞, the imaginary-time evolution typically converges to an exponentially decreasing wavefunction

from which one extracts the ground state wavefunction ψ g (z) = ψ ∞ (z)/||ψ ∞ || and the ground state energy E g = E ∞ . Inserting the ansatz (A3) into (1), one indeed recovers the time-independent SP equation of the form (2). In a practical calculation, we can choose a generic wavefunction as the initial state ψ(z, β = 0). We then have to numerically evolve ψ(z, β) in β according to (A1) until we reach convergence. This is determined by looking at the convergence of the SP energy of ψ(z, β) to a constant value. The imaginary-time evolution (as well as the following real-time one) is carried out using a split-step method: the evolution at each time-step is Trotter-split into the non-commuting kinetic and potential energy parts and each of them is sequentially implemented in the space in which it is diagonal, namely k-space for the kinetic energy and real-space for the potential and interaction energy terms. At each time-step, inter-conversion between the k-and the real-space and back is performed by Fast Fourier Transform.

The imaginary-time method is most transparent in the non-interacting limit where the nonlinear term in the evolution equation is negligible. In this case, we can decompose the initial wavefunction

on the orthonormal basis of eigenfunctions ψ o j=1,2,... (z) of the Schrödinger problem in the bare quantum well potential, with energy E o j=1,2,... . The random initial condition reflects into a random choice of the initial value a j (β = 0) of the expansion coefficients. By linearity, the imaginary-time evolution acts independently on each of them, a j (β) = a j (0) e -βE o j . At late times β → ∞, only the lowest-energy eigenvector survives (the higher ones j ≥ 1 are exponentially suppressed at least as e -β(E o 2 -E o 1 ) ) and ψ(z, β) converges to the lowest-energy eigenvector,

recovering the limiting form (A3).

From this discussion, it is immediate to see that the correct ground state is found independently on the choice of the initial wavefunction provided a 1 (β = 0) ̸ = 0, a condition which is satisfied by any randomly-chosen initial wavefunction. This independence from the initial condition, mathematically proven in the non-interacting case, has been numerically verified to also hold in the interacting case by repeating the calculation for different choices of the initial wavefunction.

Appendix B: Analytical study of the oscillation frequency in the limit of low electron density

As a further verification of the numerical calculations and a quick guiding tool for the design of new structures, it is interesting to look at the linearized equation (3) in the small electron density limit where an analytical treatment is possible. The discussion in this Appendix is inspired from the Bogoliubov theory of the weakly interacting Bose gas [START_REF] Pitaevskii | Bose-Einstein Condensation and Superfluidity[END_REF] .

We indicate with ψ o 1,2 (z) the two lowest states of the quantum well for negligible electron density. In the limit of a weak excitation, we can expand the linearized dynamics in the basis of these two states only,

where α(t) is the excitation amplitude and [u 2 , v 2 ] T is the projection of the linearized eigenmode [u(z), v(z)] T on the excited ψ 2 state to which we are restricting our attention, with the usual normalization

Plugging the ansatz (B1) into the linearized SP dynamics (3) and imposing a harmonic evolution of the excitation amplitude, α(t) = ᾱ e -iωrest , we get to the eigenvalue equation,

56 Note that a positive interaction coefficient between ISB polaritons g > 0 was assumed by two of us in Ref. 25. This erroneous choice was motivated by the analogy with the exciton-polariton case and, mutatis mutandis does not bring major prejudice to the overall conclusions of that work.