Electro-Optical Modulation of the Nonlinear Optical Response in a GaAs/AlGaAs Symmetric Multiple Quantum Well System

1. Introduction

The control of light propagation through semiconductor heterostructures with multiple quantum wells (MQWs) has become a subject of intense research over the past decade [1,2,3,4,5,6,7,8,9]. In MQW systems, the effective coupling between wells usually arises when the structural parameters of the quantum wells are carefully modified and cause discrete localized energy levels of each well to begin to interact strongly. This interaction leads to extended states that spread over several wells rather than isolated quantized levels [10]. The above can lead to changes in the band structure that significantly alter the interaction of light with matter, leading to notable non-linear phenomena reported in experiments, such as high-harmonic generation processes including non-linear optical rectification (NOR) [11], second harmonic generation (SHG) [12], and third harmonic generation (THG) [13], which play a crucial role in modulating the emission spectrum of incident light [14]. In contrast to bulk materials, electronic confinement in MQW systems can be tailored, leading to enhanced non-linearities [15,16]. Furthermore, in MQW systems, the possibility of engineering diverse configurations enables tunable nonlinear responses with potential applications in electro-optical modulators [17,18], far-infrared detectors [19], semiconductor optical amplifiers [20], and nonlinear processes such as four-wave mixing [21].

In recent years, two approaches for modifying the nonlinear optical (NLO) properties in semiconductor heterostructures have been studied: the non-resonant intense laser (nIL) field and the electric field [22,23,24,25,26,27]. The nIL field alters quantum confinement, whereas the electric field breaks the inversion symmetry, enabling precise control over the NLO behavior. Although modifying structural parameters–such as well and barrier widths–in MQW systems has successfully generated enhanced optical nonlinearities, practical applications often require fixed quantum well dimensions. To address this, external fields offer a compelling alternative: they enable tunable control of higher-order harmonic generation while preserving the system’s structural configuration. For example, in a combined theoretical and experimental study, Miller et al. [28] report the effects of electric fields—both parallel and perpendicular to the quantum well layers—on optical absorption in GaAs/AlGaAs quantum well structures. Their results demonstrate that an electric field, whether applied parallel or perpendicular to the quantum-well layers, induces significant changes in optical absorption. Kasapoglu et al. [29] study the linear and nonlinear optical properties in a confining potential modeled from the symmetric and asymmetric harmonic-Gaussian potential of double quantum wells under the effect of an applied magnetic field. Their results show that the structural parameters allow control of the coupling between the two wells of the system, and the asymmetry parameter induces changes in the system’s selection rules. In addition, introducing a magnetic field produces blue or red shifts in the optical properties of the system. Additionally, Ozturk et al. [30] have investigated the effects of a laser field and an electric field, applied in different directions, on the intersubband optical absorption in a graded ${\mathrm{Ga}}_{1-x}{\mathrm{Al}}_x$As/GaAs quantum-well structure. Their results demonstrate that the intersubband transitions depend directly on the applied external fields. Subsequent studies have exploited various approaches to modify the NLO properties, including temperature and pressure effects [31,32,33,34], magnetic fields [35], and excitonic effects [36].

It is well known that the application of an electric field on a MQW system modifies the confinement potential experienced by the electrons in the system, which can lead to the breaking of the inversion symmetry and, consequently, to an adjustment of the intersubband transition (ISBT) energies and dipole moment matrix elements (DMMEs), which are fundamental factors in the non-linear optical response of the system [37,38,39,40,41]. For example, the experimental work by Park et al. [42] presented a tunable nonlinear response resulting from the electrical modulation of an intersubband polariton metasurface composed of MQW. Other experimental studies have reported an electro-optical modulation effect of light at low voltages and speeds on the order of MHz in lithium niobate metasurfaces [43], electro-optical modulation in ring resonators [44,45], and significant second-harmonic generation as a result of a linear electro-optical effect in atomically thin 2D materials [46].

However, the nIL field induces a dressing effect on the confinement profile of the MQW system, modifying the effective potential by reducing the depth of the wells and decreasing the height of the barriers. Consequently, particle confinement is altered, leading to a shift in energy levels without changing the original symmetry of the system [47]. Sakiroglu et al. [48] demonstrated the effect of the nIL field on the non-linear optical properties of a Morse quantum well. Their results reveal shifts in the positions of the peaks in both the optical absorption coefficient and the refractive index. Furthermore, other studies indicate that the combined effects of nIL and magnetic fields can modulate the interaction forces between particles and alter the energy level [49]. Many previous studies on the NLO response in multiple quantum well systems have focused on its electrical or optical modulation. However, electro-optical modulation has received limited attention, although it is known that effective coupling between wells in MQW structures can induce strong light-matter interactions, significantly enhancing the nonlinear response.

Due to the external fields playing a determining role in tailoring the system’s NLO response, achieving independent and simultaneous tuning of both fields is highly desirable. For this purpose, it is convenient to propose a GaAs/AlGaAs semiconductor heterostructure, because of the mature fabrication techniques that allow precise control over the composition and dimensions of the wells [50]. In addition, these heterostructures feature well-defined energy levels and a bandgap that can be tuned by external fields [51].

In this contribution, we perform flexible control over the ISBT energies and DMMEs by applying the nIL and the electric fields. To this end, we theoretically investigate the effect of external fields on the electronic properties of the MQW system through optical and electrical means [Figure 1a]. We focus on a ${\mathrm{GaAs/Al}}_{0.3}{\mathrm{Ga}}_{0.7}$As MQW system, considering four electron subbands. Our results demonstrate that combining these external fields enables a unified electro-optical manipulation of the NLO response. Furthermore, we show that the external fields have a measurable impact on the ISBT energies and DMMEs. Finally, we present optimal configurations of the external field values that modulate the transition-energy crossings, and the manipulation of the position of the resonant peaks involved in the susceptibility coefficient, which determines the light conversion efficiency and allows dynamic and adjustable control of the NLO response.

This paper is organized as follows. In Section 2, we present the derivation of the laser-dressed potential energy for a MQW system under the influence of the nIL field and derive the expressions for the coefficients of the NLO response and the DMMEs. Section 3 defines the parameters used in our study. Section 4 presents our main results on the influence and control of the NLO response through external fields. In Section 5, we summarize our main findings.

2. Theoretical Framework

