Yukawa model of screening in low-dimensional excitons: diagonalization, perturbation, variation, and resummation analysis

The Yukawa interaction describes Coulomb screening in several physical systems. We investigate bound quantum states of the Yukawa potential in low-dimensional structures as a model of e.g. excitons in semiconductor nanostructures. Diagonalization, perturbation, variation, and resummation methods are all applied to the problem and their accuracy is compared. For moderate positive screening, all methods are found to be applicable and, in particular, the variational approach is highly accurate except near the threshold for bound states. In contrast, only hypergeometric resummation captures the correct behaviour in the negative (antiscreened) regime. We also determine the dimensional dependence of the critical screening, above which bound states cease to exist. As an application, the critical screening is used to compute critical doping levels for the existence of bound excitons in quantum wells.


Introduction
The Yukawa (or screened Coulomb or Debye-Hückel) potential V r e r r = -k -( ) / has played a prominent role in physics since its original introduction as a model of nucleon interactions [1]. The coefficient in the exponent κ (henceforth called the 'screening') truncates the range of the bare Coulomb potential. This situation arises whenever screening by mobile charges modifies the Coulomb interaction. Important examples include plasma interactions, electrolytes, and electron-hole pairs, i.e. excitons, in semiconductors [2][3][4][5][6][7][8][9][10][11]. Hence, although the Yukawa model as introduced in a very different context, it has been widely applied in condensed matter physics. It plays a significant role in semiconductor physics [4]. In particular, following the pioneering work by Gay [5], exciton screening due to free carriers from doping or intense excitation has frequently been treated using Yukawa potentials [5][6][7][8][9][10][11] or the closely related Hulthén model [6,12]. Fundamentally, exciton screening is a dynamical process but a static approximation is frequently adopted due to the inherent simplicity. Nonetheless, excellent agreement with density-dependent exciton absorption in GaAs [6] and ultrafast screening dynamics of ZnO excitons has been found [10]. Moreover, the static Yukawa model has been applied to screening of excitons in one-and two-dimensional semiconductors [7,8,11] as well as bulk materials [5,6,9,10].
Bound states of quantum particles interacting via the Yukawa potential have been studied for many decades and their properties in three dimensions are well understood. Several studies applying finite difference equation [13], numerical diagonalization [14], perturbation theory [15][16][17], 1/N expansion [18], variational approaches [19,20], Padé resummation [15,21], or pseudospectral methods [22] have been reported. A particularly interesting feature of the three-dimensional Yukawa model is that bound states cease to exist when the screening exceeds a critical value κ>κ crit ≈1.19061 in atomic units [17,19,21]. Similarly, a few studies have considered the mathematical properties of the Yukawa model in two dimensions [18,23]. In this case, no critical screening exists and bound states are always found. Comparatively few analytical results are known, however, and a general analysis of the low-dimensional Yukawa model is lacking. For excitonic models, such low-dimensional geometries are highly relevant. For instance, electron-hole pairs in narrow quantum wells are approximately two-dimensional. Moreover, several such geometries are best described using fractional dimensions. Thus, quantum wells have effective dimensions 2<D<3 [24] and carbon nanotubes have dimension D≈1.7 [25]. Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Consequently, it is relevant to study electron-hole pairs bound by screened Coulomb potentials in fractional dimensions. In particular, the dependence of κ crit on D is an open problem. It should be noted that an accurate model of exciton screening in low-dimensional structures such as nanotubes [26] and monolayer semiconductors [27] is itself a complicated problem. As mentioned above, dynamical effects are of importance and screening affects a range of material properties other than excitons [28]. One such property is band gap renormalization. In the present work, however, we only consider the exciton binding energy, i.e. energy relative to the band gap, and consequently band gap renormalization plays no role for the results. Hence, the statically screened Yukawa interaction model can only be seen as a simplistic approximation to the full problem. The virtue of the model, however, lies precisely in the simplicity, which allows for analytical or semi-analytical results as well as transparent physics.
In the Stark problem, (low-dimensional) excitons exposed to a constant electrostatic field shift their energy and may dissociate with a characteristic field-dependent decay rate Γ. The computation of Γ is a complicated non-perturbative problem [29] but accurate numerical values can be extracted from the complex-scaling approach [30]. However, a computationally much simpler alternative is resummation based on a few terms of a perturbation series. Recently, a highly efficient hypergeometric resummation method has been proposed [31][32][33] and applied to low-dimensional Coulomb-bound systems. Excellent agreement with complex scaling has been demonstrated for Stark-type cases but evidence from a wider class of problems is needed to establish the broad applicability of the method. In fact, the Yukawa model shares several properties of the Stark problem. Thus, if negative screening κ<0, i.e. antiscreening cf figure 1, is allowed, bound states become resonances that will decay in time. It follows that the Yukawa potential can also facilitate a testbed for hypergeometric resummation. Hence, similarly to the screened model, the antiscreened version is a highly simplistic model of e.g. exciton ionization, whose virtue lies in its simplicity and applicability as a testbed for calculations of resonances.
In this paper, we consider the low-dimensional Yukawa potential for geometries with dimension 1<D3 We will compare four approaches: (1) numerical diagonalization, (2) perturbation expansion, (3) variational estimates, and (4) Padé and hypergeometric resummation. Specifically, the screening κ dependence of the D-dimensional ground state energy will be investigated. Also, the dependence of κ crit on dimension will be obtained from both diagonalization, variational, and resummation approaches. Finally, we compare numerical complex scaling results for the unstable antiscreened case with hypergeometric resummation. In the next section, the theoretical framework behind all four approaches is established. Subsequently, in section 3, a numerical comparison of all methods is presented. Then, we apply the model to bound excitons in quantum wells in section 4. A summary and conclusions are offered in the last section.

