Morphology of wetting-layer states in a simple quantum-dot wetting-layer model

TL;DR

This paper investigates the complex morphology of wetting-layer states in quantum-dot models, revealing irregularities and chaos-like features that challenge simple plane wave approximations.

Contribution

It introduces a quantum chaos-based analysis to characterize wetting-layer states, highlighting their irregular and complex nature in a simple quantum-dot model.

Findings

Wetting-layer states exhibit irregular, complex morphology.

Standard plane wave approximation is insufficient.

Quantum chaos tools reveal spectral and localization properties.

Abstract

The excitation of semiconductor quantum dots often involves an attached wetting layer with delocalized single-particle energy eigenstates. These wetting-layer states are usually approximated by (orthogonalized) plane waves. We show that this approach is too crude. Even for a simple model based on the effective-mass approximation and containing one or a few lens-shaped quantum dots in a rectangular wetting layer, the wetting-layer states typically show a substantially irregular and complex morphology. To quantify this complexity we use concepts from the field of quantum chaos such as spectral analysis of energy levels, amplitude distributions, and localization of energy eigenstates.

Figures

Tables

Table 1. Table 1: Parameter sets of different QD-WL systems. a k subscript 𝑎 𝑘 a_{k} , e k subscript 𝑒 𝑘 e_{k} , b k subscript 𝑏 𝑘 b_{k} , L x subscript 𝐿 𝑥 L_{x} , and L y subscript 𝐿 𝑦 L_{y} are given in nm. d k subscript 𝑑 𝑘 d_{k} is given in nm -4 .

System	$N_{\text{dot}}$	$k$	$a_k$	$e_k$	$d_k$	$b_k$	$L_x$	$L_y$
$𝒜\left(t\right)$	$1$	$1$	$4.3t$	$6.3t$	$2.9$	$2.5$	$10t$	$10t\Phi$
$ℬ$	$4$	$1$	$11.575$	$18.698$	$2.9$	$2.5$	$34$	$15\Phi$
		$2$	$24.135$	$13.708$	$2.9$	$2.5$
		$3$	$25.827$	$9.446$	$2.9$	$2.5$
		$4$	$5.562$	$3.297$	$2.9$	$2.5$
$𝒞$	$1$	$1$	$4.3$	$6.3$	$1.1$	$2.5$	$10$	$10\Phi$

Equations

$$V(x,y)=\sum _{k=1}^{N_{\text{dot}}}{V}_k(x,y),$$

$$V_k(x,y)=\text{min}\left\{d_k\left[(x-a_k{)}^2+(y-e_k{)}^2-\frac{b_k^2}{4}\right],0\right\},$$

$$\hat{H}{\psi }_i(x,y)=E_i{\psi }_i(x,y)$$

$$\hat{H}=-\frac{{\hslash }^2}{2M}{\Delta }_{x,y}+V(x,y)$$

$${\psi }_i(x,y)=\sum _n{c}_n^{(i)}{\phi }_n(x,y).$$

$${\phi }_n(x,y)=⟨x,y|{\phi }_n⟩$$

$$⟨{\phi }_n|{\hat{H}}_0|{\phi }_m⟩=E_{n_x,n_y}{\delta }_{n_x,m_x}{\delta }_{n_y,m_y}$$

$$⟨{\phi }_n|V|{\phi }_m⟩=$$

$$\times \left[(x-a_k{)}^2+(y-e_k{)}^2-\frac{b_k^2}{4}\right],$$

$$N(E)=\underset{-\infty }{\overset{E}{\int }}d{E}^{\prime }\rho (E^{\prime })=\sum _i\Theta (E-E_i)$$

$$N(E)=\bar{N}(E)+N_{\text{fluc}}(E).$$

$$\bar{N}(E)\approx \frac{A}{4\pi }\frac{2ME}{{\hslash }^2}-\frac{L}{4\pi }\sqrt{\frac{2ME}{{\hslash }^2}}.$$

$$P_{\text{P}}(s)=e^{-s}.$$

$$P_{\text{W}}(s)=\frac{\pi }{2}s{e}^{-\pi s^2/4}.$$

$$P_{\text{SP}}(s)=4s{e}^{-2s}.$$

$$P_{\text{B}}(s,\beta )=(\beta +1)b{s}^{\beta }{e}^{-b{s}^{\beta +1}}$$

$$b={\left[\Gamma \left(\frac{\beta +2}{\beta +1}\right)\right]}^{\beta +1}$$

$$P_{\text{BR}}(s,\beta )$$

$$+\left(2\alpha \beta +\frac{\pi }{2}{\beta }^{3}s\right)e^{-\alpha s-\pi {\beta }^2{s}^2/4},$$

$${\Sigma }^2(L)=⟨(n(L,e)-L{)}^2⟩,$$

$${\Sigma }_{\text{SP}}^2(L)=\frac{L}{2}+\frac{1}{8}(1-e^{-4L}),$$

$${\Sigma }_{\text{GOE}}^2(L)=\frac{2}{{\pi }^2}\left[\ln (2\pi L)+\gamma +1-\frac{{\pi }^2}{8}\right]+\mathcal{O}(L^{-1}),$$

$$P(\psi )={\left(\frac{A}{2\pi }\right)}^{1/2}{e}^{-A{\psi }^2/2}.$$

$$W_i=\sum _{k=1}^{N_{\text{dot}}}\underset{B_k}{\iint }dydx|{\psi }_i(x,y){|}^2,$$

$$\frac{A_{\circ }}{A}=\frac{\pi }{4L_x{L}_y}\sum _{k=1}^{N_{\text{dot}}}{b}_k^2.$$