First, we outline the dressing effect of the nILF field on the confining potential. Next, we detail the derivation of the NLO response using the nonlinear susceptibility coefficients. In addition, we provide the relevant expressions for the DMMEs. We consider a symmetric MQWs system with GaAs (${\mathrm{Al}}_{0.3}{\mathrm{Ga}}_{0.7}$As) well (barrier) regions of $L_w$ and $L_b$ widths, respectively. The total length (L) of the system is the sum of 6 barriers and 5 wells. We will calculate the corresponding subband energy levels, density probability, and NLO response, such as the NOR, SHG, and THG, derived from the intersubband transitions of the MQWs under a nIL field and an electric field.

2.1. Laser Dressed Potential Energy

The time-dependent Schrödinger equation in the effective-mass approximation, under the effect of an applied external electric field and a non-resonant intense laser field, the confinement potential V changes according to the Kramers-Henneberger transformation as follows [52]:

$$i\hslash {\displaystyle \frac{\partial }{\partial t}}\psi (\overrightarrow{r},t)=\left[-{\displaystyle \frac{{\hslash }^2}{2}}\nabla ·\left({\displaystyle \frac{1}{m^{\ast }\left(\overrightarrow{r}\right)}}\nabla \right)+V(\overrightarrow{r},\overrightarrow{\alpha }\left(t\right))-e\overrightarrow{F}·\overrightarrow{r}\right]\psi (\overrightarrow{r},t).$$

(1)

where $m^{\ast }\left(\overrightarrow{r}\right)$ is the position-dependent electron effective mass, e is the absolute value of electron charge, $\overrightarrow{F}$ is the electric field, and $V(\overrightarrow{r},\overrightarrow{\alpha }\left(t\right))$ is the laser-dressed potential of the system. The quantum wells are grown along the x-axis and the external electric field lies along the growth direction, and we assume that the nIL field is polarized in the same direction. The reference frame origin is set at the geometrical center of the heterostructure.

We use the separation of variables method to obtain a 1D differential equation for the x-coordinate. In this approach, the wave function is written in terms of the wave vector ${\overrightarrow{k}}_{\perp }$ and the electron coordinate $\overrightarrow{\rho }$ in the $yz$ plane (perpendicular to the growth direction of the heterostructure) $\psi \left(\overrightarrow{r}\right)=e^{i{\overrightarrow{k}}_{\perp }·\overrightarrow{\rho }}{\mathrm{\Psi }}_n\left(x\right)$. Assuming the bottom of all energy subbands, then ${\overrightarrow{k}}_{\perp }=0$. Therefore, by applying a Fourier expansion to the potential and using Floquet’s theory [53], while considering only a high-frequency laser field, we obtain the time-independent Schrödinger equation as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill [-{\displaystyle \frac{{\hslash }^2}{2}}{\displaystyle \frac{d}{dx}}\left({\displaystyle \frac{1}{m^{\ast }\left(x\right)}}{\displaystyle \frac{d}{dx}}\right)+V(x,{\alpha }_0)& -eF(x-L/2)]{\mathrm{\Psi }}_n\left(x\right)=E_n{\mathrm{\Psi }}_n\left(x\right),\hfill \end{array}\end{array}$$

(2)

where $E_n$ is the electron energy eigenvalue associated with the eigenfunction ${\mathrm{\Psi }}_n\left(x\right)$ of the n-th subband and F is the electric field strength. A shift of $-L/2$ is introduced in the x-coordinate to set the zero of the potential energy at that point of the structure. The second term of Equation (2) accounts for the laser-dressed potential and is written as:

$$V(x,{\alpha }_0)={\displaystyle \frac{{\omega }^{\prime }}{2\pi }}{\int }_0^{2\pi /{\omega }^{\prime }}V\left[x+{\alpha }_{0}sin\left({\omega }^{\prime }t\right)\right]dt.$$

(3)

In the above equation, we have replaced the explicit form of the term associated with a classical displacement of the electron under the influence of an electric field $\overrightarrow{\alpha }\left(t\right)={\alpha }_{0}sin\left({\omega }^{\prime }t\right)\widehat{x}$. Where ${\omega }^{\prime }$ is the non-resonant frequency, $\widehat{x}$ indicate the propagation direction of the field and ${\alpha }_0\equiv e{A}_0/\left(m^{\ast }\left(x\right){\omega }^{\prime 2}\right)$, being $A_0$ the field strength and ${\alpha }_0$ the laser-dressing parameter [47].

The intensity of the nIL field is determined by the laser-dressing parameter ${\alpha }_0$, which can be defined in terms of the time-averaged intensity $I={\displaystyle \frac{1}{2}}{\left|A_0\right|}^2{ϵ}_{0}c$ of the laser field as [20]:

$$I={\displaystyle \frac{{\alpha }_0^2{\omega }^{\prime 4}{m}_{eff}^{\ast 2}{ϵ}_{0}c}{2e^2}},$$

(4)

where c is light speed and ${ϵ}_0$ is the permittivity of the vacuum. Following a similar procedure to the one described in Ref [54]: (i) We choose a specific amplitude in the high-intensity regime of the laser field (i.e. the laser-dressing parameter) ${\alpha }_0=4$ nm. (ii) For our GaAs/AlGaAs heterostructure, the characteristic intensity at high frequencies is defined by $I_c=m_{eff}^{\ast 2}{\alpha }_B^{\ast 2}{\omega }^{4}c{ϵ}_0{ϵ}_r^{1/2}/\left(2e^2\right)$ [55], where ${\alpha }_B^{\ast }={\hslash }^2{ϵ}_r/\left(m_{eff}^{\ast }k{e}^2\right)=1.4$ nm is the effective Bohr radius, $\omega /2\pi =100$ THz is the laser frequency, and ${ϵ}_r=10.89$ is the relative dielectric constant at high frequencies and k is the Coulomb constant. (iii) To deal with the problem of non-uniform effective masses, we use the average effective mass approximation [56]:

$${\displaystyle \frac{1}{m_{eff}^{\ast }}}={\displaystyle \frac{L_w}{L}}·{\displaystyle \frac{1}{m_w^{\ast }}}+{\displaystyle \frac{L_b}{L}}·{\displaystyle \frac{1}{m_b^{\ast }}},$$

(5)