Theoretical framework
Our starting point is the Schrödinger equation in natural units using the reduced mass μ of the two-particle system as the unit of mass: e 4 1 r 0  pe e m = = = = with ε r the relative dielectric constant of the ambient medium. Hence, V r with L the angular momentum operator [24], and the Yukawa potential given by We have provided the Taylor series of the potential for later reference. Note that the radial part of the Laplacian interpolates linearly between the known 1D, 2D, and 3D limits. Moreover, the radial integration measure is r dr D 1 - [24]. Throughout, we restrict the analysis to rotationally-symmetric states with vanishing angular momentum, which include the ground state. In the following, we describe four distinct methodologies applied to compute the ground state.

Numerical diagonalization
The eigenvalue problem equation (1) can be solved to arbitrary precision using an expansion in a complete basis. Since we focus on the ground state, we need only include states with vanishing angular momentum, i.e. spherically symmetric ones. To this end, we apply a radial Gaussian basis i ar exp Here, U is the confluent hypergeometric function. In the numerical results below, exponents are selected according to a 3 10 2

Perturbation theory
For small values of κ, a low-order expansion of the ground state energy may be sufficient. This expansion is partially captured by first order perturbation theory. Thus, we use the hydrogen ground state for the unscreened potential [24] r D D e This state is radially normalized r r dr 1.
This contribution to the full expansion of the energy is exact to third order but the second order perturbation correction leads to contributions of order κ 4 and higher. Looking at the Taylor series in equation (2), the relevant perturbation is V r 2.
2 k D = -/ We apply Dalgarno-Lewis perturbation theory [35] to compute the exact second order contribution. Hence, writing 0 1 j j j = + + ¼with the subscript indicating order in κ 2 , we proceed by collecting all terms of order κ 2 in the exact eigenvalue problem leading to Here, the right-hand side represents the κ 2 term in E . 1 k ( ) The solution is readily found to be D r 1 . . 7 In total, collecting zeroth, first, and second order corrections, The first of these agrees with the perturbation series in [15].

Variational ansatz
The variational approach has previously been applied to the three-dimensional problem with success [19]. We therefore attempt to adapt it to the D-dimensional case. We choose the same ansatz as the one applied in [19], i.e.
with a and b variational parameters. This form resembles the one applied successfully to the two-dimensional exciton problem in [36]. Note that the hydrogen solution equation (4) is contained in the ansatz since From the variational ansatz, we estimate the energy using The parameters a and b are determined by a numerical binary search for the minimum of the variational estimate E.