$${\chi }_{\nu }^2(z)=\frac{(\nu /2{)}^{\nu /2}}{\Gamma \left(\nu /2\right)}{z}^{\nu /2-1}{e}^{-\nu z/2},$$

$${\nu }_{\text{max}}=e_{\text{max}}\frac{4}{\pi }\frac{A_{\circ }}{A}.$$

$$⟨x,y|{\phi }_{\text{dot}}⟩=C\exp \left(-\frac{M\omega }{2\hslash }[(x-a_1{)}^2+(y-e_1{)}^2]\right),$$

$$\{u_i{\}}_{i=1,\dots ,N}=\{u_1={\phi }_{\text{dot}},u_2={\phi }_2,\dots ,u_N={\phi }_N\}.$$

$$|{\psi }_1^{(\text{OPW})}⟩$$

$$|{\psi }_i^{(\text{OPW})}⟩$$

$$⟨{\psi }_1^{(\mathrm{OPW})}|{\psi }_i^{(\mathrm{OPW})}⟩$$

Full text

I Introduction

Self-assembled semiconductor quantum dots (QDs) have attracted a lot of attention due to their potential for fundamental research, such as superradiance Jahnke et al. (2016) and cavity-quantum electrodynamics Lodahl et al. (2004), as well as device applications, such as single-photon sources Michler et al. (2000) and low-threshold microlasers Wiersig et al. (2009). QDs allow for a carrier confinement in all three dimensions on the nanometer scale with a discrete atomic-like density of states. The corresponding single-particle energy eigenstates are spatially strongly localized at the respective QD position. The morphology of such QD states has been studied comprehensively in the literature with pseudopotential theory Bester, Nair, and Zunger (2003); Bester and Zunger (2003), ${\bf k\cdot p}$ models Rodt et al. (2005), tight-binding models Sheng, Cheng, and Hawrylak (2005); Baer et al. (2005), and density functional theory Nair and Masumoto (2000).

In self-assembled QDs grown in the Stranski-Krastanow growth mode, the QD states are located energetically below a quasicontinuum of delocalized states, which corresponds to the two-dimensional (2D) motion in a thin wetting layer (WL). Because of carrier-carrier and carrier-phonon scattering the WL is of great importance for QD devices if the excitation involves the quasicontinuum Nielsen, Gartner, and Jahnke (2004). Computing the above many-particle processes requires the knowledge not only of the QD states but also of the spatially-extended WL states. It is common that the latter are described by plane waves for their in-plane part Magnusdottir et al. (2002); Chang, Chuang, and Holonyak Jr. (2004), thereby omitting the QDs’ influence on the WL states. As a consequence, the QD and WL states are no longer orthogonal to each other, which potentially adulterates the QD-WL Coulomb coupling Schneider, Chow, and Koch (2001). The orthogonalized plane wave (OPW) method is an attempt to remedy this shortcoming by a supplementary orthogonalization Bockelmann and Egeler (1992); Ferreira and Bastard (1999); Schneider, Chow, and Koch (2001); Nielsen, Gartner, and Jahnke (2004). However, the resulting WL states are still approximated by plane waves, but with local modifications at the position of the QDs. Calculations without the above approximations exist for single QDs Markussen et al. (2006); Popescu, Bester, and Zunger (2009); Michael et al. (2018), but the considered spatial domain around the QD is small, so the WL states are treated only locally and their properties are not in the focus.

In the field of quantum chaos Gutzwiller (1990); Stöckmann (2000); Wimberger (2014), plane wave-like energy eigenstates with a regular morphology appear in integrable systems. Examples of such analytically solvable systems are the rectangular and the circular billiard. The term billiard refers to a potential-free region enclosed by reflecting walls, in which a point mass moves along a straight line until it hits the boundary and bounces back with the angle of reflection equal to the angle of incidence. Most billiards are non-integrable and exhibit chaos, i.e. there exist trajectories with a sensitive dependence on initial conditions. Some geometries exhibit full chaos, i.e. almost all trajectories exhibit a sensitive dependence on initial conditions. The quantum properties of chaotic billiards have been studied experimentally in, e.g. optical microcavities (for a review see Ref. Cao and Wiersig (2015)) and ballistic semiconductor microstructures Marcus et al. (1992). The latter are sometimes also called QDs but are very different from the above self-assembled QDs, for instance because of the size (compared to the wavelength) and the lack of a WL.

The aim of this article is to demonstrate that even in the simple model of a rectangular-shaped WL with one or a few lens-shaped self-assembled QDs in the effective-mass approximation, the WL states are typically rather complex and do not resemble (orthogonalized) plane waves, as illustrated in Fig. 1. We use methods and concepts from the field of quantum chaos to study the complexity of the states.

The paper is organized as follows. Section II introduces the model. Subsequently, Sec. III describes the numerical method to calculate energy eigenstates and eigenvalues. The morphology of these states is qualitatively discussed in Sec. IV. Section V characterizes the system via spectral statistics. Sections VI and VII provide a quantitative study of the global and local spatial properties of the WL states. The results are compared to the OPW method in Sec. VIII. Finally, the conclusion is given in Sec. IX.

II The quantum-dot wetting-layer model