$m_w^{\ast }$ and $m_b^{\ast }$ represent the effective mass corresponding to GaAs wells and AlGaAs barriers. Thus, we obtain an approximate value of the characteristic intensity at high frequencies, which is $I_c=7.3\times {10}^{11}$ W/${\mathrm{cm}}^2$. (iv) Using Equation (4) for ${\alpha }_0=4$ nm, we calculate a maximum intensity value, $I_{max}=1.8\times {10}^{12}{\mathrm{W}/\mathrm{cm}}^2$. Although the maximum value of the intensity exceeds the characteristic intensity at high frequencies, usually intensities of about ${10}^{12}{\mathrm{W}/\mathrm{cm}}^2$ are accepted in experimental procedures [57,58,59].

Finally, we solve the differential equation corresponding to the Schrödinger equation [Equation (3)] using the finite element method (FEM) implemented in the COMSOL-Multiphysics semiconductor module [60]. After determining the energy levels and corresponding wave functions, we computed the NLO properties and the relative shifts in the ISBT energies and DMMEs.

2.2. Nonlinear Optical Response

To study the NLO response, we employ a semi-classical approach. This method combines a quantum description of the interacting particle with a classical representation of the radiation field. The radiation field is modeled as monochromatic optical radiation, characterized by the electric field $E\left(t\right)$, which is polarized along the growth direction of the quantum wells, i.e.

$$E\left(t\right)=\tilde{E}{e}^{i\omega t}+{\tilde{E}}^{\ast }{e}^{-i\omega t},$$

(6)

where $\tilde{E}$ is the complex amplitude. Using the density matrix formalism and assuming a phenomenological approach to the dissipative process, the Liouville–von Neumann equation explicitly describes its time evolution via the following master equation [61]:

$${\displaystyle \frac{\partial {\rho }_{ij}}{\partial t}}={\displaystyle \frac{1}{i\hslash }}{\left[{\widehat{H}}_0+{\widehat{H}}_I,\widehat{\rho }\right]}_{ij}-{\mathrm{\Gamma }}_{ij}{\left(\widehat{\rho }-{\widehat{\rho }}^{\left(0\right)}\right)}_{ij},$$

(7)

where $\widehat{\rho }$ is the density matrix operator, $\widehat{M}$ is the dipole moment operator and ${\mathrm{\Gamma }}_{ij}$ are the matrix terms of the phenomenological operator $\widehat{\mathrm{\Gamma }}$. The steady-state density matrix ${\widehat{\rho }}^{\left(0\right)}$ is diagonal and its elements ${\rho }_{ii}^{\left(0\right)}$ represent the electron populations of the i-th energy levels $E_i$. ${\widehat{H}}_0={\sum }_{k=1}^4\hslash {\omega }_k|k⟩⟨k|$ is the free Hamiltonian of the system in the absence of the electric field for the four subbands, and ${\widehat{H}}_I=-\widehat{M}E\left(t\right)=-e\widehat{x}E\left(t\right)$ is the perturbation term that accounts for the light-matter interaction. The phenomenological operator $\widehat{\mathrm{\Gamma }}$ incorporates damping effects caused by interactions between electrons and phonons and electron-electron collisions. We assume that $\widehat{\mathrm{\Gamma }}$ is a diagonal matrix associated with the relaxation process and its elements are the inverse of the relaxation time that involves the |m⟩ state, ${\mathrm{\Gamma }}_{mm}=1/{\tau }_m$ with $m=1,2,3$ [62,63]. To simplify, we use the notation for the term ${\mathrm{\Gamma }}_{mm}={\mathrm{\Gamma }}_m$.

The master equation [Equation (7)] is solved using a perturbative method, in which the system is divided into an unperturbed part ${\widehat{H}}_0$ and a perturbative part ${\widehat{H}}_I$ that introduces slight corrections to the system dynamics. This approach yields a series solution, each perturbative order adding corrections to the total evolution. Consequently, the density matrix is expanded perturbatively in terms of the electric field as:

$$\widehat{\rho }\left(t\right)=\sum _{n=0}^{\infty }{\widehat{\rho }}^{\left(n\right)}\left(t\right).$$

(8)

Replacing into the Equation (7) and using the completeness relation ${\sum }_k|k⟩⟨k|=\mathbb{1}$, we have:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{\partial {\rho }_{ij}^{(n+1)}}{\partial t}}& ={\displaystyle \frac{1}{i\hslash }}\sum _{n=0}^{\infty }\left[\sum _{k=1}^4\right\{\hslash {\omega }_k\left({\delta }_{ik}^{\prime }{\rho }_{kj}^{(n+1)}-{\delta }_{kj}^{\prime }{\rho }_{ik}^{(n+1)}\right)\hfill \\ \hfill & -E\left(t\right)\left(M_{ik}{\rho }_{kj}^{\left(n\right)}-M_{kj}{\rho }_{ik}^{\left(n\right)}\right)\}-i\hslash {\mathrm{\Gamma }}_{ij}{\rho }_{ij}^{(n+1)}].\hfill \end{array}\end{array}$$

(9)

Here ${\delta }^{\prime }$ is the Kronecker delta (0 for $i\ne j$ and 1 for $i=j$). We can expand the different perturbative orders of $\rho$ as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{\partial {\rho }_{ij}^{\left(1\right)}}{\partial t}}& ={\displaystyle \frac{1}{i\hslash }}\sum _{k=1}^4\left\{\hslash {\omega }_k\left({\delta }_{ik}^{\prime }{\rho }_{kj}^{\left(1\right)}-{\delta }_{kj}^{\prime }{\rho }_{ik}^{\left(1\right)}\right)-E\left(t\right)\left(M_{ik}{\rho }_{kj}^{\left(0\right)}-M_{kj}{\rho }_{ik}^{\left(0\right)}\right)\right\}-i\hslash {\mathrm{\Gamma }}_{ij}{\rho }_{ij}^{\left(1\right)},\hfill \end{array}\end{array}$$

(10a)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{\partial {\rho }_{ij}^{\left(2\right)}}{\partial t}}& ={\displaystyle \frac{1}{i\hslash }}\sum _{k=1}^4\left\{\hslash {\omega }_k\left({\delta }_{ik}^{\prime }{\rho }_{kj}^{\left(2\right)}-{\delta }_{kj}^{\prime }{\rho }_{ik}^{\left(2\right)}\right)-E\left(t\right)\left(M_{ik}{\rho }_{kj}^{\left(1\right)}-M_{kj}{\rho }_{ik}^{\left(1\right)}\right)\right\}-i\hslash {\mathrm{\Gamma }}_{ij}{\rho }_{ij}^{\left(2\right)},\hfill \end{array}\end{array}$$

