Electroabsorption by confined excitons with Gaussian interaction potential
Yuriy D. Sibirmovsky, Ivan S. Vasil'evskii, Nikolay I. Kargin
aa r X i v : . [ phy s i c s . c o m p - ph ] J a n E LEC TROABSORPTION B Y C ONFINED E XC ITONS WITH G AUSSIAN I NTER ACTION P OTENTIAL
Yuriy D. Sibirmovsky, Ivan S. Vasil’evskii, Nikolay I. Kargin
Nanoengineering in Electronics, Spintronics and Photonics InstituteNational Nuclear Research University "MEPhI"Moscow 115409, Russian Federation [email protected]
January 18, 2021 A BSTRACT
We consider the effects of electron-hole interaction, 2D confinement and applied electric field on di-rect allowed transitions in III-V semiconductors, with InGaAs as a study case. Instead of Coulombinteraction, we use Gaussian potential. It is finite at the origin and has a finite effective range, whichallows for a more efficient numerical solution of Schrödinger equation. Yet, we can expect elec-troabsorption phenomena to remain qualitatively similar to the ones observed for Coulomb excitons.Moreover, we use variation of parameters to fit both position and magnitude of the first absorptionpeak to the Coulomb case.We combine and compare several numerical and approximate methods, including spectral expansion,finite differences, separation of variables and variational approximation.We find that separation of variables approach works only for quantum well widths smaller than theexciton radius. After separation of variables, finite difference solution of the resulting interactionequation gives a much better agreement with the full spectral solution than naive variational approx-imation.We observe that electric field has a critical effect on the magnitudes of exciton absorption peaks,suppressing previously allowed transitions and enhancing forbidden ones. Moreover, for excitedstates, initially suppressed transitions are enhanced again at higher field strengths. K eywords Schrödinger equation · absorption · electroabsorption · excitons · Stark effect · quantum wells · spectralmethods · finite difference · variational approximation Electroabsorption and electrorefraction phenomena play a crucial role in microwave photonics and other areas. Tocreate compact and efficient devices, such as electro-optic modulators, one needs to utilize the properties of semicon-ductor nanostructures, such as quantum wells (QWs) and superlattices. The research on the influence of electric fieldon the optical properties of such low dimensional structures had started more than 30 years ago [1, 2, 3, 4, 5, 6] withthe term "Quantum-Confined Stark Effect" being coined by Miller et al. in 1984 [1] to emphasize the wide range ofphenomena arising there, compared to the ordinary atomic Stark Effect.Molecular beam epitaxy allows one to produce a wide range of QW-based heterostructures, combining dozens of lay-ers of varying thickness and chemical composition. For microwave photonics and fiber communication systems oneof the most popular QW materials is In x Ga − x As with barriers made of In y Al − y As, while InP is used as substrate.Increasingly often, more complicated InGaAsP compounds are used to achieve independent control of bandgap and lat-tice constant [7, 8], while special attention is paid to the role of asymmetry in increasing efficiency of electroabsoprtionor electrorefraction [9, 10].
PREPRINT - J
ANUARY
18, 2021Table 1: Constants and material parameters used in calculation
Name Parameter Value Units - ~ m . eV · nm Fine structure constant α = e πε ~ c / . -Bohr radius r = ~ m · πε e . nmHomogeneous broadening Γ 0 . eVIn . Ga . As parameters [8]:Bandgap E g . eVOptical matrix parameter E p . eVBackgr. dielectric const. ε b = n b . -Electron effective mass m e m . -Heavy hole effective mass m h m . -Reduced exciton mass m eh m = m e m h m ( m e + m h ) . -Exciton Bohr radius r ex = r ε b m m eh . nmExciton Rydberg energy Ry = ~ m eh r ex . eVFor such complicated multilayered heterostructures, electro-optic effects depend on the layer composition in a non-trivial way, and design of more efficient devices necessarily involves numerical modeling and simulation. While therehave been a multitude of studies dedicated to modeling absorption and refraction (with or without applied electric field)of quantum-confined semiconductor structures [11, 12, 13, 14], this is still an on-going research field [15, 16, 17]. Thephotonics industry requires a combination of different methods: fast algorithms involving various approximations areneeded for real time multi-parameter optimization of heterostructures, while slower, but more accurate algorithmsshould be used to verify the results quantitatively.However, there are several challenges involved with the task of modeling electroabsorption. The most important effectsare electron-hole interaction and quantum confinement. Secondary effects include non-parabolicity, scattering, bandmixing or coupling, strain etc. While complex methods which take most of these features into account have alreadybeen developed, they are usually quite complicated, require a lot of computational power, and are not straightforwardto implement. In fact, the influence of Coulomb interaction even on bulk Franz–Keldysh effect is still a subject ofinvestigation both theoretical [18] and applied [19].Which is why in this paper, we consider a rather basic problem: a single exciton, confined in a quantum well andsubjected to electric field in the direction, perpendicular to the QW plane. In addition, we use a modified interactionpotential, which, while qualitatively similar to the Coulomb one, is finite at the origin, and has a finite effectiverange. This leads to a more efficient numerical solution, which in turn allows us to compare different methods andapproximations, as well as study the electroabsorption phenomena, and the role played by quantum confinement. Thispotential can also be useful in modeling the phenomenon of static screening, which becomes important at higher lightintensities and/or doping levels. It should be noted that simplified interaction potentials are a common approach toprobe complex phenomena, where analytical solutions are not readily available otherwise, see for example [20] wheredelta-function potential is used to model exciton interaction.The main goals of this work are as follows. The first one is to test the widely used separation of variables approachand its particular case, variational approximation for the electron-hole interaction wavefunction. The second is to testthe spectral method, which, compared to grid methods (finite differences or finite elements) leads to much smaller andmore manageable matrices. The third is to obtain a qualitative picture of electroabsorption in a quantum well, whichshould be quite similar to the actually observed Coulomb exciton case.We used R language [21] for most of the calculations and for plot generation, while Wolfram Mathematica [22] wasused for the spectral expansion method and for evaluating certain special functions.We will present all the absorption spectra in absolute units of µ m − , which will allow us to quantitatively compare themagnitude for different cases. The important constants and parameters are summarized in Table 1. We use a two-bandparabolic effective mass envelope function approximation in position space. For clarity, in his work we don’t take intoaccount polarization dependence of QW absorption and only consider the heavy-hole exciton.2 PREPRINT - J
ANUARY
18, 2021
For bulk (3D) exciton with no applied electric field, the absorption spectrum (due to direct allowed transitions) can beexpressed analytically [23, 24]: α ( ~ ω ) = 43 α n b ~ m E p r ex · ℑ " + ∞ X n =1 /r ex · ~ wE n [ E n − ( ~ w ) ] n + m eh ~ w G ( ~ w ) E n = E g − Ryn G ( ~ w ) = g s RyE g − ~ w ! + g s RyE g + ~ w ! − g s RyE g ! (1) g ( ξ ) = ln ξ + φ ( ξ ) + 12 ξ ~ w = ~ ω + i Γ Where φ ( ξ ) - digamma function. For a hypothetical 2D exciton confined in a quantum well of with L , there’s a similarexpression [25, 24]: α ( ~ ω ) = 43 α n b ~ m E p L · ℑ + ∞ X n =1 /r ex · ~ w ˜ E n h ˜ E n − ( ~ w ) i ( n − / + m eh ~ w F ( ~ w ) ˜ E n = ˜ E g − Ry ( n − / F ( ~ w ) = f s Ry ˜ E g − ~ w ! + f s Ry ˜ E g + ~ w ! − f s Ry ˜ E g ! (2) f ( ξ ) = ln ξ + φ (cid:18) ξ + 12 (cid:19) Where ˜ E g takes into account the increase of transition energy due to quantum size effect for electrons and holes.We may notice the similarity between (1) and (2), which highlights the fact that strongly confined excitons have groundstate binding energy times greater and peak absorption r ex /L times greater than their bulk counterparts.Using these expressions and taking the real part instead of imaginary, we could also calculate refractive index, with noneed to directly apply Kramers-Kronig relations.If we discard the electron-hole interaction, and further take Γ = 0 , Eq. (1) reduces to the well known square rootdependence [27]: α ( ~ ω ) = 2 √ α n b (cid:18) ~ m (cid:19) − / (cid:18) m eh m (cid:19) / E p ~ ω p ~ ω − E g (3)While (2) reduces to the step function θ ( ξ ) dependence (in the following expression we only take into account theground state) [27]: α ( ~ ω ) = 2 π L α n b m eh m E p ~ ω θ ( ~ ω − ˜ E g ) (4)These expressions are important as the way to check the numerical results for modified interaction potentials.3 PREPRINT - J
ANUARY
18, 2021
When no confinement is present (3D) or in the case of very strong confinement (2D) and no applied electric field,Schrödinger equation reduces to a single variable case, with either spherical or cylindrical symmetry. ˆ Hψ = − ~ m eh ∆ r ψ ( r ) − RyU ( r ) ψ ( r ) = Eψ ( r ) (5) U ( r ) = r ex r Coulomb interaction a exp (cid:18) − r b r ex (cid:19) Gaussian interactionFor the Gaussian interaction potential, introduced here, a is the depth parameter and b the width parameter. Thiscombination of parameters allows us to vary both the interaction strength and its range.Reduction to (effectively) one dimension allows us to apply a finite difference (FD) scheme with small and uniformstep size, which is the most simple way to account for any interaction potential. However, in these cases, one needsto be careful when deriving the FD approximations for both kinetic and potential energy operators. In general, theHamilton operator matrix will not be Hermitian.Let the grid be defined by r j = j ∆ r with j = 0 , , . . . , N r and N r ∆ r = R . We summarize our schemes in Table2. The 2D scheme mostly follows [12], except for the Coulomb potential approximation. The 3D scheme is derivedindependently here.Let’s introduce ˆ S as a diagonal matrix with elements s j and ˆ H as a tridiagonal matrix with elements defined by Table2. Then equation (5) transforms into a generalized matrix eigenvalue problem: ˆ S ( ˆ H + E g ) ~ψ = E ˆ S ~ψ (6)It’s easy to see that for the schemes in Table 2, the matrix ˆ S ˆ H is Hermitian for the 2D exciton, but non-Hermitian inthe 3D case. However, the latter asymmetry is not significant, and only requires us to take the real parts of eigenvaluesand eigenvectors as a precaution.After solving (6) and obtaining eigenvalues E n and eigenvectors ~ψ ( n ) where n = 1 , . . . , N r + 1 , it’s convenient torenormalize the eigenvectors in the following way: ~φ = (cid:8) √ s ψ , √ s ψ , . . . , √ s N r ψ N r (cid:9)qP N r j =0 s j ψ j (7)Note that the FD approximations at j = 0 for operators with a singularity at the origin are obtained here by averagingwith the appropriate weight function from to ∆ r/ .The resulting expressions for the absorption spectrum will be: α D ( ~ ω ) = 43 α n b ~ m E p L · ℑ " N r +1 X n =1 / (∆ r ) · | φ ( n )0 | · ~ wE n [ E n − ( ~ w ) ] (8) α D ( ~ ω ) = 43 α n b ~ m E p ∆ r · ℑ " N r +1 X n =1 / (∆ r ) · | φ ( n )0 | · ~ wE n [ E n − ( ~ w ) ] (9)The continuous part of the spectrum is modeled by the contribution of discrete FD states with E n ≥ E g , providedthat R is large enough so that E n +1 − E n < Γ and ∆ r is small enough so that E N r +1 > ~ ω max , where ~ ω max isthe largest photon energy we are considering. Verification of the proposed schemes and formulas (8-9) is provided inFigure 1(a,b), where the numerical results are compared with the exact expressions (1-2). For the strongly confinedcase we take ˜ E g = E g for convenience.To compare relative magnitudes of the exciton absorption peaks, we introduce a dimensionless parameter, proportionalto the transition oscillator strength: 4 PREPRINT - J
ANUARY
18, 2021Table 2: Finite difference approximations to operators in 2D and 3D
Operator 2D 3D k ψ k R + ∞ πr | ψ ( r ) | dr R + ∞ πr | ψ ( r ) | dr ∆ r ψ r ∂∂r (cid:18) r ∂ψ∂r (cid:19) r ∂∂r (cid:18) r ∂ψ∂r (cid:19) FD approx. 2D 3D (∆ r ψ ) r ) ( − ψ + ψ ) 6(∆ r ) ( − ψ + ψ )(∆ r ψ ) j ≥ r ) (cid:20)(cid:18) − j (cid:19) ψ j − − ψ j + (cid:18) j (cid:19) ψ j +1 (cid:21) r ) "(cid:18) − j (cid:19) ψ j − − (cid:18) j (cid:19) ψ j + (cid:18) j (cid:19) ψ j +1 ψr (cid:19) r ψ r ψ (cid:18) ψr (cid:19) j ≥ ψ j j ∆ r (cid:16) e − β r ψ (cid:17) j ≥ e − β j (∆ r ) ψ j k ψ k π (∆ r ) P N r j =0 s j | ψ j | π (∆ r ) P N r j =0 s j | ψ j | s
18 124 s j ≥ j j + 112 PREPRINT - J
ANUARY
18, 2021 −5 0 5 10 . . . . . . Photon energy − E g , Ry A b s o r p t i on , m m - Exact, G = 2 meVFinite difference, R = 40 r ex , D r = 0.01 r ex −10 −5 0 5 10 Photon energy − E g , Ry A b s o r p t i on , m m - Exact, G = 2 meVFinite difference, R = 40 r ex , D r = 0.01 r ex (a) (b)Figure 1: Excitonic absorption spectra for (a) 3D (bulk) and (b) 2D (strongly confined) excitons ( L is taken equal to r ex ), as given by exact formulas (1-2) and finite difference calculations. −5 0 5 10 . . . . . Photon energy − E g , Ry A b s o r p t i on , m m - Coulomb, G = 0.5 meVGaussian, a = 0.286 , b = 0.964No interaction −10 −5 0 5 10 Photon energy − E g , Ry A b s o r p t i on , m m - Coulomb, G = 0.5 meVGaussian, a = 0.0683 , b = 0.276No interaction (a) (b) − . − . − . − . . Distance from the origin, r ex E ne r g y , R y Coulomb U(r)Gaussian U(r)Coulomb E , y Gaussian E , y − − − − − Distance from the origin, r ex E ne r g y , R y Coulomb U(r)Gaussian U(r)Coulomb E , y Gaussian E , y (c) (d)Figure 2: Comparison between Coulomb and Gaussian (with parameters given in Table 3) excitons. Absorptionspectra for (a) 3D and (b) 2D excitons. Potentials near the origin and ground state wavefunctions for (c) 3D and (d)2D excitons. 6 PREPRINT - J
ANUARY
18, 2021Table 3: Ground state parameters for 3D and 2D excitonsa b Binding energy E , Ry Oscillator strength Q Coulomb 3D - - -1 1Coulomb 2D - - -4 8Gaussian 3D 0.286 0.964 -0.999 1.034Gaussian 2D 0.0683 0.276 -3.99 7.96 Q Dn = 4 r ex (∆ r ) · | φ ( n )0 | (10) Q Dn = 6 r ex (∆ r ) · | φ ( n )0 | By varying the depth and width parameters a, b of the Gaussian interaction potential, it is possible to obtain differentnumbers of bound states, as well as binding energy and amplitude for the ground state. Empirically, we found thata reasonably good agreement with both 3D and 2D Coulomb exciton absorption edges is achieved for the valuespresented in Table 3. Figures 2(a,b) show the absorption spectra comparison, this time with a small broadening
Γ = 0 . meV to reduce the contribution from the above-bandgap states, which is present in the 3D case.In both cases, Gaussian potential supports a single bound state with binding energy and oscillator strength provided inTable 3. The values were calculated with R = 40 r ex and ∆ r = 0 . r ex , and the provided digits are correct.From Figures 2(c,d) we can see that the wavefunction for the Gaussian potential has zero slope at the origin, whichis due to the fact that the potential is finite everywhere. However, for the 2D case, the exponential Coulomb solutionagrees very well with the numerical one for the Gaussian potential.As for the continuous part of the spectrum, due to the short range of the Gaussian potential, there’s no significantSommerfeld enhancement, and this part strongly resembles the no interaction dependence in 3D and (especially) 2Dcases. Let us consider the problem of interacting electron-hole pair under applied electric field F , and confined in the z direction in a QW of width L . Then Schrödinger equation in relative coordinates will have the form (where we haveexplicitly added the bandgap energy): − ~ m eh (cid:20) ∂ ψ∂z + 1 r ∂∂r (cid:18) r ∂ψ∂r (cid:19)(cid:21) + eF zψ − RyU ( z, r ) ψ + E g ψ = Eψ (11) − L ≤ z ≤ L, ≤ r < + ∞ This is a simplified picture, which ignores the contribution of center-of-mass confinement. However, the latter is notaffected by electric field and thus is not of interest for the present work. It will be considered in further publications.Let us expand the wavefunction of equation (11), using the following basis: ψ ( z, r ) = N z X n z =1 N r X n r =1 C n z n r f n z ( z ) g n r ( r ) (12) f n z ( z ) = 1 √ L sin h πn z (cid:16) − zL (cid:17)i g n r ( r ) = 1 √ πRJ ( γ n r ) J (cid:16) γ n r rR (cid:17) PREPRINT - J
ANUARY
18, 2021Where J , J are Bessel functions, γ n r - zeros of J ( ξ ) , and we now assume ≤ r ≤ R where R ≫ r ex . This choiceof basis functions allows us to describe both the discrete and continuous parts of the absorption spectrum, similar tothe finite difference method. Normalization conditions for the basis functions are: Z L − L dzf n z ( z ) f n ′ z ( z ) = δ n z ,n ′ z Z R πrdrg n r ( r ) g n ′ r ( r ) = δ n r ,n ′ r Let us now find the matrix elements of the operators in equation (11) in the chosen basis: (cid:28) f n z (cid:12)(cid:12)(cid:12)(cid:12) − ~ m eh ∂ ∂z (cid:12)(cid:12)(cid:12)(cid:12) f n ′ z (cid:29) = ~ π n z m eh L δ n z ,n ′ z (13) (cid:28) g n r (cid:12)(cid:12)(cid:12)(cid:12) − ~ m eh r ∂∂r (cid:18) r ∂∂r (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) g n ′ r (cid:29) = ~ γ n r m eh R δ n r ,n ′ r (14) (cid:10) f n z | eF z | f n ′ z (cid:11) = 8 eF Lπ (cid:16) − ( − n z + n ′ z (cid:17) ( n z − n ′ z ) n z n ′ z (15) (cid:28) f n z (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) − z b r ex (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) f n ′ z (cid:29) = √ π ν (cid:2) V n z − n ′ z ( ν ) − V n z + n ′ z ( ν ) (cid:3) (16) ν = Lbr ex V n ( ν ) = cos πn (cid:18) − π n ν (cid:19) h erf (cid:16) ν + i πn ν (cid:17) + erf (cid:16) ν − i πn ν (cid:17)i (17) (cid:28) g n r (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) − r b r ex (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) g n ′ r (cid:29) = 2 J ( γ n r ) J ( γ n ′ r ) Z tdte − η t J ( γ n r t ) J ( γ n ′ r t ) (18) t = rR , η = Rbr ex The integral in (18) can’t be evaluated in analytical form, however, using the condition R ≫ r ex , we can claim that η ≫ , then it follows: Z tdte − η t J ( γ n r t ) J ( γ n ′ r t ) ≈ Z + ∞ tdte − η t J ( γ n r t ) J ( γ n ′ r t ) From [26] 6.633.2 the latter integral evaluates to (where I is modified Bessel function): Z + ∞ tdte − η t J ( γ n r t ) J ( γ n ′ r t ) = 12 η exp − γ n r + γ n ′ r η ! I (cid:18) γ n r γ n ′ r η (cid:19) Then it follows: (cid:28) g n r (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) − r b r ex (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) g n ′ r (cid:29) ≈ η J ( γ n r ) J ( γ n ′ r ) exp − γ n r + γ n ′ r η ! I (cid:18) γ n r γ n ′ r η (cid:19) (19)8 PREPRINT - J
ANUARY
18, 2021Gaussian form of the interaction potential allows us to express its matrix element as a product of two previously foundmatrix elements: (cid:10) f n z g n r | U ( z, r ) | g n ′ r f n ′ z (cid:11) = 1 a (cid:28) f n z (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) − z b r ex (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) f n ′ z (cid:29) (cid:28) g n r (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) − r b r ex (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) g n ′ r (cid:29) (20)Collecting all the matrix elements, we again obtain a matrix eigenvalue problem, with matrix dimensions being N z N r × N z N r . To be able to evaluate special functions with good accuracy and precision, we used Wolfram Mathe-matica for this calculation. After solving the eigenvalue problem, we obtain the set of eigenvectors and eigenvalues: n E λ , ~ψ λ o , λ = 1 , , . . . , N z N r Then the absorption spectrum can be expressed as: α ( ~ ω ) = 43 α n b ~ m E p L · ℑ " N z N r X λ =1 /r ex · Q λ · ~ wE λ [ E λ − ( ~ w ) ] (21)Where we again introduce a dimensionless oscillator strength, which now has the form: Q Dλ = r ex R N z X n z =1 N r X n r =1 C ( λ ) n z n r J ( γ n r ) sin πn z ! (22) Q Dλ = r ex L Q ( λ )2 D Note that nothing prevents us from making the parameters a and b functions of L to fit the Coulomb absorption edge.However, lacking the full numerical solution for the latter case, we will use the separation of variables approach forthe fit, which will be presented in the following section. Let’s assume that variables in equation (11) can be separated. This could be achieved by using the average of interac-tion potential over the z variable. Then we have: ψ ( z, r ) = ξ ( z ) ζ ( r ) − ~ m eh ∂ ξ∂z + eF zξ + E g ξ = E z ξ (23) − ~ m eh r ∂∂r (cid:18) r ∂ζ∂r (cid:19) − RyU [ ξ ] exp (cid:18) − r b r ex (cid:19) ζ = E r ζU [ ξ ] = 2 a Z L − L exp (cid:18) − z b r ex (cid:19) | ξ ( z ) | dz We can use the same basis for ξ ( z ) as in (12). When there’s no electric field, this basis is itself the solution to the firstequation in (23). Under electric field, the coefficients are found as usual by numerical matrix diagonalization.Now let us find the average of interaction potential U [ ξ ] explicitly: U [ ξ ] = √ π aν N z X n z ,n ′ z =1 C n z C n ′ z (cid:2) V n z − n ′ z ( ν ) − V n z + n ′ z ( ν ) (cid:3) (24)9 PREPRINT - J
ANUARY
18, 2021We will use two possible ways to solve the second equation in (23), which can be called the electron-hole interactionequation. The first is finite difference method, as presented in 2.3. It’s straightforward to implement in this case. Theabsorption spectrum then has the form: α ( ~ ω ) = 43 α n b ~ m E p L · ℑ " N z X λ z =1 N r +1 X λ r =1 /r ex · Q Dλ z λ r · ~ w ( E λ z + E λ r ) [( E λ z + E λ r ) − ( ~ w ) ] (25)Where: Q Dλ z λ r = 4 r ex (∆ r ) | φ ( λ z λ r )0 | (26)We need to numerically solve the interaction equation for each state λ z . The second way to solve (23) is variational method, which was one of the first employed for the case of Coulombexcitons [2, 28] and is still being used to this day [17]. In the present case, the most simple choice of ansatz forthe ground state is a Gaussian function (this choice is also supported by the shape of numerically found ground statewavefunction in Figure 2(c)). ζ ( r ) = r π ce − c r (27) Z + ∞ πrdr | ζ ( r ) | = 1 Let’s find the energy in this state: E ( c ) = h ζ | H r | ζ i = ~ c m eh − RyU b r ex c b r ex c (28)Minimization of E ( c ) leads to a simple quadratic equation with exact solution: c = 1 √ br ex q b p U − (29) E ( c ) = − Ryb (cid:16) b p U − (cid:17) (30)Then the variational absorption spectrum will have the following form: α ( ~ ω ) = 43 α n b ~ m E p L · ℑ " N z X λ z =1 /r ex · Q λ z · ~ w ( E λ z + E λ z ) [( E λ z + E λ z ) − ( ~ w ) ] (31)Where: Q λ z = 1 b (cid:16) b p U [ ξ λ z ] − (cid:17) N z X n z =1 C ( λ z ) n z sin πn z ! (32)When F = 0 the expressions simplify to: 10 PREPRINT - J
ANUARY
18, 2021 U [ ξ λ z ] = √ π aν [ V ( ν ) − V λ z ( ν )] Q λ z = 1 b (cid:16) b p U [ ξ λ z ] − (cid:17) (33) E λ z + E λ z = E g + ~ π λ z m eh L − Ryb (cid:16) b p U [ ξ λ z ] − (cid:17) Note that for λ z = 1 we have the following relations: b = s − E RyQ (34) a = b ( b Q + 1) √ π ν [ V ( ν ) − V ( ν )] Which allow us to fit the interaction potential parameters to any fixed values of E , Q .This ansatz, however, doesn’t accurately describe the strong confinement case, as can be seen from Figure 2(d), wherethe ground state wavefunction doesn’t have Gaussian shape. Thus, we will use a more advanced choice: ζ ( r ) = r π (cid:18) p + 4 sp + q + s q (cid:19) − / (cid:16) e − p r + se − q r (cid:17) (35) E ( p, q, s ) = 2 Ryr ex (cid:18) p + 4 sp + q + s q (cid:19) − ·· (cid:20) s + 8 sp q ( p + q ) − b U (cid:18) b r ex p + 1 + 2 sb r ex ( p + q ) + 1 + s b r ex q + 1 (cid:19)(cid:21) (36)Instead of explicitly calculating partial derivatives of E ( p, q, s ) and solving the resulting system of non-linear alge-braic equations, we will use the numerical minimization routine implemented in Wolfram Mathematica. Q ( p, q, s ) = 2 r ex (1 + s ) (cid:18) p + 4 sp + q + s q (cid:19) − (37)In Section 5 we will denote the results obtained from ansatz (27) as "variational 1" and from (35) as "variational 2". In this case instead of the last two equations in (23) we have the following: − ~ m eh r ∂∂r (cid:18) r ∂ζ∂r (cid:19) − RyU [ ξ ]( r ) ζ = E r ζ (38) U [ ξ ]( r ) = 2 r ex Z L − L | ξ ( z ) | √ r + z dz In this case the exact explicit expression for the averaged potential in the chosen basis is not straightforward to derive,however there’s a way to avoid oscillatory integrals, which we will present here. U [ ξ ]( r ) = 2 r ex L N z X n z ,n ′ z =1 C n z C n ′ z (cid:2) W n z − n ′ z ( µ ) − W n z + n ′ z ( µ ) (cid:3) (39)11 PREPRINT - J
ANUARY
18, 2021Where µ = r/L and: W m ( µ ) = cos πm Z cos (cid:0) πm t (cid:1)p t + µ dt (40)For odd m it follows that W m ( µ ) = 0 , and for even m we have: W m ( µ ) = K (cid:18) π | m | µ (cid:19) − cos πm Z + ∞ uJ ( µu ) e − u u + π m / du (41)For the special case m = 0 : W ( µ ) = asinh µ (42)For finite difference approximation at j = 0 we again use the averaging with weight function πr from r = 0 to r = ∆ r/ , which leads to the following expressions: h W m (cid:16) rL (cid:17)i j =0 = 4 L ∆ r · · " π | m | (cid:18) Lπ | m | ∆ r − K (cid:18) π | m | ∆ r L (cid:19)(cid:19) − cos πm Z + ∞ J ( ∆ r L u ) e − u u + π m / du (43) h W (cid:16) rL (cid:17)i j =0 = asinh (cid:18) L ∆ r (cid:19) + 11 + p r ) / (4 L ) (44)Where J , J , K , K are Bessel functions and modified Bessel functions. Integration w.r.t. u was performed usingGauss-Laguerre quadrature, which allowed to achieve the accuracy of at least digits using quadrature points.We did not attempt variational approximation for this case, as it offers no advantage over finite difference method, dueto the complexity of integrals involved for any choice of variational ansatz. The results of previous sections were used to fit the parameters a, b of the Gaussian potentials to the Coulomb case,but only for
L < r ex , as previous calculations showed that separation of variables approach only works in this range.To fit the parameters, we first calculated the L dependence of the Coulomb absorption spectrum, as described in theprevious section. In particular, we were interested in the ground state binding energy E and oscillator strength param-eter Q . Then we used (34) to obtain the initial fit of the parameters a, b . However, since the first variational ansatzdoesn’t give very accurate results, we corrected them for a better fit, comparing them with the pure 2D parametersfrom Table 3. Finally, we utilized the second variational ansatz to check the results and further correct the parameters.Our goal is to obtain a good agreement between the energies and oscillator strengths (absorption peak magnitudes) inthe whole range of quantum well widths, but fitting them using approximate solution for the Coulomb potential leadsto incorrect results for L > r ex . In the end, we had to use a piecewise function for a ( L/r ex ) and b ( L/r ex ) which isdescribed below: a ( y ) = a D + P m =1 A m y m / ( v m + y m ) , y ≤ / A y + A y + A , / < y ≤ a D − A / ( v + y ) , y > . (45) b ( y ) = b D + P m =1 B m y m / ( v m + y m ) , y ≤ / B y + B y + B , / < y ≤ b D − B / ( v + y ) , y > . (46)12 PREPRINT - J
ANUARY
18, 2021Table 4: Parameters for (45) and (46) A . B . A − . B − . A . B . A − . B − . A . B . A . B . A . B . Where a D , a D , b D , b D are taken from Table 3, y = L/r ex , v = 3 / , the parameters A − , B − were obtained byleast squares fit to the Coulomb results for E , Q at eF = 0 , while A − , B − were used to obtain continuous andmonotonic functions, increasing from a D ( b D ) at L = 0 to a D ( b D ) at L → ∞ . The values for all the parametersare provided in Table 4. Note that when presenting the calculation results, we subtract the following expression from the energy (which in thefigures’ labels we denote simply by E g ): ˜ E g = E g + ~ π m eh L This allows us to discard the strong L dependence of the QW ground state and focus instead on the binding energychange, when comparing the results for QWs of different widths.The implementation details are as follows. All the FD calculations were performed using R language (base packageonly) on an Intel ® Core™ i5-3210M CPU @ 2.50GHz × z direction, as we have to solve the interaction equation for each of them. Overall, depending on the above conditions,and for grid sizes of − the whole spectrum evaluation procedure took from around to minutes.As for the spectral method, the calculations were performed using Intel® CoreTM i7-4930MX CPU @ 3.00GHzlaptop with 32 Gb RAM. In this case the most time consuming task was evaluation of matrix elements. We haveutilized the symmetry, which required us to only evaluate the main diagonal and the upper triangular part of theHamiltonian matrix, which for two variables is equivalent to the condition N r ( n z −
1) + n r ≥ N r ( n ′ z −
1) + n ′ r .Wolfram Mathematica’s Parallelize function was used to speed up the process. When computing the spectrum itself,the summation was restricted only to the states with the energy between ~ ω min and ~ ω max . For N z = 30 and N r = 150 the whole procedure took about an hour. For N z = 50 and N r = 150 approximately . hours, and for N z = 50 and N r = 180 up to − hours.However, it should be noted that computation of spectra for different values of electric field doesn’t require evaluatingthe Hamiltonian matrix more than once, as the electric field strength is just a scalar factor in the equation. Whichis why in this case, after initial matrix evaluation had been completed, each spectrum took less than minutes tocompute, making this method very efficient for modeling electro-optic phenomena. We have calculated the absorption spectra without electric field for QW widths ranging from effectively (the 2Dstrong confinement regime) to L = 5 r ex , which, as we can see from Figure 3(c,e,f), gives almost identical results tothe bulk 3D case.We have used the following parameters: R = 500 nm for all the methods, for the spectral method, N z = 30 and N r = 150 ; for the finite difference method (within the separation of variables approach) ∆ r = 0 . nm.As can be seen from Figure 3(a-e), there’s excellent agreement for L ≤ r ex between all the calculation methods for theGaussian interaction potential, as well as the first absorption peak for Coulomb excitons (calculated with separation ofvariables approach together with the finite difference method for the interaction equation).13 PREPRINT - J
ANUARY
18, 2021 −5 0 5 10 . . . . . . Photon energy − E g , Ry A b s o r p t i on , m m - L = 0.5 r ex Full spectralFinite differenceCoulomb FD −5 0 5 10 . . . . . Photon energy − E g , Ry A b s o r p t i on , m m - L = 1 r ex Full spectralFinite differenceCoulomb FD (a) (b) − − − − QW width, r ex B i nd i ng ene r g y , R y Full SpectralFinite differenceVariational 2Variational 1Coulomb FD
QW width, r ex Q D Full SpectralFinite differenceVariational 2Variational 1Coulomb FD (c) (d) . . . . . QW width, r ex Q D - Full SpectralFinite differenceVariational 2Variational 1Coulomb FD −5 0 5 10 . . . . . . Photon energy − E g , Ry A b s o r p t i on , m m - L = 1 r ex L = 1.5 r ex L = 5 r ex Bulk (3D) (e) (f)Figure 3: Absorption spectra, calculated with different methods for (a) L = 0 . r ex ≈ nm and (b) L = 1 r ex ≈ nm. (c) Exciton binding energy. (d) 2D dimensionless oscillator strength. (e) Reciprocal 3D dimensionless oscilla-tor strength, proportional to peak magnitude. (f) Absorption spectra, calculated with the spectral method for differentQWs compared with the bulk (3D) case, calculated with FD.14 PREPRINT - J
ANUARY
18, 2021Simple variational ansatz with a single parameter (we denote results "Variational 1") gives a reasonably good agree-ment for both the energy and magnitude of the first peak, but a more elaborate three-parameter choice ("Variational 2")provides a much better accuracy, practically indistinguishable from numerically exact finite difference calculations.When it comes to
L > r ex , the spectral solution of the full excitonic equation with no separation of variables assump-tion (we denote it "Full spectral") gives results qualitatively different from the approximations.It should be noted, that piecewise fit of the Gaussian potential parameters to the Coulomb case does not change thequalitative results, only provides the necessary quantitative agreement in the range of applicability of the separationof variables approximation. Keeping the parameters a, b constant, for example, fitting them to the 3D Coulomb caseonly, still gives the same qualitative outcome when it comes to the L − dependence, namely, monotonic decrease ofthe binding energy (Figure 3c) and monotonic decrease of the peak magnitude (Figure 3e) from pure 2D to pure3D values. The separation of variables approach only works for L ≤ r ex for either constant or variable interactionpotential parameters.Dimensionless quantities proportional to oscillator strength Q D and Q D provide a useful comparison between QWsof different width. We can see that for L > r ex the full spectral solution gives a physically correct picture where Q D experiences linear growth (Figure 3d), while Q D approaches a constant value for free bulk excitons (Figure 3e). Figure 4(a) shows the field dependence of the absorption spectrum for a QW with L = 0 . r ex which is approximatelyequal to nm for the material considered in this work. The color describes the magnitude of the absorption, darkbands are the absorption peaks. Note, that we used dimensionless excitonic units for the field, while the range inabsolute units is F = [0 , kV/cm.We can observe that both magnitude and position of the peaks change with the field, and while the ground statetransition is suppressed by the field, the previously forbidden second transition is, to the contrary, enhanced.Due to the small width, all calculation methods are more or less in agreement, including the results for the Coulombexciton (obtained, as mentioned above, using the separation of variables approach). There’s a better agreement betweenpeak positions than peak magnitudes (Figures 4(b,c)).The main difference between the spectra for the Coulomb and Gaussian excitons is the Sommerfeld enhancement,which is much stronger for the former (Figures 4(d-f)).Figures 5(a,b) for wider QWs provide more information about the effects of the electric field on the absorption. If wenumber the QW states in order of increasing energy as λ = 1 , , , . . . , then in the absence of the electric field ( F = 0 )every odd numbered state corresponds to an allowed transition, while every even numbered state - to a forbidden one.This is because in our initial problem statement we only considered the relative motion of the electron-hole pair. Whichmeans that when z = 0 both electron and hole have the same absolute coordinate ( z e = z h ). If the relative motionwavefunction has a node at this point (which is true for all the even numbered states), then absorption is forbidden.Application of the electric field displaces both particles in the opposite directions, which leads to a deformation ofthe relative motion wavefunction. Depending on the number of nodes, the optical transitions are either suppressed orenhanced at small fields. For multi-node states (with λ > ) we observe non-monotonic behavior of the absorptionpeak magnitude. This is caused by the fact that with increasing field strength several maxima of electron and holerespective wavefunctions are displaced in the opposite directions and periodically overlap.At very large fields, all the absorption peaks are eventually displaced in the direction of smaller energies (largerwavelengths) and their magnitude decreases until they are suppressed completely.In addition, we have plotted field dependence of the integrated absorption spectra for three different QWs (Figure 5(c)).There’s a so called "sum rule" [3] which states that this quantity should stay constant. In our case, it’s not fulfilledexactly, but the change is only about − even at relatively large fields. The discrepancy can be easily explained bythe finite range of photon energies ~ ω we considered. Which means some exciton peaks, which remained out of rangein our calculations, should affect the integrated absorption and provide a better agreement with the sum rule.Finally, we have calculated the electroabsorption effect (the absolute change of absorption coefficient due to appliedfield) at ~ ω = 0 . eV (Figure 5(d)). This energy is below the bulk bandgap when no field is applied. In this calculation,unlike all the previous ones, we have accounted for the full quantum size effect for each quantum well, by adding thefollowing quantity to the transitions energies: 15 PREPRINT - J
ANUARY
18, 2021 − Electric field, Ry/r ex P ho t on ene r g y − E g , R y − − Electric field, Ry/r ex P ea k ene r g y − E g , R y Full spectralFinite differenceCoulomb FD (a) (b) . . . . Electric field, Ry/r ex P ea k m agn i t ude , m m - Full spectralFinite differenceCoulomb FD Peak 1Peak 3Peak 2 −20 0 20 40 . . . Photon energy − E g , Ry A b s o r p t i on , m m - eF = 1 kV/cmFull spectralFinite differenceCoulomb FD (c) (d) −20 0 20 40 . . . . . . Photon energy − E g , Ry A b s o r p t i on , m m - eF = 40 kV/cmFull spectralFinite differenceCoulomb FD −20 0 20 40 . . . . Photon energy − E g , Ry A b s o r p t i on , m m - eF = 80 kV/cmFull spectralFinite differenceCoulomb FD (e) (f)Figure 4: L = 0 . r ex ≈ nm: (a) Absorption spectrum dependence on the electric field. Color describes magnitudeof the absorption coefficients. (b) Positions of absorption peaks. (c) Magnitudes of absorption peaks. Absorptionspectra for different electric field strengths: (d) F = 1 kV/cm, (e) F = 40 kV/cm, (f) F = 80 kV/cm.16 PREPRINT - J
ANUARY
18, 2021 − − Electric field, Ry/r ex P ho t on ene r g y − E g , R y − − Electric field, Ry/r ex P ho t on ene r g y − E g , R y (a) (b) . . . . . . Electric field, Ry/r ex I n t eg r a t ed ab s o r p t i on , e V / m m L = 0.75 r ex L = 1 r ex L = 1.75 r ex − − − − Electric field, Ry/r ex A b s o r p t i on c hange , m m - L = 0.75 r ex L = 1 r ex L = 1.75 r ex (c) (d)Figure 5: Absorption spectra dependence on the electric field. Color describes magnitude. (a) L = 1 r ex ≈ nm, (b) L = 1 . r ex ≈ nm. (c) Dependence of integrated absorption spectrum on the electric field. (d) Electroabsorptioneffect at ~ ω = 0 . eV, log scale. ∆ E = π ~ m eh L This ensured that the relative positions of the absorption peaks are correct. The results (presented on a logarithmicscale) indicate initial exponential growth below the bandgap, which is stronger for larger QWs. However, due to thedisplacement of the ground state absorption peak, there’s a critical field strength at which the absorption increaseis replaced by decrease. This critical value is greater for smaller QWs, which makes them more convenient forelectroabsorption modulator applications.
In conclusion, we have numerically investigated excitonic electroabsorption in a quantum well, using a modifiedinteraction potential, namely a Gaussian function. The purpose of this work was to utilize the favorable properties ofthe latter (namely, the lack of singularity at r = 0 and finite range) to simplify the numerical solution of Schrödingerequation for confined excitons under external electric field.We have developed a finite difference scheme for the radial Schrödinger equation with either axial or central symmetry,using weighted averaging for Coulomb potential at r = 0 . The FD calculations are in excellent agreement with knownanalytical expressions for the absorption spectrum. 17 PREPRINT - J
ANUARY
18, 2021We have managed to fit the free parameters of the Gaussian potential (magnitude and effective range) to the usualCoulomb case (for the first absorption peak) for 3D and 2D excitons. We then used the separation of variables approachto obtain an approximate solution for both cases applicable to thin quantum wells (
L < r ex ). This allowed us to achievea good agreement between the Coulomb and Gaussian excitons for arbitrary QW widths. Thus, we can expect thatelectroabsorption modeling conducted in the present work can be used to analyze the behavior of actual excitons notonly qualitatively, but in some cases quantitatively as well, especially when it comes to the absorption edge.As for the numerical methods used in this work, the most effective way to solve the full exciton Schrödinger equationin 2 variables proved to be the spectral expansion method. We have managed to find either exact or approximateexpressions for all the matrix elements and find the absorption spectra for a wide range of QW widths. The results forvery wide wells agree with bulk (3D) calculations performed by finite differences.As a promising way to decrease the calculation time, we have also used another approach - the separation of variables.In this case, the interaction potential, which can’t be separated exactly, is replaced by its average over the extra vari-able(s). To solve the resulting equation, one can either use a single parameter variational ansatz, a multi-parametervariational ansatz, or a finite difference scheme. We found that the first approach, while being the most simple and lead-ing in this case to a semi-analytical solution, doesn’t give the best agreement with the original two-variable equation.Meanwhile, a more advanced variational function or numerical finite difference solution provide an almost perfectagreement for L < r ex .Regarding the electric field dependence of the absorption spectrum, we find that for small fields the results agreewell with the perturbation theory treatment, i.e. the ground state energy decreases, while all the rest of the energylevels increase. The displacement of the electron and hole wavefunctions to the opposite sides of the quantum wellresults in either decreasing or increasing magnitude of the absorption peaks, which is related to the symmetry andrelative position of the wavefunction extrema. Which is why for excited states with several such extrema we observenon-monotone field dependence of both the magnitude and transition energy.As for electroabsorption below the edge, the magnitude of the effect varies with QW width, and for certain fields wefind that quantum wells of intermediate thickness perform better than those that are either too thin or too thick.The main conclusions of this work are in good agreement with previously published results on Coulomb excitons,which means that the numerically stable Gaussian interaction potential may be useful for modeling more complexphenomena as well, especially with many interacting particles or fast-changing fields. Acknowledgments
The research was carried out within the state assignment of FASO of Russia N o References [1] D.A.B. Miller, D. S. Chemla, T. C. Damen et al., Band-Edge Electroabsorption in Quan-tum Well Structures: The Quantum-Confined Stark Effect, Phys. Rev. Lett., 53 (1984), 2173. https://doi.org/10.1103/PhysRevLett.53.2173 [2] D.A.B. Miller, D.S. Chemla, T.C. Damen et al., Electric Field Dependence of Optical Absorp-tion near the Band Gap of Quantum-Well Structures, Phys. Rev. B, 32, 2 (1985), 1043–60. https://doi.org/10.1103/PhysRevB.32.1043 [3] D. Miller, J. Weiner, D. Chemla et al., Electric-field dependence of linear optical properties in quantum wellstructures: Waveguide electroabsorption and sum rules, IEEE J. Quantum Electronics, 22, 9 (1986), 1816-1830. https://doi.org/10.1109/JQE.1986.1073167 [4] P.C. Klipstein and N. Apsley, A theory for the electroreflectance spectra of quantum well structures. J. Phys. C:Solid State Phys., 19 (1986), 6461. https://doi.org/10.1088/0022-3719/19/32/020 [5] Kan, Y., H. Nagai, M. Yamanishi et al., Field Effects on the Refractive Index and Absorption Coefficient inAlGaAs Quantum Well Structures and Their Feasibility for Electrooptic Device Applications, IEEE J. Quant.Electronics, 23, 12 (1987), 2167–80. https://doi.org/10.1109/JQE.1987.1073288 [6] J.S. Weiner, D.A.B. Miller and D.S. Chemla, Quadratic electro-optic effect due to the quantum-confined Starkeffect in quantum wells, Appl. Phys. Lett., 50 (1987), 842-844. https://doi.org/10.1063/1.98008 [7] H. Mohseni, H. An, Z. A. Shellenbarger, M. H. Kwakernaak, and J. H. Abeles, Enhanced Electro-Optic Effect in GaInAsP-InP Three-Step Quantum Wells, Appl. Phys. Lett., 84, 11 (2004), 1823–25. https://doi.org/10.1063/1.1682699 PREPRINT - J
ANUARY
18, 2021[8] T. Arakawa, T. Toya, M. Ushigome et al., InGaAs/InAlAs five-layer asymmetric coupled quan-tum well exhibiting giant electrorefractive index change, Jap. J. Appl. Phys., 50 (2011), 032204. https://doi.org/10.1143/JJAP.50.032204 [9] D.P. Dave, Enhanced Refractive Index Change in Asymmetrical Quantum Wells with an Applied Electric Field,J. Appl. Phys., 74, 11 (1993), 6872–75. https://doi.org/10.1063/1.355089 [10] H. Feng, J.P. Pang, M. Sugiyama et al., Field-Induced Optical Effect in a Five-Step Asymmetric Cou-pled Quantum Well with Modified Potential, IEEE J. Quant. Electronics, 34, 7, (1998), 1197-1208. https://doi.org/10.1109/3.687863 [11] S.L. Chuang, S. Schmitt-Rink, D.A. Miller, and D.S. Chemla, Exciton Green s-function approach to op-tical absorption in a quantum well with an applied electric field, Phys. Rev. B, 43, 2 (1991), 1500-1509. https://doi.org/10.1103/PhysRevB.43.1500 [12] S. Glutsch, D.S. Chemla, and F. Bechstedt., Numerical Calculation of the Optical Absorption in SemiconductorQuantum Structures. Phys. Rev. B, 54, 16 (1996), 11592. https://doi.org/10.1103/PhysRevB.54.11592 [13] A. Ahland, D. Schulz, E. Voges et al., Efficient Modeling of the Optical Properties of MQW Modula-tors on InGaAsP with Absorption Edge Merging, IEEE J. Quant. Electronics, 34, 9 (1998), 1597–1603. https://doi.org/10.1109/3.709576 [14] J. Hader, J.V. Moloney, and S.W. Koch, Microscopic Theory of Gain, Absorption, and Refractive Index in Semi-conductor Laser Materials—Influence of Conduction-Band Nonparabolicity and Coulomb-Induced IntersubbandCoupling, IEEE J. Quant. Electronics, 35, 12 (1999), 1878-1886. https://doi.org/10.1109/3.806602 [15] A.V. Maslov and D.S. Citrin., Numerical Calculation of the Terahertz Field-Induced Changes in theOptical Absorption in Quantum Wells, IEEE J. Sel. Topics Quant. Electronics, 8, 3 (2002), 457-463. https://doi.org/10.1109/JSTQE.2002.1016348 [16] K. Sivalertporn, L. Mouchliadis, A. L. Ivanov et al., Direct and indirect excitons in semicon-ductor coupled quantum wells in an applied electric field, Phys. Rev. B, 85 (2012), 045207. https://doi.org/10.1103/PhysRevB.85.045207 [17] V.V. Zolotarev, I.S. Shashkin, V.S. Golovin et al., Enhancement of the refractive index modulationin a modulator based on GaAs/AlGaAs quantum wells, Semicond. Sci. Technol., 34 (2019), 095005. https://doi.org/10.1088/1361-6641/ab3131 [18] F. Duque-Gomez and J.E. Sipe, The Franz–Keldysh Effect Revisited: Electroabsorption In-cluding Interband Coupling and Excitonic Effects, J. Phys. Chem. Solids, 76 (2015), 138–52. https://doi.org/10.1016/j.jpcs.2014.07.023 [19] Yue Hu, Curtis R. Menyuk, Meredith N. Hutchinson et al. Impact of the Coulomb interactionon the Franz–Keldysh effect in high-current photodetectors. Opt. Lett., 41, 3 (2016), 456-459. https://doi.org/10.1364/OL.41.000456 [20] S. Schmitt-Rink, D.S. Chemla, W.H. Knox, and D.A.B. Miller. How fast is excitonic electroabsorption? Opt.Lett., 15, 1 (1990), 60-62. https://doi.org/10.1364/OL.15.000060 [21] R Core Team (2020). R: A language and environment for statistical computing. R Foundation for StatisticalComputing, Vienna, Austria. URL [22] Wolfram Research, Inc., Mathematica, Version 11.2, Champaign, IL (2017).[23] C. Tanguy. Optical Dispersion by Wannier Excitons, Phys. Rev. Lett., 75 (1995), 4090. https://doi.org/10.1103/PhysRevLett.75.4090 [24] T.Garm Pedersen. Excitons in Bulk and Two-Dimensional Semiconductors. Lectures for PhD course "Optics atthe Nanoscale", Aalborg University. http://homes.nano.aau.dk/tgp/ [25] C. Tanguy, Complex dielectric constant of two-dimensional Wannier excitons, Solid State Comm., 98 (1996), 65. https://doi.org/10.1016/0038-1098(95)00750-4 [26] I.S. Gradshteyn, I.M. Ryzhik, Y.V. Geronimus, and M.Y. Tseytlin, Table of Integrals, Series, and Products. (inRussian) (4 ed.). Moscow: Gosudarstvennoe Izdatel’stvo Fiziko-Matematicheskoy Literatury (Fizmatgiz), 1963.[27] S.L. Chuang, Physics of Photonics Devices, Wiley-Interscience, New York, 2009.[28] T. Ishikawa, S. Nishimura and K. Tada, Quantum-Confined Stark Effect in a Parabolic-Potential Quantum Well,Jap. J. Appl. Phys., 29, 8 (1990), 1466–73. https://doi.org/10.1143/JJAP.29.1466https://doi.org/10.1143/JJAP.29.1466