We describe shallow lens-shaped QDs sitting on a thin, quasi-2D WL in a one-band effective-mass approximation for the conduction-band electrons, ignoring the effects of band mixing, mass anisotropy, and spin-orbit coupling. A constant effective mass $M$ is used assuming that the QDs and the WL are of the same material and embedded in a different bulk material. For the thin WL we consider only the first subband exhibiting no nodes in $z$-direction. We determine the single-particle energy eigenstates in the envelope function approximation, and consider only the in-plane component $\psi(x,y)$ of the wave functions. The QD confinement is treated by a potential $V(x,y)$ which naturally also influences the spatially-extended states in the WL. The WL is considered as an area of finite size, also called Mesa (see, e.g. Gschrey et al. (2013)), which is chosen to be a region of rectangular shape $\Omega=[0,L_{x}]\times[0,L_{y}]$. Following Ref. Schneider, Chow, and Koch (2001) we describe the confinement in the WL by an infinitely high potential well in 2D.

It has been clearly demonstrated in theoretical Wojs et al. (1996); Riel (2008) and experimental Park et al. (1998); Teichmann,Wenderoth,Prüser, Pierz,Schumacher, and Ulbrich (2013) investigations that the single-particle energy spectrum of a shallow lens-shaped QD can be accurately reproduced by a truncated parabolic confinement potential with rotational symmetry in the $x,y$-plane. Thus the QD is a finite size version of an isotropic 2D harmonic oscillator. The energy levels of the isolated QD are therefore given by $\hbar\omega(n+1)$ with degree of degeneracy $n+1$, where $n=0,1,2,\ldots,n_{\text{max}}$. For a single dot the potential is depicted in Fig. 1. Since the depositioning process in general leads to a spatial distribution of QDs on the WL we have to incorporate the $k$-th QD confinement potential to the overall WL-potential by summing over all $N_{\text{dot}}$ non-overlapping QDs

[TABLE]

where the set of parameters $\{a_{k},e_{k},d_{k},b_{k}\}$ characterizes the $k$-th QD confinement potential with its minimum at the coordinates $(a_{k},e_{k})$, width $b_{k}$, and depth $d_{k}b_{k}^{2}/4$.

Table 1 summarizes the sets of parameters studied in the next sections. To avoid nongeneric behavior, we choose the aspect ratio $L_{y}/L_{x}$ to be an irrational number by using the reciprocal golden ratio $\Phi=(1+\sqrt{5})/2$. In our numerics we express the parameters $a_{k}$, $e_{k}$, $b_{k}$, $L_{x}$, and $L_{y}$ in nm and, using $\hbar=1=M$, express energies in nm*-2*. First, we study the case of a single QD ($N_{\text{dot}}=1$) with three QD states. Three confined states are realistic for several applications such as QD lasing, see, e.g. Gies et al. (2007). Within the model, for fixed $b_{1}$ the maximum number of QD states can be adjusted by the parameter $d_{1}$. Even though the situation with only one QD does not serve as a generic type of a semiconductor with QD densities around $10^{10}\text{cm}^{-2}$ it is experimentally feasible in a Mesa, as in Ref. Gschrey et al. (2013), and offers the opportunity to study the impact on the system properties on the most simplified level. The first sample system is referred to as $\mathcal{A}(1)$. The QD position $(a_{1},e_{1})$ is chosen such that the QD is not exactly at the center to avoid nongeneric effects due to symmetry. Leaving $d_{1}$ and $b_{1}$ unchanged we extend the area $\Omega$ by multiplying $L_{x}$ and $L_{y}$ of system $\mathcal{A}(1)$ by $t=2$, $3$, and $4$. Also the QD position is changed by the same factor such that $L_{x}/a_{1}$ and $L_{y}/e_{1}$ remain the same for all $t$. In the same order, these systems are labeled by $\mathcal{A}(2)$, $\mathcal{A}(3)$, and $\mathcal{A}(4)$. Moreover, another system containing an ensemble of four identical, non-overlapping QDs is denoted by $\mathcal{B}$. System $\mathcal{C}$ possesses the same parameter settings as system $\mathcal{A}(1)$ but with an adjusted parameter $d_{1}$ reducing the number of QD states to one.

The restriction to the first subband implies a cut-off energy $E_{\text{cut-off}}$ below which the system behaves quasi-2D. This is analog to the case of microwave billiards, which are popular experimental systems in the field of quantum chaos, see, e.g. Stöckmann (2000). The cut-off energy is determined by the condition that the second subband appears which happens in our case at $E_{\text{cut-off}}=3\pi^{2}\hbar^{2}/(2Md_{\text{WL}}^{2})$ where $d_{\text{WL}}$ is the thickness of the WL. Assuming a thin WL with $d_{\text{WL}}=0.4\,$nm we get $E_{\text{cut-off}}\approx 92\,$nm*-2*. In the following we will consider only energy levels below this cut-off energy.

We mention that the model introduced above is related to the model in Ref. Markussen et al. (2006) with, in the context of the present paper, two important differences: (i) in Ref. Markussen et al. (2006) only a single QD is studied and (ii) this dot is sitting on the center of a circular-shaped WL. This composed system is nongeneric as it possesses a rotational symmetry. The resulting integrability is convenient for computations but it leaves no room for a complex morphology of energy eigenstates. In strong contrast, our model is nonintegrable and allows to study generic features of QD-WL systems.

At first glance at the lower part of Fig. 1 there is also a superficial similarity to the Sinai billiard Sinai (1970), a square billiard with a circular wall located at its center, and the Šeba billiard Šeba (1990), a rectangular billiard containing a point scatterer (see also the application to ballistic QDs in Ref. Röder and Hentschel (2010)). However, our QD-WL model is very different as the potential is attracting, nonuniform, continuous, and finite.