(10b)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{\partial {\rho }_{ij}^{\left(3\right)}}{\partial t}}& ={\displaystyle \frac{1}{i\hslash }}\sum _{k=1}^4\left\{\hslash {\omega }_k\left({\delta }_{ik}^{\prime }{\rho }_{kj}^{\left(3\right)}-{\delta }_{kj}^{\prime }{\rho }_{ik}^{\left(3\right)}\right)-E\left(t\right)\left(M_{ik}{\rho }_{kj}^{\left(2\right)}-M_{kj}{\rho }_{ik}^{\left(2\right)}\right)\right\}-i\hslash {\mathrm{\Gamma }}_{ij}{\rho }_{ij}^{\left(3\right)}.\hfill \end{array}\end{array}$$

(10c)

We assume, in the weak probe approximation, all electron population at $t=0$ is in the ground state, which requires that ${\rho }_{11}^{\left(0\right)}=1$, ${\rho }_{22}^{\left(0\right)}={\rho }_{33}^{\left(0\right)}={\rho }_{44}^{\left(0\right)}=0$ and ${\rho }_{ij}^{\left(0\right)}=0$ for $i\ne j$. The analytical solution for the matrix elements at different perturbative orders is obtained using a procedure similar to that described in Refs. [61,62].

The electronic polarization and susceptibility coefficients are related through the light polarization degree, induced by the electric field $E\left(t\right)$. This relation is derived by equating terms proportional to $exp\left(\pm {\omega }_{n}t\right)$ in the polarization expansion with those obtained from the perturbative expansion of the density matrix in the steady-state response ${\widehat{\rho }}^{\left(n\right)}\left(t\right)={\tilde{\rho }}^{\left(n\right)}\left(\omega \right)e^{-i{\omega }_{n}t}+{\tilde{\rho }}^{\left(n\right)}(-\omega )e^{i{\omega }_{n}t}$ with ${\omega }_n=\{0,\omega ,2\omega ,3\omega \}$, through the following equation [64]:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{P}}^{\left(n\right)}\left(t\right)=& {\displaystyle \frac{1}{V}}{\widehat{\rho }}^{\left(n\right)}\widehat{M}\hfill \\ \hfill =& {ϵ}_0\left({\chi }_{\omega }^{\left(1\right)}\tilde{E}{e}^{-i\omega t}+{\chi }_0^{\left(2\right)}{\tilde{E}}^2+{\chi }_{2\omega }^{\left(2\right)}{\tilde{E}}^2{e}^{-2i\omega t}+{\chi }_{\omega }^{\left(3\right)}{\tilde{E}}^2\tilde{E}{e}^{-i\omega t}+{\chi }_{3\omega }^{\left(3\right)}{\tilde{E}}^3{e}^{-3i\omega t}+\mathrm{c}.\mathrm{c}.\right),\hfill \end{array}\end{array}$$

(11)

where V is the volume of the system, ${\chi }_{\omega }^{\left(1\right)}$, ${\chi }_0^{\left(2\right)}$, ${\chi }_{2\omega }^{\left(2\right)}$, ${\chi }_{\omega }^{\left(3\right)}$ and ${\chi }_{3\omega }^{\left(3\right)}$ are the linear optical susceptibility, nonlinear optical rectification, second harmonic generation, third order linear susceptibility and third harmonic generation, respectively. This work focuses on the coefficients associated with optical rectification, second harmonic generation, and third harmonic generation.

Finally, we find the corresponding expressions per unit surface as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\chi }_0^{\left(2\right)}& ={\displaystyle \frac{4e^3{\sigma }_V}{{ϵ}_0}}{M}_{21}^2{\delta }_{21}\times {\displaystyle \frac{E_{21}^2\left(1+\frac{{\mathrm{\Gamma }}_2}{{\mathrm{\Gamma }}_1}\right)+{\hslash }^2\left({\omega }^2+{\mathrm{\Gamma }}_2^2\right)\left(\frac{{\mathrm{\Gamma }}_2}{{\mathrm{\Gamma }}_1}-1\right)}{\left({\left(E_{21}-\hslash \omega \right)}^2+{(\hslash {\mathrm{\Gamma }}_2)}^2\right)\left({\left(E_{21}+\hslash \omega \right)}^2+{(\hslash {\mathrm{\Gamma }}_2)}^2\right)}},\hfill \end{array}\end{array}$$

(12)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\chi }_{2\omega }^{\left(2\right)}& ={\displaystyle \frac{e^3{\sigma }_V}{{ϵ}_0}}\times {\displaystyle \frac{M_{21}{M}_{32}{M}_{31}}{\left(\hslash \omega -E_{21}-i\hslash {\mathrm{\Gamma }}_3\right)\left(2\hslash \omega -E_{31}-i\hslash {\mathrm{\Gamma }}_3/2\right)}},\hfill \end{array}\end{array}$$

(13)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\chi }_{3\omega }^{\left(3\right)}& ={\displaystyle \frac{e^4{\sigma }_V}{{ϵ}_0}}\times {\displaystyle \frac{M_{21}{M}_{32}{M}_{43}{M}_{41}}{\left(\hslash \omega -E_{21}-i\hslash {\mathrm{\Gamma }}_3\right)\left(2\hslash \omega -E_{31}-i\hslash {\mathrm{\Gamma }}_3/2\right)\left(3\hslash \omega -E_{41}-i\hslash {\mathrm{\Gamma }}_3/3\right)}},\hfill \end{array}\end{array}$$

(14)

where $E_{ij}=E_i-E_j$ is the intersubband transition energy, and ${\sigma }_V$ is the carrier density. In the above equations, the DMMEs are defined as:

$$M_{fi}=\int {\mathrm{\Psi }}_f^{\ast }\left(x\right)x{\mathrm{\Psi }}_i\left(x\right)dx(i,f=1,2,3,4),$$

(15)

where ${\delta }_{21}=\left|M_{22}-M_{11}\right|$, $M_{ij}$ ($M_{ii}$) are the off-diagonal (on-diagonal) DMMEs, respectively. The eigenfunctions ${\mathrm{\Psi }}_i\left(x\right)$ and ${\mathrm{\Psi }}_f\left(x\right)$ are obtained from Equation (2).

3. System and Parameters