Padé and hypergeometric resummation
Resummation is a technique designed to convert a finite perturbation series, such as equation (8) However, Padé resummation suffers from a number of weaknesses in problems characterized by branch cuts, such as the antiscreened regime of the Yukawa potential model. Fortunately, for these problems several alternative schemes exist. In [32], hypergeometric resummation was applied to the Stark effect in lowdimensional hydrogen. The choice in that case was a Gauss hypergeometric function F 2 1 of the form ¼ are undetermined parameters. The integer l>4 is a constant adjusted to the particular problem at hand, for which l=30 was chosen in [32]. In the present case, it is found that the limit l  ¥ is well suited, which leads to a confluent hypergeometric function The Taylor expansion of this form agrees by construction with the first two terms of equation (8) and h h 1 3 ¼ are fixed by matching to orders 2 4 k k ¼ in equation (8) In the next section, we numerically compare the four methodologies described in sections 2.1-2.4.

Numerical results
We now compare diagonalization, perturbation, variational, and resummation approaches to the Yukawa problem for various characteristic dimensions. We wish to establish the regime of validity of these different methods. For κ0, all four methods are, in principle, applicable. In contrast, for κ<0 the variational approach fails because the problem is not bounded from below. In figure 2, the 'exact' numerical diagonalization, perturbative, and variational approaches are compared for D=2.0, 2.5, and 3.0. In the unscreened Coulomb case κ=0, the energy in all methods is −2/(D−1) 2 . The fourth-order perturbation series deviates markedly from the exact result around κ≈0.5 for D=3.0 and fares slightly better in lower dimensions. Conversely, the variational estimates are barely distinguishable from the exact curves except for the largest values of κ. Thus, the variational approach apparently provides accurate results for a wide range of dimensions and screenings. However, close to the threshold κ=κ crit the relative accuracy is rather poor and, in particular, variational predictions for κ crit itself are not generally reliable, as described below. In contrast to perturbation and variational estimates, the diagonalization approach also provides accurate results for excited states. As an example, the inset in figure 2 shows diagonalization results for the 2s excited states. These are seen to follow a dimensional dependence similar to the ground state. For the 2s state, however, a critical screening exists even for D=2.
The numerical diagonalization approach also applies to the unstable (antiscreened) κ<0 case. However, it is necessary to adjust for the fact that real-valued eigenvalues become complex resonances if κ switches from positive to negative. The imaginary part E 2Im G =provides the decay rate corresponding to 'tunneling' into the ionized state. A highly efficient way of dealing with such resonances is through complex scaling [30,33,37], in which one rotates the Schrödinger equation into the complex plane r e r i  q   (we take θ=0.5 in practice).
This complexified problem can, of course, also be solved in the stable κ0 situation and one then finds Γ=0 and a real part equal to the usual eigenvalue.
In figure 3, complex scaling results are compared to resummation approximants for a range of characteristic dimensions. Both Padé (dotted lines) and hypergeometric (dashed lines) approximants are considered. Several important conclusions can be drawn from this comparison. The Padé approximant is generally excellent for κ0. In the unstable regime κ<0, however, it completely fails to capture the correct behavior. It has an identically vanishing imaginary part and, due to the rational form equation (12), a pole at a characteristic negative κ. For instance, a pole at κ≈−1.0 is observed in the three-dimensional case. In contrast, the hypergeometric approximant equation (14) provides a good approximation for both positive and negative κ. In particular, the correct decay rate behavior is observed in the unstable regime. Mathematically, this is a consequence of the branch cut z , 0 ] This property ensures a correct behavior for negative κ provided the real part agrees with the perturbation series near κ=0.
An important feature of the Yukawa potential is the existence of a critical screening κ crit , beyond which no bound states exist. A critical screening can be defined for both ground and excited bound states but here we exclusively consider the ground state. It is well-established that κ crit ≈1.19061 for D=3.0 [19]. Numerically, κ crit can be determined from the first zero-crossing of the energy. For accurate results using the diagonalization approach, it is necessary to add very small exponents to the basis set, however. The variational ansatz also provides a critical screening that, incidentally, coincides with a vanishing exponent b in equation (10) so that κ crit ≈1.19021K in excellent agreement with the exact value. In figure 4, we plot the dimensional dependence of the critical screening obtained from diagonalization, variation, and Padé resummation. Again, excellent agreement between the former two is noted near D=3.0 Significant deviations, however, are seen for lower dimensions. In fact, for sufficiently low dimensions, the Padé result is superior. None of the two approximate method correctly capture the divergence crit k  ¥ in the limit D→2, though. The variational ansatz at D=2.0 leads to 4 2 ln 2 1 10.355 whereas the Padé approximant incorrectly yields κ crit =16 and diverges at D≈1.893. For the intermediate dimension D=2.5, the 'exact' value is κ crit ≈2.976 while both approximations underestimate this.