Finally, our model is related to the quasi-2D microwave billiard with randomly placed spherical caps which has been studied in Ref. Barkhofen, Metzger, Fleischmann, Kuhl, and Stöckmann (2013) in the context of branch flows in random potentials. These caps remind strongly on our lens-shaped QDs positioned on the WL but the caps are placed inside the resonator. The latter results in repulsive local potentials whereas in our case the local potentials are attractive. The energy eigenstates and eigenvalues have not been discussed in Ref. Barkhofen, Metzger, Fleischmann, Kuhl, and Stöckmann (2013).

III Computation of eigenstates and eigenvalues

We consider the quantum-mechanical eigenvalue problem

[TABLE]

with the single-particle Hamiltonian

[TABLE]

and the eigenstates $\psi_{i}$ fulfilling zero Dirichlet boundary conditions on the WL boundary $\partial\Omega$. We solve it numerically by expressing ${\hat{H}}$ in a suitable orthonormal basis $\{\varphi_{n}\}$ and diagonalizing the associated Hamiltonian matrix $H=(H_{n,m})$. This procedure supplies the energy eigenvalues $E_{i}$ which, due to considering a finite region $\Omega$, are discretized, and the coefficients $c_{n}^{(i)}=\langle\varphi_{n}|\psi_{i}\rangle$ required for the series expansion of the eigenstates

[TABLE]

Here, it is convenient to choose the eigenbasis $\{\varphi_{n}\}$ of $\hat{H}_{0}$ resulting from removing the QD potential from the Hamiltonian in Eq. (4). According to the boundary conditions the basis is given by the plane wave-like (checkerboard-like) states

[TABLE]

with wave numbers $k_{x}=n_{x}\pi/L_{x}$, $k_{y}=n_{y}\pi/L_{y}$, and $n_{x},n_{y}\in\mathbb{N}^{+}$. The multi-index $n=\{n_{x},n_{y}\}$ can be assembled by using an appropriate pairing function facilitating a bijectively mapping of two integers onto one and thus providing a certain ordering scheme. Because $n_{x},n_{y}\in\mathbb{N}^{+}$ we are using a bijective map $h:\mathbb{N}^{+}\times\mathbb{N}^{+}\rightarrow\mathbb{N}^{+}$ introduced by Hopcroft and Ullman Hopcroft and Ullman (1979).

By construction, the matrix elements of the first term in Eq. (4), recognized as $\hat{H}_{0}$, can be evaluated analytically

[TABLE]

with $E_{n_{x},n_{y}}=\hbar^{2}\left(k_{x}^{2}+k_{y}^{2}\right)/2M$. The matrix elements of the second term in Eq. (4) can be written as

[TABLE]

with $B_{k}$ constituting the region where $V_{k}\neq 0$; cf. Fig. 1. The numerical evaluation of the double integral (8) is done in local polar coordinates by using scipy.integrate.nquad from the SciPy package Jones et al. (01). The upper bound on the number of subintervals used in the adaptive algorithm has been set to 100 to ensure convergence. For the computation of eigenvectors and eigenvalues of the real symmetric matrices we use numpy.linalg.eig and numpy.linalg.eigvalsh from the same package with default settings. The resulting eigenstates and eigenvalues are ordered such that $E_{1}\leq E_{2}\leq\ldots$.

IV Morphology of energy eigenstates

This section provides a qualitative discussion of the numerically obtained eigenstates. For a quantitative analysis we refer to Secs. VI and VII. As representatives we choose systems $\mathcal{A}(1)$ and $\mathcal{B}$, cf. Table 1. For the single-QD system $\mathcal{A}(1)$ the numerically computed QD states are shown in Fig. 2. We define the QD states as those states having negative energy; see Eqs. (1)-(2). The QD states are localized in the QD and obey the expected behavior of an isotropic 2D harmonic oscillator. The ground state ($s$-shell) is concentrated in the center of the QD and its energy is well approximated by a 2D harmonic oscillator, shifted by $V_{\text{min}}=-d_{1}b_{1}^{2}/4$. The next two consecutive states ($p$-shells) possess the typical appearance with probability densities leaking into the classically prohibited region. The eigenenergies are slightly splitted due to the presence of the WL boundary.

WL states have nonnegative energies and are therefore spatially extended. Figure 3 shows that these states are strongly affected by the presence of the QD. While some states like $\psi_{8},\psi_{50}$, and $\psi_{1356}$ display an irregular or nonuniformly spatial structure, some others like $\psi_{475},\psi_{866}$, and $\psi_{1197}$ display regular and almost uniform checkerboard-like densities. Occasionally, for system $\mathcal{A}(1)$ we also find characteristic states like $\psi_{429}$ which might be classified as “bouncing-ball states” Bäcker, Schubert, and Stifter (1997).

Furthermore, by taking a closer look at the WL states in Fig. 3 we notice that some of them exhibit a higher probability in the region of the QD, e.g. $\psi_{7},\psi_{8}$, and $\psi_{13}$. Some others like $\psi_{11}$ and $\psi_{37}$ show a diminished probability. The latter states seem to tend to “elude” the QD. The examination of this effect on a statistical level is the subject of Sec. VII.

Whereas the WL states of system $\mathcal{A}(1)$ constitute a set of several types, for system $\mathcal{B}$ Fig. 4 shows that increasing the number of QDs implies a loss of diversity with respect to the spatial morphology. Checkerboard-like states have disappeared and generally the spatial structure of the higher-excited states looks more irregular. Additionally, due to the QD configuration in system $\mathcal{B}$ we come across the opportunity to observe coupled QD states when the inter-dot distance is small enough, like in Figs. 4(a) and (b).

To conclude the qualitative discussion, we can say that the morphology of the WL states in this simple QD-WL model is much more complex than that of plane waves.