The system consists of a sequence of five coupled GaAs quantum wells and six ${\mathrm{Al}}_x{\mathrm{Ga}}_{1-x}$As barriers, as illustrated in Figure 1b. Because of the small separation between subbands, close to optical phonon energies ($\hslash {\omega }_{LO}\sim 36$ meV in GaAs), this structure exhibits short excited states lifetimes ($\tau \sim 1$ ps) [65]. This short lifetime of electrons in the excited subbands primarily arises from electron-phonon interactions, the dominant mechanism for energy relaxation in this system.

For the calculations, the following parameters have been used from experimental results: an Al concentration of $x=0.3$, based on the $60\%:40\%$ distribution rule for the band offset between the conduction and valence bands, respectively. The barrier height is expressed as $V_0=\Xi (E_g^{AlGaAs}-E_g^{GaAs})=228$ meV, where $\Xi =0.60$ [66]. We consider well (barrier) width of 4 nm (3 nm). The effective electron mass for GaAs (AlGaAs) is $m_w^{\ast }=0.067m_0$ ($m_b^{\ast }=0.09m_0$), ${\sigma }_V=4\times {10}^{20}1/{\mathrm{m}}^3$, ${\mathrm{\Gamma }}_1=1/\left(2\mathrm{ps}\right)$, ${\mathrm{\Gamma }}_2=1/\left(1.5\mathrm{ps}\right)$, and ${\mathrm{\Gamma }}_3=1/\left(1.0\mathrm{ps}\right)$ [13,64,67,68].

4. Results and Discussion

In this section, we discuss the external field-mediated electro-optical control of the ISBT and DMMEs. To proceed systematically, we first analyze the influence of individual fields (electric and nIL) separately, followed by their combined effect. Next, we discuss cases where ISBT and wavefunction overlap are manipulated using purely electrical or optical means. Finally, we delve into the impact of simultaneous electric and laser fields on the NLO response of the MQWs system.

4.1. Influence of the External Fields on the DMMEs and the ISBT

Figure 1b shows the profile of the conduction band for a ${\mathrm{GaAs/Al}}_{0.3}{\mathrm{Ga}}_{0.7}$As symmetric MQW system composed of five wells and six barriers. The probability density associated with the first four states is also depicted, corresponding to eigenenergies $E_1$ to $E_4$. Without any external field ($F=0$ and ${\alpha }_0=0$), the initial symmetric confinement potential profile results in alternating symmetric and antisymmetric parity for the wave functions. The electrons in the ground state are located to a greater extent at the central quantum well center (see blue curve), with the probability density being maximum at this point and different from zero in all the wells, which indicates that there may be tunneling of the electrons. The maximum probability density localizes in the central well due to the electrons favoring regions of minimal confinement, in agreement with theoretical expectations.

Figure 2 shows the modifications to the potential due to the presence of the nIL field and the electric field. Figure 2a shows the potential modified by a nIL field without an electric field ($F=0$, black line) and the first four probability densities of the confined states. Under the presence of the nIL field, a new symmetric potential arises in the system of wells and barriers, which directly depends on the changes in the laser-dressing parameter ${\alpha }_0$, as shown in the figure (compare Figure 1b with Figure 2a). For ${\alpha }_0=4$ nm, the number of wells is modified in the system, going from five to four effective wells; note that in the initial system, at the center position ($x=0$), there was one well and the laser effect transforms this well into a barrier. Additionally, all wells become shallower and wider at the top and narrower at the bottom, leading to a noticeable separation in the energy levels. This behavior has also been previously reported by Lima et al. [69] for values of the laser parameter larger than the well width ${\alpha }_0>L/2$. As a result, the electrons are distributed mainly in the two central wells (see blue curve).

Figure 2b shows the potential profile modified by the effects of the electric field without the nIL field (${\alpha }_0=0$) and its corresponding probability densities. The electric field induces polarization in the system, tilts the potential, and shifts the charge distribution, breaking the inversion symmetry of the system. Due to this, the electrons experience a force that displaces them parallel to the fields and locates them in the wells on the left. Figure 2c shows the potential under the presence of both fields, the nIL (${\alpha }_0=4$ nm) and the electric field ($F=20$ kV/cm). We observe a combined effect in the modified potential: alterations in the shapes of the wells and barriers caused by the nIL field, and a tilt induced by the electric field. Note that, with the action of both fields, the MQW system undergoes a significant change in the electron confinement, which could alter its electronic and optical properties.

Figure 3a-d show the ISBT energies between the first four states as a function of the laser-dressing parameter and the electric field strength for the symmetric ${\mathrm{GaAs/Al}}_{0.3}{\mathrm{Ga}}_{0.7}$As MQW system illustrated in Figure 1. Note that the ISBT energies are scaled (i.e., divided by a factor of two and three $E_{31}/2,E_{41}/3$, respectively), which is a direct consequence of the quantization of energy levels in confined semiconductor structures. This scaling facilitates the analysis of the optical properties discussed later in this investigation. In each panel, the curves include the first four states as functions of the aforementioned external fields, which allows for a direct observation of the ISBT energy modification induced. Figure 3a and Figure 3d show the ISBT energies as a function of the electric field. Without the nIL field [Figure 3 a], a monotonic increase in the behavior of the ISBT energies is observed with increasing electric field strength (see Table 1). Furthermore, three ISBT energy intersections occur: $E_{21}=E_{31}/2,E_{21}=E_{41}/3$ and $E_{41}=3E_{31}/2$, in the electric field range $F\sim 7.5-12.5$, kV/cm, these factors come from the denominators of the equations. () and () to obtain a maximum in the corresponding susceptibilities. The intersections result from the coincidence between ISBT energies in a specific range of electric field values. This coincidence influences the dynamics of the optical response of the system, as will be discussed later. The insets in all panels of Figure 3 show that the energies do not present crossings or anticrossings. In contrast to Figure 3a, in Figure 3c the $E_{21}=E_{31}/2$ and $E_{21}=E_{41}/3$ intersections are observed, in this case for an electric field value of $F\sim 15$ kV/cm, instead of the $F\sim 7$ kV/cm of panel (a). This shift towards higher energy values arises because of the modifications in the potential profile induced by the nIL field. Figure 3b shows the ISBT as a function of the laser-dressing parameter without an electric field. The increase in the laser-dressing parameter induces a monotonic increasing behavior of the ISBT energies without intersections between them. In contrast, the application of an electric field, Figure 3 d, is again responsible for the occurrence of intersections between ISBT energy pairs: $E_{31}=2E_{31}/3$ and $E_{21}=E_{41}/3$, for laser-dressing parameter values of ${\alpha }_0\sim 2.5-4$ nm.