Application to quantum wells
We will now apply the above theoretical framework to quantum wells with a finite thickness L that is allowed to vary between small and large values. We assume ideal confinement corresponding to hard-wall boundary conditions. Thus, by varying L we probe an electronic system, which effectively varies between two-and threedimensional limiting cases. Specifically, we consider excitons confined to such structures. This setting has previously been described using fractional dimensions with success [24]. By restricting the analysis to quantum wells, i.e. D 2 3 ,   we cover the range, for which a critical screening exists, cf figure 4. Note that, in the present natural units, the unit of length corresponds to the exciton Bohr radius. Hence, two-and threedimensional behavior is expected whenever L 1  and L 1,  respectively. Below, we investigate structures with L varying in the range L 0, 10 Î [ ]in order to cover the transition between the limiting cases. The confining direction is taken as the z-axis. Due to translational symmetry, exciton states are described by the in-plane radial coordinate ρ as well as electron and hole z-coordinates z e and z . ] We stress that this model is significantly more complicated than the fractional-dimensional one equation (1). In particular, equation (1) is effectively one-dimensional for s-states, whereas equation (18) implies a truly threedimensional problem. We consequently obtain the eigenstates of equation (18) with n, m integer and p b a set of suitably chosen exponents. Tractable expressions for the required overlap and Coulomb matrix elements are readily derived, whereas the kinetic energy leads to a rather complicated result. In the numerical applications, we take n, m between 1 and 6 and 2 p p 1 b =with integer p between 0 and 4. We then determine the critical screening by numerical identification of the first zero-crossing of the ground state binding energy.
We aim to compare the finite quantum well model to the much simpler fractional-dimensional one. To this end, we therefore need to find the effective quantum well dimension as a function of thickness L. There is no unique way of achieving this, but one useful and simple approach consists in matching the ground state energies of effective and full models [25]. In the case of carbon nanotubes, it is known that the full model [38] agrees with experiments [39]. Thus, the energy matching approach ensures accurate exciton energies in the effective model as well. In the D-dimensional model, the ground state exciton binding energy in the absence of screening is By identifying this expression with the numerically computed energy E 0 k= from the In three-dimensional materials, the usual result n E 6 is recovered. By setting κ equal to the critical value, a corresponding critical Fermi level for a given dimension D can be found from equation (23). We prefer to give the Fermi level rather than the critical density, as n D above is a mathematic construct (carriers per D-dimensional volume) while the Fermi level is a physical quantity. The result is shown in the bottom panel in figure 5  » exceeds the Fermi level by a sizable factor [10]. In the strict 2D limit L=0, the exciton ground state is stable, as mentioned above. However, for L≈2 the critical Fermi level E F ≈1 is already increased significantly over the bulk value. Hence, the intermediate case D=2.5 corresponding to L≈1.64 is stabilized to a large degree compared to the 3D case. It may be checked that the general plasma frequency n L 4 p D D 3 w p = exceeds the Fermi level in these cases, thus validating the static approximation. This analysis shows that the low-dimensional Yukawa model can be applied to provide concrete predictions for the properties of actual semiconductor nanostructures.

Summary
In summary, we have considered the Yukawa potential in fractional dimensional space in order to assess various approximations for the ground state energy and the critical screening. Both perturbation and variation methods have been attempted as well as analytical resummation based on Padé and hypergeometric approximants. For screening well below the critical value κ crit , a variational ansatz designed for the three-dimensional problem generally provides an excellent approximation for lower dimensions as well. However, this method fails to describe the ground state near κ crit and, in particular, fails to provide an accurate estimate of κ crit itself. In the negative or antiscreened regime, the bound states become unstable and eigenvalues are replaced by complex resonances. We compare numerically exact complex scaling results to the two resummation approximants and find that only hypergeometric resummation is able to capture the correct behaviour for the decay rate describing tunnelling ionization. Finally, by applying the model to bound excitons in quantum wells, we calculate the dimensional dependence of the critical doping level.