V Spectral analysis

Here, we investigate the information contained in the energy spectrum utilizing statistical measures developed in the field of quantum chaos. Two of them are the so-called level spacing or nearest-neighbor spacing distribution $P(s)$ which reveals the short-range correlations of the energy levels and secondly the number variance $\Sigma^{2}(L)$ revealing the long-range level correlations, see, e.g. Stöckmann (2000); Wimberger (2014). Both quantities are dealing with the fluctuations within a discrete level sequence $\{E_{i}\}$ with $E_{1}\leq E_{2}\leq\ldots$. We start with the level density $\rho(E)=\sum_{i}\delta(E-E_{i})$ and the integrated level density

[TABLE]

giving the number of levels up to the energy $E$. Here $\Theta$ denotes the Heaviside step function. $N(E)$ can be segregated into two parts, a smooth part ${\bar{N}}(E)$ and a fluctuating one

[TABLE]

Determining the smooth part ${\bar{N}}(E)$ is in general a nontrivial task Gómez et al. (2002). However, for 2D quantum billiards of arbitrary shape with area $A$ and circumference $L$, the smooth part of the integrated level density, for not too small energies $E$, can be very well approximated by the generalized Weyl’s law Gutzwiller (1990)

[TABLE]

We assume (and verify below) that this law can be also applied to our QD-WL model with $A=L_{x}L_{y}$ and ${L}=2(L_{x}+L_{y})$ if the levels of the QD states are omitted. This is justified as Weyl’s law is only valid for not too small energies and, moreover, a few events are insignificant to the overall statistics. The validity of Weyl’s law is demonstrated for system $\mathcal{A}(1)$ in Fig. 5. The spectrum is calculated for a Hamiltonian matrix of size $3000\times 3000$. Up to the $2100$th energy level Weyl’s law constitutes a marvelous approximation to the integrated level density. Deviations beyond this level can be attributed to truncating the Hilbert space.

With the aid of ${\bar{N}}(E)$, the spectrum $\{E_{i}\}$ can be mapped onto the so-called unfolded spectrum $\{e_{i}\}$ through $e_{i}={\bar{N}}(E_{i})$. This procedure separates the local level fluctuations from an overall energy dependence. Owing to the unfolding the sequence of nearest-neighbor energy spacings $\{s_{i}=e_{i+1}-e_{i}\}$ constitutes a set for which the mean level spacing $\langle s\rangle\approx 1$. The corresponding probability distribution is called $P(s)$. Uncorrelated energy levels follow the Poisson statistics with nearest-neighbor spacing distribution

[TABLE]

In the field of quantum chaos the Poisson statistics describe the level statistics of generic integrable systems. This is also true for the special case of a rectangular billiard with a generic aspect ratio Marklof (1998) which is nothing else than our QD-WL model without QDs having plane wave-like energy eigenstates (6). In the presence of a QD deviations from Poisson statistics can be observed, e.g. in Fig. 6(a) for system $\mathcal{A}(1)$. These deviations therefore reveal correlations between adjacent energy levels, i.e. short-range level correlations. In particular, $P(s)$ is significantly reduced for small values of $s$. Hence, it is less probable to observe level crossings when a parameter is varied.

Figure 6(a) also shows that $P(s)$ cannot be described by the Gaussian orthogonal ensemble (GOE) of random-matrix theory which applies to classically fully chaotic systems with time-reversal symmetry. The nearest-neighbor spacing distribution of the GOE is well approximated by the Wigner surmise (see, e.g. Wimberger (2014))

[TABLE]

A better agreement for system $\mathcal{A}(1)$ is achieved with the so-called Semi-Poisson statistics Wiersig (2002); Bogomolny, Gerland, and Schmit (2001) with nearest-neighbor spacing distribution

[TABLE]

Like the Wigner surmise $P_{\text{SP}}(s)$ shows a linear behavior at small $s$. But for large $s$ it exhibits an exponential fall-off like in the Poisson case. The Semi-Poisson statistics is a model for intermediate spectral statistics Bogomolny, Gerland, and Schmit (2001) and appears, e.g. in pseudointegrable systems Wiersig (2002); Gorin and Wiersig (2003); Wiersig and Carlo (2003) and multi-walled carbon nanotubes Ahn et al. (2003). Furthermore, it has also been controversially discussed in the context of the Šeba billiard Tudorovskiy, Kuhl, and Stöckmann (2010).

Next, we check the Brody distribution Brody et al. (1981)

[TABLE]

with

[TABLE]

and the gamma function $\Gamma$. The Brody distribution is a heuristic interpolation between Wigner surmise ($\beta=1$) and Poisson statistics ($\beta=0$), see Ref. Stöckmann (2000). Figure 6(a) shows a very good agreement with the numerical data for $\beta=0.35$, which is even slightly better than for the Semi-Poisson statistics.

We have also tested the Berry-Robnik distribution Berry and Robnik (1984)

[TABLE]

with the complementary error function erfc and $\alpha=1-\beta$. This distribution also interpolates between Wigner surmise ($\beta=1$) and Poisson statistics ($\beta=0$). In contrast to the Brody distribution, it has a theoretical foundation for systems with a mixed classical phase space where regular/integrable and chaotic dynamics coexist. The spectra associated with the corresponding regions in phase space are assumed to be statistically independent, and their mean level spacing is determined by the size of the regions. Even though our QD-WL model should be in the generic class of systems with mixed phase space, the Berry-Robnik distribution does not give a satisfactory fit to the data for $s\leq 1$ (not shown). This failure of the Berry-Robnik distribution is consistent with other numerical studies Prosen and Robnik (2001); Zhou et al. (2010). It results from neglecting dynamical tunneling between regular and chaotic regions Bäcker et al. (2011).