Taking into account the effects of external fields on the ISBT energies, we can now explore the manipulation of the DMMEs through these fields. The DMMEs as a function of the electric field strength without a nIL field (${\alpha }_0=0$) and with a nIL field (${\alpha }_0=4$) are shown in Figure 4a and Figure 4c, respectively. In Figure 4a, without the nIL field, a monotonic decrease is observed in the matrix elements $|M_{21}|,|M_{32}|,|M_{41}|$ and $|M_{43}|$ as the strength of the electric field increases. For $|M_{31}|$ and $|M_{22}-M_{11}|$, a non-zero value is obtained, due to the breaking of the symmetry inversion induced by the electric field, which enables transitions that were previously forbidden. With increasing electric field strength, the wave functions shift and deform, reaching at some point an optimal overlap ($F\sim 10$kV/cm) that maximizes the DMMEs $|M_{31}|$ and $|M_{22}-M_{11}|$. For higher field values, distortion and excessive shifting reduce the overlap, causing the matrix elements to decrease again. In contrast, the DMMEs $|M_{21}|,|M_{32}|,|M_{41}|$ and $|M_{43}|$ shows a monotonic decrease as the electric field strength increases, suggesting a lower overlap of the wave functions. In Figure 4c, the nIL field causes the matrix elements $|M_{21}|$ and $|M_{32}|$ to decrease slightly as the electric field increases (compared to Figure 4a). For DMMEs involving the fourth state ($|M_{41}|$ and $|M_{43}|$), a non-monotonic behavior is appreciable because this state, for high nIL field strengths (${\alpha }_0=4$ nm), can be coupled to the continuum states (cyan line in Figure 2c). Furthermore, the maximum value of the matrix elements $|M_{31}|$ and $|M_{22}-M_{11}|$ shifts towards higher values of the electric field ($F\sim 15$kV/cm) and remains almost constant from that value onward. All this suggests that the presence of the nIL field alters the overlapping wave functions and shifts the values of the electric field where maxima occur in the DMMEs.

Figure 4b and Figure 4d show the DMMEs as a function of the laser-dressing parameter, without an electric field ($F=0$) and with an electric field ($F=20$ kV/cm), respectively. In Figure 4b, with $F=0$, a slight decrease of the matrix elements $|M_{21}|,|M_{32}|,|M_{41}|$ and $|M_{43}|$ is observed as the nIL field increases. Furthermore, $|M_{31}|$ and $|M_{22}-M_{11}|$ remain constant at zero (the associated lines are superimposed on the plot), indicating that the nIL field slightly alters the overlapping of the wave functions and does not break the inversion symmetry, preserving its initial selection rules. In Figure 4d, the electric field $F=20$ kV/cm breaks the inversion symmetry of the system and the DMMEs increase as a function of the laser-dressing parameter, reaching a maximum value for ${\alpha }_0\sim 3$ nm, and slightly decreasing from that value onwards. The reason behind this is the displacement and overlapping of the wave functions produced by the effects of both fields.

Before calculating the susceptibility coefficients, it is convenient to study how these external fields alter the NLO response in the system. These changes are, in principle, associated with the so-called geometric factors, which refer to the terms involving the product between the DMMEs in each susceptibility coefficient. Figure 5a and Figure 5b show the geometric factor $M_{21}^2{\delta }_{21}$, where ${\delta }_{21}=|M_{22}-M_{11}|$, associated with the nonlinear optical rectification NOR [Equation ()], as a function of the electric field strength and the laser-dressing parameter, respectively. In Figure 5a, without the nIL field (solid line), it is observed that the geometric factor reaches its maximum value for an electric field of $F\sim 5$ kV/cm, decreasing from that point onward as the field increases (see Table 2). With an nIL field of ${\alpha }_0=4$ nm (dashed line), the maximum value shifts to higher field intensities and is reached for an electric field of $F\sim 7.5$ kV/cm. The reason for this lies in the non-monotonic variation of the DMMEs $M_{21}^2$ and ${\delta }_{21}$ with the electric field, and their product reaches a maximum value in the ranges where both are relatively large (see Figure 4a). Figure 5b demonstrates that, in the absence of an electric field (solid line), the geometric factor remains zero when varying the laser-dressing parameter because the nIL field does not break the inversion symmetry of the system so that the term ${\delta }_{21}$ involving the matrix elements $M_{11}$ and $M_{22}$ is zero (see Figure 4b). In contrast, in the presence of an electric field (dashed line), a non-monotonic behavior is observed in the curve associated with the geometric factor, reaching its maximum value for a nIL field of ${\alpha }_0\sim 3$ nm.

Figure 6a and Figure 6b show the geometric factor $M_{21}{M}_{32}{M}_{31}$, associated with second harmonic generation (SHG) [Equation ()], as a function of the electric field strength and the laser-dressing parameter, respectively. In Figure 6a, in the absence of the nIL field (solid line), it is observed that the geometric factor reaches its maximum value for an electric field of $F\sim 5$ kV/cm. On the other hand, when applying an nIL field to the system (dashed line), a reduction is observed for the maximum values, in addition to a shift towards an electric field value of $F\sim 10$ kV/cm. Figure 6b manifests that in the absence of an electric field, the geometric factor remains constant at zero by varying the laser-dressing parameter because the matrix element $M_{31}$ requires the inversion symmetry breaking induced by the electric field (see Figure 4b).

Figure 7a and Figure 7b show the geometric factor $M_{21}{M}_{32}{M}_{43}{M}_{41}$, associated with third harmonic generation (THG) [Equation (15)], as a function of the electric field strength and the laser-dressing parameter, respectively. In comparison to the geometric factors discussed above, this does not require the electric field-induced inversion symmetry breaking, due to the alternating symmetry and antisymmetric parity of the wavefunctions; therefore, the matrix elements involved are non-zero as a function of the external fields, as shown in Figure 4. In Figure 7a, the geometric factor exhibits a monotonic decreasing behavior as the intensity of the electric field increases, both in the absence of the nIL field (solid line) and in the presence of nIL field (dashed line). The main difference lies in the maximum values reached in $F=0$ kV/cm by this factor in the presence of the nIL field, while without this field, the maximum value is lower, with a more pronounced decrease as the electric field strength increases. Figure 7b shows the curves of the geometric factor as a function of the nIL field; in the absence of an electric field (solid line), this factor remains nearly constant for a laser-dressing parameter ${\alpha }_0$ ranging from 0 to $\sim 1$ nm. For larger values of the laser-dressing parameter, the geometric factor decreases monotonically. When an electric field of $F=20$ kV/cm is applied, the opposite behavior is observed, i.e., the geometric factor monotonically increases for ${\alpha }_0>1$ nm, being rather constant for lower values.