We observe that the agreement of the numerical data and the Brody distribution remains reasonable when the WL area is increased by going through the systems $\mathcal{A}(t)$ with $t=2,3,4$ (keeping the size of the Hamiltonian matrix constant), see Table 1. For $t=3$ (not shown for $t=2$ and $4$) this can be clearly observed in Fig. 6(b). It can also be seen that in this case the Semi-Poisson statistics gives a slightly better agreement.

Following the discussion in Refs. Ahn et al. (2003); Tudorovskiy, Kuhl, and Stöckmann (2010) we also consider the spectral statistics on various energy windows. Dividing the overall spectrum of $2000$ energies into smaller subsets reveals that $P(s)$ of system $\mathcal{A}(1)$ gets different contributions from separate energy windows as shown in Figs. 7(a)-(c). In the lowest energy window $P(s)$ is equally well described by the Semi-Poisson statistics and by a Brody distribution with fitting parameter $\beta=0.5$. For increasing energies the distribution is slightly better described by Brody distributions with $\beta$ decreasing to $0.35$ and $0.3$, i.e. the distribution slowly shifts closer to the Poisson statistics. This is reasonable as the influence of the QD potential is reduced for high energies. If the QD potential can be ignored then the system is effectively a rectangular billiard obeying Poisson statistics.

The tendency to Poisson statistics for higher energies is also observable when more than one QD is present. Figures 7(d)-(f) show the case of four QDs in system $\mathcal{B}$ with altogether twelve QD states using a $4000\times 4000$ Hamiltonian matrix. For low energies, the nearest-neighbor spacing distributions [Fig. 7(d)] can be characterized by the Brody distribution with $\beta=1$ which equals the Wigner surmise. By increasing the energy a slow transition to $\beta=0.9$ [Fig. 7(e)] and $\beta=0.6$ [Fig. 7(f)] is observed. In the high-energy regime also the Semi-Poisson statistics describes the numerical data very well. Note that system $\mathcal{B}$ is significantly closer to the statistics of random-matrix theory than system $\mathcal{A}(1)$ even though the QD density differs just by a factor of $1.28$. This is consistent with our findings on the energy eigenstates in Sec. IV, see Figs. 3 and 4.

To study long-range level correlations we consider the number variance

[TABLE]

measuring the local variance of the number $n(L,e)=N(e+L/2)-N(e-L/2)$ of energy levels in the interval $\left[e-L/2,e+L/2\right]$. The bracket $\langle\ldots\rangle$ denotes the average over all spectral positions $e$. For the Poisson statistics one gets $\Sigma^{2}_{\text{P}}(L)=L$, for the Semi-Poisson statistics Wiersig (2002); Bogomolny, Gerland, and Schmit (1999)

[TABLE]

and for the GOE (see, e.g. Wimberger (2014))

[TABLE]

with the Euler constant $\gamma\approx 0.57722$.

In the case of the single-QD system $\mathcal{A}(1)$ Fig. 8 shows a remarkable agreement with the number variance of the Semi-Poisson statistics in Eq. (19). For $L>6$ we observe deviations which might be the generic saturation behavior due to a finite number of energy levels used, see, e.g. Ref. Bäcker, Steiner, and Stifter (1995). For system $\mathcal{A}(3)$ with one QD on a three times larger WL, Fig. 8 shows a shift towards the GOE in Eq. (20). For the 4-QD system $\mathcal{B}$ we find that $\Sigma^{2}(L)$ is even closer to $\Sigma^{2}_{\text{GOE}}(L)$. To summarize, the short- and long-range level correlations are intermediate between Poisson statistics and GOE with a tendency towards GOE when the size of the WL or the number of QDs is increased. This shows that the WL states in our model show features known from irregular energy eigenstates in nonintegrable systems.

VI Amplitude distribution

After the qualitative discussion on the morphology of the WL states in Sec. IV, we now provide quantitative measures for the observed complexity. According to Berry Berry (1977), the (real-valued) amplitude $\psi$ of an individual “chaotic” or irregular wave function behaves like a Gaussian random variable. Hence, the probability of finding a value of it at any arbitrary given point in a region of area $A$ is McDonald and Kaufman (1988)

[TABLE]

To determine $P(\psi)$ for a given WL state we evaluate its wave function $\psi(x,y)$ at $128\times 128$ points in the $x,y$-plane. The outermost points of the WL region $\Omega$ are excluded to avoid an artificial peak at $\psi=0$ due to the zero Dirichlet boundary conditions. Then $P(\psi)$ is generated as a histogram with $90$ bins in $[-0.5,\ldots,0.5]$, normalized to unity.

For the single-QD system $\mathcal{A}(1)$ Fig. 9(a) shows two examples, $\psi_{503}$ and $\psi_{1356}$, which are fairly well described by Eq. (21), with small deviations in the center and in the left and right wings. However, there are also many examples in this system that do not follow Eq. (21) (not shown), for instance the bouncing-ball state $\psi_{429}$ and the checkerboard-like states $\psi_{866},\psi_{1197}$; cf. Fig. 3. It is well known that such regular states are not described by a Gaussian random variable McDonald and Kaufman (1988); Shigehara (1994).

For the 4-QD system $\mathcal{B}$ we find that nearly all among twenty randomly chosen sample states between $\psi_{492},\ldots,\psi_{1581}$ show a remarkably good agreement with Eq. (21). Two examples are shown in Fig. 9(b).

Hence, as in the previous section, we find that the QD-WL model shows features known from nonintegrable and chaotic systems. Again, the system with more than one QD appears to be more “chaotic”, in the sense that its WL states behave more like Gaussian random variables.

VII Localization

In this section we take up a question raised in Sec. IV about the WL states’ tendency to avoid or seek the region of the QD(s). For this purpose, we define the localization of the $i$-th eigenstate as

[TABLE]

where, again, $B_{k}$ is the region where the confinement potential $V_{k}$ of the $k$-th QD is nonzero. These double integrals are summed over all QDs and, like Eq. (8), are evaluated by using scipy.integrate.nquad with upper bound on the number of subintervals being set to 100. For comparison we consider the ratio of the area of QDs ($A_{\circ}$) and the area of the WL ($A$) which can be written as

[TABLE]

For system $\mathcal{A}(1)$ we compute the localization $W_{i}$ for the WL states with $i=4,\ldots,1000$. We consider here a smaller amount of states than shown in Fig. 3 due to the high sensitivity of the $W_{i}$ on the convergence of the series expansion [Eq. (5)], which here again contains $3000$ terms. For comparison, we calculate $W_{i}$ also for system $\mathcal{A}(1)$ under withdrawal of the QD having the plane wave-like energy eigenstates in Eq. (6). The resulting probability distributions $P(W)$ are depicted in Fig. 10. In the system without QD the distribution essentially shows one sharp peak at $\langle W\rangle\approx 0.03\approx A_{\circ}/A$ [see Eq. (23)] quickly declining to zero and two flattened side lobes to its left and right. Inserting a QD considerably broadens the distribution $P(W)$ without changing the mean $\langle W\rangle$ which is therefore again very close to the area ratio in Eq. (23). That means a QD does not change the average localization but it does increase the fluctuations resulting in more WL states that tend to avoid the QD region but also more WL states that tend to seek it.

To explain the behavior of $P(W)$ in the QD-WL model we first note that if we would introduce a uniform absorption within the QD region, the resulting decay rate $\gamma_{i}$ of the $i$-th eigenstate would be in first-order perturbation theory for weak coupling proportional to the localization $W_{i}$. Hence, $\gamma_{i}/\langle\gamma\rangle=W_{i}/\langle W\rangle$. It has been shown that for weakly open systems (see, e.g. Grabsch, Savin, and Texier (2018)) the lifetime statistics is given by a $\chi^{2}$ distribution

[TABLE]

with normalized decay rate $z=\gamma_{i}/\langle\gamma\rangle$ and number of decay channels $\nu$. Equation (24) is the many-channel generalization of the Porter-Thomas distribution Porter and Thomas (1956). The $\chi^{2}$ distribution applies to integrable and chaotic systems provided that the coupling to the environment is weak and that the decay channels can be treated as independent Gaussian variables, see Gorin et al. (1997) and references therein. We therefore assume $P(W)\approx\chi^{2}_{\nu}(W/\langle W\rangle)$ and use the number of decay channels $\nu$ as a fitting parameter. An upper bound for $\nu$ can be crudely estimated in the following way. First, we assume that the size of a decay channel in 2D real space is on average $(\lambda/2)^{2}$ with the wavelength $\lambda$ corresponding to the considered energy. Second, we count how many of such areas fit into the area $A_{\circ}$ of the QD region. With the largest (unfolded) energy $e_{\text{max}}$ used in the statistics and Eq. (23) we find

[TABLE]

It can be seen in Fig. 10 that the numerically obtained $P(W)$ of $\mathcal{A}(1)$ is rather well fitted by the $\chi^{2}$ distribution for $\nu=12$. This is consistent with the upper bound $\nu_{\text{max}}\approx 39$ using $e_{\text{max}}=1000$.

The system $\mathcal{A}(1)$ without QD, see Fig. 10, cannot be well fitted by a $\chi^{2}$ distribution. We explain this by the fact that the decay channels here are not independent as the amplitudes of a regular state within a small QD region are strongly correlated. In particular, it means that a regular state cannot avoid to overlap with the QD region which implies that $P(W)$ is significantly reduced for small $W$.

Figure 11 shows $P(W)$ for the 4-QD system $\mathcal{B}$. Here, the agreement with the $\chi^{2}$ distribution is even slightly better than for system $\mathcal{A}(1)$. The fitted value of $\nu=21$ is consistent with the upper bound $\nu_{\text{max}}\approx 30$, using $e_{\text{max}}=1000$ in Eq. (25).

Based on the so far observed statistical properties of the WL states it is natural to ask whether the localization in the QD(s) is different from the localization properties outside the QD(s). To answer this question we choose now $B_{1}$ in Eq. (22) for the single-QD system $\mathcal{A}(1)$ to be a spatially shifted region of the same shape and size, still in $\Omega$, but non-overlapping with the QD. Figure 12 shows the resulting localization probability distribution with the shifted region centered at $(a_{1},e_{1})=(3.18,13.24)$. One can see that the general shape of $P(W)$ is as before with, however, an increased $\nu=18$. Similar observations we have made for regions centered at $(a_{1},e_{1})=(7.52,11.74)$, $(2.64,2.57)$, and $(6.98,2.36)$. We therefore conclude that the statistical properties of the WL states are in first approximation independent of the position in the WL which is a feature shared with eigenstates in chaotic systems.

VIII Orthogonalized plane waves

In this section we contrast the WL states with the OPWs. The OPW method tries to incorporate the impact of a QD on the WL states by orthogonalizing plane waves to the QD state(s). Following Ref. Schneider, Chow, and Koch (2001), we consider the single-QD system $\mathcal{C}$ in Table 1 which possess only one QD state due to a reduced parameter $d_{1}$. This state is well approximated by the ground state of the isotropic 2D harmonic oscillator given by