4.2. External Fields Control Nonlinear optical response

Figure 8 shows the NOR coefficient (${\chi }_0^{\left(2\right)}$) as a function of the incident photon energy for different values of the electric field strength and the nIL field. Panel (a) of this figure shows that in the absence of the nIL field (${\alpha }_0=0$), the peak associated with the NOR coefficient reaches its maximum value when $F=5$ kV/cm. As the electric field increases, this peak experiences a blueshift, and its maximum value decreases. The energy shift is due to the energy of the intersubband transition $E_{21}$, involved in the NOR coefficient [Equation (12)], which exhibits a monotonic increasing behavior (black line in Figure 3a). The above indicates that the increase in the electric field generates a separation between the intersubband levels, shifting the resonance peak towards higher energies. Additionally, the decrease in the maximum value occurs because the electric field alters the shape and position of the quantum state wave functions, thereby modifying their overlap. This effect is illustrated by the solid line in Figure 5a, which shows that the geometric factor reaches a maximum value at approximately $F=5$ kV/cm before decreasing as the field increases. Finally, it is important to note that for $F=0$, the NOR peak is zero since the electric field is required to generate the non-linear response of the system.

For comparison with panel (a), panel (b) of Figure 8 illustrates the NOR coefficient as a function of the incident photon energy for the same four electric field strength values, but in the presence of the nIL field (${\alpha }_0=4$ nm). The first remark when looking at panels (a) and (b) is that the magnitude of the maximum value reached by the NOR coefficient in $F=5$ kV/cm is lower in the presence of the nIL field. Under this setup, a blueshift in the NOR peak and an increase in its maximum value for $F=10$ kV/cm; however, from that value, the maximum value decreases. The reason for this blueshift is that, under the influence of the nIL field, the geometric factor associated with the NOR (see Figure 5a, dashed line) experiences a shift in its maximum value and a decrease in its magnitude. Finally, Figure 8c presents the NOR peak as a function of the photon energy for different values of the laser-dressing parameter ${\alpha }_0$ in the presence of an electric field of $F=20$ kV/cm. In contrast to Figure 8a and Figure 8b, here the blueshift of the peak is less pronounced, and its magnitude increases as ${\alpha }_0$ grows. The nIL field-induced confinement produces a slight blueshift in the $E_{21}$ intersubband energy as ${\alpha }_0$ increases (black line in Figure 3d), influencing, therefore, on the separation between levels, which turns into a less marked shift towards higher energies. As previously discussed, the curve associated with the geometric factor (Figure 5b, dashed line) determines the behavior of the NOR peak maximum value. The combined effect of the electric field and the nIL field offers a tunable mechanism that modifies both the energy level separation and the overlap of the wave functions. These variations dynamically tune the magnitude and spectral position of the nonlinear optical rectification.

Figure 9 shows the SHG coefficient (${\chi }_{2\omega }^{\left(2\right)}$) as a function of the incident photon energy for different values of the electric field and the nIL field. In Figure 9a, when ${\alpha }_0=0$ (absence of nIL field), the SHG coefficient exhibits a blueshift with the photon energy, and when it reaches its maximum value ($F=7$ kV/cm, pink curve) the magnitude decrease as the electric field intensity increases. This blueshift is due to the monotonic increasing behavior of the intersubband transition energies $E_{21}$ and $E_{31}$ with the electric field strength (see Figure 3a, black and red curves). Furthermore, we observe that at $F=7$ kV/cm, the two characteristic peaks of the SHG coefficient are merged into a single peak with the highest magnitude; however, as the electric field increases, these peaks move apart. This separation influences the maximum value of the SHG coefficient, which describes the efficiency of the process and depends on both the intersubband transition energies and the incident photon energy. In contrast to Figure 8, the maximum value of the SHG coefficient is not governed exclusively by the behavior of the geometric factor $M_{21}{M}_{32}{M}_{31}$, but depends mainly on the approach or intersections between the ISBT energies. When the peaks move apart, the energies $E_{21}$ and $E_{31}$ become more distinguishable, indicating a change in the electronic structure of the multiple quantum well system induced by the electric field. This phenomenon is clearly shown from the variations in the level spacing (see Figure 3a, black and red curves): at $F=5$ kV/cm, the transition energies are close (even crossing $E_{21}=E_{31}/2$ at $F=7$ kV/cm, pink circle Figure 3a), as the electric field increases, they start to move apart. The moving apart of the peaks reduces the efficiency of SHG since the mismatch between resonances prevents the incident photon energy from being transferred efficiently, producing a notable decrease in the maximum coefficient value. Nevertheless, as the peaks move apart, SHG becomes more selective concerning the photon incident frequencies, implying that the system will respond strongly only at specific frequencies corresponding to the individual ISBT energies.

In Figure 9b, the presence of the nIL field (${\alpha }_0=4$ nm) causes a blueshift and an increase in the magnitude of the SHG coefficient as the electric field strength increases, reaching its maximum value at $F=16$ kV/cm (pink curve), from this point on, the coefficient decreases. The blueshift is associated with an increase in the ISBT energies with the electric field strength (see Figure 3c). On the other hand, in contrast to Figure 9a, the two characteristic peaks of the SHG coefficient are initially separated and closer as the electric field strength increases. When an approach or intersection occurs between the intersubband transition energies $E_{21}$ and $E_{31}$ ($F=16$ kV/cm, see Figure 3c pink circle), the two characteristic peaks of the SHG coefficient merge into a single, more intense peak. Physically, this indicates that the incident photon energy coincides, at some point, with the intersection between the ISBT energies, which maximizes the SHG efficiency under this resonance condition. In this scenario, the incident photon energy is transferred more effectively to the second harmonic generation process. The induced coincidence, through the two external fields, between the ISBT energies involved in the SHG leads to a significant increase in the magnitude of the second-harmonic signal. Finally, Figure 9c shows that under a fixed electric field of $F=20$ kV/cm, both a blueshift and an increase in the maximum value of the SHG coefficient occur as the laser-dressing parameter increases. Initially, the characteristic peaks of the SHG coefficient are separated. However, as the laser-dressing parameter increases, these peaks become closer. In this case, the ISBT energies $E_{21}$ and $E_{31}$ do not cross (see Figure 3d), but there is an approach between them as the nIL field increases; for that reason, we observe at ${\alpha }_0=4$ nm the highest magnitude of the SHG coefficient.

Figure 10a and Figure 10b show the THG coefficient (${\chi }_{3\omega }^{\left(3\right)}$) as a function of the incident photon energy for different values of the electric field strength in the absence (${\alpha }_0=0$) and presence (${\alpha }_0=4$ nm) of the nIL field, respectively. In Figure 10a, for electric field intensities ranging in $0-10$ kV/cm, the THG coefficient shows a blueshift together with a coalescence of the three characteristic peaks and increase in the magnitude. The coefficient for $F=20$ kV/cm breaks the previous trend, and undergoes a notable decrease in magnitude and a splitting of the peaks. The blueshift is due to the monotonic increase of the ISBT energies $E_{21}$, $E_{31}$, and $E_{41}$ (see Figure 3a). The coalescence and splitting of the three peaks is attributed to the fact that, in the range of $F\sim 7.5-12.5$ kV/cm, these transition energies approach to each other considerably (even presenting intersections between pairs such as $E_{21}=E_{31}/2$, $E_{21}=E_{41}/3$ and $E_{31}=2E_{41}/3$, but no coincidence between the three energies in a single point is found), whereas for $F>12.5$ kV/cm $E_{31}$ and $E_{41}$ remain close and $E_{21}$ separates, originating the appearance of two peaks close together and one further away. As previously discussed, the behavior of the geometrical factor $M_{21}{M}_{32}{M}_{43}{M}_{41}$ as a function of F (see Figure 7a, dashed line) does not explain the maximum values reached by the THG coefficient, which are mainly related to the proximity and intersections between the ISBT energies. In particular, the maximum value observed of the THG coefficient at $F=10$ kV/cm indicates that for this specific electric field value, the ISBT energies $E_{21}$, $E_{31}$ and $E_{41}$ are closer, which favors the efficiency of the THG process. Whereas for $F=20$ kV/cm, the separation of the peaks reduces the magnitude of the THG coefficient and confers selectivity concerning the incident photon frequencies. In comparison with Figure 10a, in Figure 10b the presence of the nIL field causes the characteristic THG peaks to be initially more separated, preserving the afore discussed blueshift but causing the occurrence of the maximum value at $F=20$ kV/cm. These results indicate that, by applying the nIL field, it is possible to adjust the electric field values at which approaches or intersections occur between ISBT energies and, in this way, determine the optimal combination of external fields that improves the THG efficiency. Indeed, these results demonstrate that the efficiency of the THG process can be improved when external fields induce coincidences between the ISBT energies. Moreover, external fields can also separate the ISBT energies and make the THG process highly selective with incident photon frequencies.

Figure 10c and Figure 10d show the THG coefficient as a function of the incident photon energy for different values of the laser-dressing parameter in the absence ($F=0$) and presence ($F=20$ kV/cm) of the electric field, respectively. In Figure 10c, as the nIL field increases, a blueshift is observed in the THG coefficient, accompanied by a progressive separation of the three characteristic peaks and a decrease in their magnitude. This blueshift is attributed to the monotonic increase of the ISBT energies $E_{21}$, $E_{31}$ and $E_{41}$ as a function of the nIL field (see Figure 3c). The separation of the peaks significantly influences the magnitude of the THG coefficient, suggesting that, in the multiple quantum well system, increasing the nIL field reduces the THG efficiency, although it improves the selectivity concerning the incident photon frequencies. This behavior could be useful in applications that require selective interaction with specific frequencies. On the other hand, Figure 10d shows that, in the presence of an electric field, the THG coefficient undergoes a weaker blueshift and its magnitude increases considerably as the nIL field increases. This small shift is due to the slight growth of the ISBT energies with the nIL field (see Figure 3d). For an electric field of $F=20$ kV/cm, the THG coefficient’s characteristic peaks align with the ISBT energies’ dependence on the laser-dressing parameter. We observe that $E_{31}$ and $E_{41}$ are very close, while $E_{21}$ is separated (see Figure 3d), resulting in two closely spaced peaks and one well-defined peak. This trend persists until the laser-dressing parameter reaches approximately $1.5$ nm. Beyond this value, an intersection occurs between the energies $E_{31}$ and $2E_{41}/3$, causing the two adjacent peaks to merge into a single peak. However, the peak corresponding to $E_{21}$ remains isolated. Notably, under these conditions, the magnitude of the THG coefficient increases substantially (see blue curve). A similar effect emerges at ${\alpha }_0=4$ nm, where another intersection arises between $E_{21}$ and $E_{41}/3$. All three ISBT energies converge at this point, leading to a higher enhancement in the magnitude of the THG coefficient.

5. Conclusions

In this research, we theoretically study the modulation of the non-linear optical response in a ${\mathrm{GaAs/Al}}_{0.3}{\mathrm{Ga}}_{0.7}$As symmetric multiple quantum well system under the simultaneous application of an IL field and an electric field. On the one hand, the nIL field modifies the electron confinement potential without altering the initial symmetry of the system; the so-induced changes affect the overlap of the wave functions and the electronic structure, leading to variations in the spacing between subbands and, consequently, in the ISBT energies. On the other hand, the electric field tilts the confinement potential. It alters the electronic distribution, breaking the inversion symmetry and enabling transitions between states previously forbidden by the selection rules. In addition, this field modifies, analogously to the nIL field, both the wave functions’ overlap and the separation between subbands.

To analyze the non-linear optical response, we determine the impact of each field on the ISBT energies and the wave function overlap, quantified by the geometrical factor. The results show that the simultaneous application of both external fields enables the manipulation of both the overlap of the wave functions and the coupling or intersection between the ISBT energies, which are determining factors in the nonlinear optical response of the system. The behavior of the geometrical factor as a function of the nIL field and the electric field determines the efficiency of the non-linear optical rectification; the system reaches a high efficiency when this factor reaches its maximum value. Additionally, accidental coincidence, induced by external fields, between pairs of ISBT energies improves the efficiency in the second and third harmonic generation. Instead, the ISBT energies separation confers selectivity to the system concerning the frequencies of the incident photon. Our results show the possibility of tuning the wavefunction overlap and ISBT energies by applying both external fields, demonstrating a dynamically tunable control of the nonlinear optical response in the multiple quantum well system. These findings might be employed in the design of optoelectronic devices with tunable optical properties.