[TABLE]

with $C=\left(M\omega/\pi\hbar\right)^{1/2}$. In Eq. (26) the parameter $d_{1}$ is incorporated via the frequency $\omega=\sqrt{2d_{1}/M}$. The set of functions that need to be orthogonalized are those of the infinitely high potential well sorted in ascending order by their energy, with its ground state substituted by $\varphi_{\text{dot}}$:

[TABLE]

It should be noted that in a set $\{u_{i}\}$ containing $u_{1}=\varphi_{1}$ the QD state $\varphi_{\text{dot}}$ could also constitute an additional state.

The approximated states are determined as

[TABLE]

where $i=2,\ldots,N$ and $N_{j}^{(\text{OPW})}$ with $j=1,\ldots,N$ is a normalization constant of the $j$-th state, such that the state is normalized to unity in $\Omega$. But unlike the exact eigenstates the OPWs $\psi_{i}^{(\mathrm{OPW})}$ are just orthogonal to the ground state $\psi_{1}^{(\mathrm{OPW})}$

[TABLE]

but not necessarily to each other.

Two examples of OPWs are shown in Figs. 13(a) and (b). As expected from the above orthogonalization procedure, these states are plane waves but with a modification in the QD region. It is obvious that such OPWs cannot approximate the complex morphology of WL states in Figs. 3 and 4. Figures 13(c) and (d) show two exact WL states with a roughly regular, checkerboard-like structure. In these specific examples, it seems that the OPWs in Figs. 13(a) and (b) almost can approximate the WL states. However, the wave function inside the respective QD region, which is important for the carrier capture rate calculations Schneider, Chow, and Koch (2001), is not correctly predicted.

There is another troublesome effect becoming noticeable for increasing energy. While in the exact treatment of the QD-WL model even for single-QD systems the QD does not lose its impact on the WL states even for rather high energies, in the OPW approach the impact quickly lessens until it finally ceases. This results from a decreasing overlap of $\langle\psi_{1}^{(\text{OPW})}|u_{i}\rangle$ in Eq. (29) due to stronger spatial oscillations of the $u_{i}$ for higher energies.

As the OPW method fails here to reproduce the complex pattern of WL states, we expect that this approach does not describe the carrier capture relaxation process into QDs better than the simpler approach of using just plane waves. This statement is consistent with the findings in Ref. Markussen et al. (2006).

IX Conclusion

In this paper we have studied the properties of wetting-layer states in a simple quantum-dot-wetting-layer model. The spectral analysis via the nearest-neighbor spacing distribution and the number variance showed that both the short- and long-range level correlations are intermediate between Poisson statistics and the Gaussian orthogonal ensemble of random-matrix theory, with a tendency towards the latter when the size of the wetting layer or the number of quantum dots is increased. By dividing the spectra into distinct energy windows we demonstrated that the nearest-neighbor spacing distribution shows a slow transition towards Poisson statistics when the energy is increased.

Moreover, we demonstrated numerically that the spatial structure of the wetting-layer states shows a complex morphology. Depending on the number of dots we observed distinct classes of wetting-layer states. In the simplest case of a rectangular-shaped wetting layer with one quantum dot the possible behavior ranges from checkerboard-like states up to states with unequally or more irregular pattern. By increasing the number of quantum dots we observed that checkerboard-like states disappear. We showed that the irregular looking states exhibit probability amplitude distributions very similar to that of a Gaussian random variable.

Furthermore, we investigated the localization behavior of wetting-layer states in the quantum-dot region(s). Due to the presence of a dot, the probability distribution of the localization is broadened without changing the mean value. Hence, the average localization does not change but there are more states that tend to avoid the quantum-dot region but also more states that tend to seek it. Studying the localization in some distance from the quantum dot revealed a similar behavior. Interestingly, the localization distribution in the presence of quantum dot(s) is well fitted by a $\chi^{2}$ distribution.

The observed localization behavior can be important for capture time calculations, such as in Ref. Schneider, Chow, and Koch (2001). In this context, we demonstrated that the orthogonalized plane waves are not capable of accurately approximating the wetting-layer states.

In the studied parameter regime we have not seen any signatures of Anderson localization Anderson (1958) which could be observable in our model for a very large wetting-layer area with many randomly-placed quantum dots.

In this work we considered quantum dots on a wetting layer of finite rectangular shape. Even this simple shape, which in the absence of any quantum dot would have plane wave-like energy eigenstates, give rise to wetting-layer states with a complex morphology. We showed that increasing the size of the wetting layer even increases the complexity of the wetting-layer states in terms of their level statistics. We expect that more complicated shapes of the wetting layer that might appear in experiments do not change the qualitative behavior as the energy eigenstates in the absence of quantum dots already have an irregular spatial structure, see, e.g. Stöckmann (2000). Following the same line of reasoning, we also expect that deviations of the quantum-dot confinement potential from the truncated parabolic shape do not change the qualitative behavior of the wetting-layer states.

Based on the insights presented in this work, we surmise that an improved approach to approximate wetting-layer states may consist in using random superpositions of plane waves. The statistical properties of such states, see, e.g. O’Connor, Gehlen, and Heller (1987), are closer to those of the wetting-layer states. We therefore expect that random superpositions of plane waves are more tailored to the needs of the actual wetting-layer states than (orthogonalized) plane waves.

Acknowledgments

We would like to thank J. Kullig, B. Gulyak, and G. Kasner for helpful discussions.