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

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 on 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.


Introduction
Self-assembled semiconductor quantum dots (QDs) have attracted a lot of attention due to their potential for fundamental research, such as superradiance [1] and cavity-quantum electrodynamics [2], as well as device applications, such as single-photon sources [3] and low-threshold microlasers [4]. 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 [5,6], k · p models [7], tight-binding models [8,9], and density functional theory [10].
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 [11]. 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 [12,13], 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 [14]. The orthogonalized plane wave (OPW) method is an attempt to remedy this shortcoming by a supplementary orthogonalization [11,[14][15][16]. 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 [17][18][19], 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 [20][21][22], 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 nonintegrable 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 [23]) and ballistic semiconductor microstructures [24]. The latter are sometimes also called QDs but are very different from the above selfassembled 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 figure 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 2 introduces the model. Subsequently, section 3 describes the numerical method to calculate energy eigenstates and eigenvalues. The morphology of these states is qualitatively discussed in section 4. Section 5 characterizes the system via spectral statistics. Sections 6 and 7 provide a quantitative study of the global and local spatial properties of the WL states. The results are compared to the OPW method in section 8. Finally, the conclusion is given in section 9.

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 singleparticle energy eigenstates in the envelope function approximation, and consider only the in-plane component ψ(x, y) of the wave functions. The QD confinement is treated by a potential V(x, y) which naturally also influences the spatiallyextended states in the WL. The WL is considered as an area of finite size, also called Mesa (see, e.g. [25]), which is chosen to be a region of rectangular shape Ω = [0, Following [14] we describe the confinement in the WL by an infinitely high potential well in 2D.
It has been clearly demonstrated in theoretical [26,27] and experimental [28,29] 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 ω(n + 1) with degree of degeneracy n + 1, where n = 0, 1, 2, . . . , n max . For a single dot the potential is depicted in figure 1. Since the depositioning process in general leads to a spatial distribution of QDs on the WL we have to incorporate the kth QD confinement potential to the overall WL-potential by summing over all N dot non-overlapping QDs where the set of parameters {a k , e k , d k , b k } characterizes the kth QD confinement potential with its minimum at the coordinates (a k , e k ), width b k , and depth d k b 2 k /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 Φ = (1 + √ 5)/2. In our numerics we express the parameters a k , e k , b k , L x , and L y in nm and, using = 1 = M, express energies in nm −2 . First, we study the case of a single QD ( N dot = 1) with three QD states. Three confined states are realistic for several applications such as QD lasing, see, e.g. [30]. 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 cm −2 it is experimentally feasible in a Mesa, as in [25], 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 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 Ω by multiplying L x and L y of system 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 A(2), A (3), and A(4). Moreover, another system containing an ensemble of four identical, non-overlapping In-plane part of a wetting-layer state in a small rectangular Mesa with side lengths L x and L y containing a single quantum dot at (x, y) = (a 1 , e 1 ). The quantum-dot confinement potential V is shown as a surface plot and as contours in the x, y-plane. Outside the quantum dot but inside the Mesa the potential V is zero.
QDs is denoted by B. System C possesses the same parameter settings as system 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 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. [21]. The cut-off energy is determined by the condition that the second subband appears which happens in our case at E cut-off = 3π 2 2 /(2Md 2 WL ) where d WL is the thickness of the WL. Assuming a thin WL with d WL = 0.4 nm we get E cut-off ≈ 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 [17] with, in the context of the present paper, two important differences: (i) in [17] 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 figure 1 there is also a superficial similarity to the Sinai billiard [31], a square billiard with a circular wall located at its center, and the Šeba billiard [32], a rectangular billiard containing a point scatterer (see also the application to ballistic QDs in [33]). 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 [34] 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 [34].

Computation of eigenstates and eigenvalues
We consider the quantum-mechanical eigenvalue problem with the single-particle HamiltonianĤ and the eigenstates ψ i fulfilling zero Dirichlet boundary conditions on the WL boundary ∂Ω. We solve it numerically by expressing Ĥ in a suitable orthonormal basis {ϕ 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 Ω, are discretized, and the coefficients c (i) n = ϕ n |ψ i required for the series expansion of the eigenstates Here, it is convenient to choose the eigenbasis {ϕ n } of Ĥ 0 resulting from removing the QD potential from the Hamiltonian in equation (4). According to the boundary conditions the basis is given by the plane wave-like (checkerboard-like) states with wave numbers k x = n x π/L x , k y = n y π/L y , and n x , n y ∈ 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 ∈ N + we are using a bijective map h : N + × N + → N + introduced by Hopcroft and Ullman [35]. By construction, the matrix elements of the first term in equation (4), recognized as Ĥ 0 , can be evaluated analytically with E nx,ny = 2 k 2 x + k 2 y /2M. The matrix elements of the second term in equation (4) can be written as with B k constituting the region where V k = 0; see figure 1. The numerical evaluation of the double integral (8) is done in local polar coordinates by using scipy.integrate.nquad from the SciPy package [36]. 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 Table 1. Parameter sets of different QD-WL systems. a k , e k , b k , L x , and L y are given in nm. d k is given in nm −4 .

Morphology of energy eigenstates
This section provides a qualitative discussion of the numerically obtained eigenstates. For a quantitative analysis we refer to sections 6 and 7. As representatives we choose systems A(1) and B, see table 1. For the single-QD system A(1) the numerically computed QD states are shown in figure 2. We define the QD states as those states having negative energy; see equations (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 min = −d 1 b 2 1 /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 ψ 8 , ψ 50 , and ψ 1356 display an irregular or nonuniformly spatial structure, some others like ψ 475 , ψ 866 , and ψ 1197 display regular and almost uniform checkerboard-like densities. Occasionally, for system A(1) we also find characteristic states like ψ 429 which might be classified as 'bouncing-ball states' [37].
Furthermore, by taking a closer look at the WL states in figure 3 we notice that some of them exhibit a higher probability in the region of the QD, e.g. ψ 7 , ψ 8 , and ψ 13 . Some others like ψ 11 and ψ 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 section 7.
Whereas the WL states of system A(1) constitute a set of several types, for system B figure 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 higherexcited states looks more irregular. Additionally, due to the QD configuration in system B we come across the opportunity to observe coupled QD states when the inter-dot distance is small enough, like in figures 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.

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 Σ 2 (L) revealing the longrange level correlations, see, e.g. [21,22]. Both quantities are dealing with the fluctuations within a discrete level sequence giving the number of levels up to the energy E. Here Θ denotes the Heaviside step function. N(E) can be segregated into two parts, a smooth part N (E) and a fluctuating one Determining the smooth part N (E) is in general a nontrivial task [38]. 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 [20] 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 A(1) in figure 5. The spectrum is calculated for a Hamiltonian matrix of size 3000 × 3000. Up to the 2100th 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 N (E), the spectrum {E i } can be mapped onto the so-called unfolded spectrum {e i } through e i =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 s ≈ 1. The corresponding probability distribution is called P(s). Uncorrelated energy levels follow the Poisson statistics with nearestneighbor spacing distribution 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 [39] 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 figure 6(a) for  Assignments of colors are the same as in figure 2 and the QD is indicated as a circle. x and y are measured in nm.
system 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. [22]) A better agreement for system A(1) is achieved with the socalled Semi-Poisson statistics [40,41] with nearest-neighbor spacing distribution   Like the Wigner surmise P 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 [41] and appears, e.g. in pseudo integrable systems [40,42,43] and multi-walled carbon nanotubes [44]. Furthermore, it has also been controversially discussed in the context of the Šeba billiard [45].
Next, we check the Brody distribution [46] and the gamma function Γ. The Brody distribution is a heuristic interpolation between Wigner surmise (β = 1) and Poisson statistics (β = 0), see [21]. Figure 6(a) shows a very good agreement with the numerical data for β = 0.35, which is even slightly better than for the Semi-Poisson statistics. We have also tested the Berry-Robnik distribution [47] P BR (s, β) with the complementary error function erfc and α = 1 − β. This distribution also interpolates between Wigner surmise (β = 1) and Poisson statistics (β = 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 1 (not shown). This failure of the Berry-Robnik distribution is consistent with other numerical studies [48,49]. It results from neglecting dynamical tunneling between regular and chaotic regions [50].
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 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 figure 6(b). It can also be seen that in this case the Semi-Poisson statistics gives a slightly better agreement.
Following the discussion in [44,45] 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 A(1) gets different contributions from separate energy windows as shown in figures 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 β = 0.5. For increasing energies the distribution is slightly better described by Brody distributions with β 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 B with altogether twelve QD states using a 4000 × 4000 Hamiltonian matrix. For low energies, the nearest-neighbor spacing distributions (figure 7(d)) can be characterized by the Brody distribution with β = 1 which equals the Wigner surmise. By increasing the energy a slow transition to β = 0.9 (figure 7(e)) and β = 0.6 (figure 7(f)) is observed. In the high-energy regime also the Semi-Poisson statistics describes the numerical data very well. Note that system B is significantly closer to the statistics of random-matrix theory than system 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 section 4, see figures 3 and 4.
To study long-range level correlations we consider the number variance and for the GOE (see, e.g. [22]) with the Euler constant γ ≈ 0.577 22. In the case of the single-QD system A(1) figure 8 shows a remarkable agreement with the number variance of the Semi-Poisson statistics in equation (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. [52]. For system A(3) with one QD on a three times larger WL, figure 8 shows a shift towards the GOE in equation (20). For the 4-QD system B we find that Σ 2 (L) is even closer to Σ 2 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.

Amplitude distribution
After the qualitative discussion on the morphology of the WL states in section 4, we now provide quantitative measures for the observed complexity. According to Berry [53], the (realvalued) amplitude ψ 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 [54] To determine P(ψ) for a given WL state we evaluate its wave function ψ(x, y) at 128 × 128 points in the x, y-plane. The outermost points of the WL region Ω are excluded to avoid an artificial peak at ψ = 0 due to the zero Dirichlet boundary conditions. Then P(ψ) is generated as a histogram with 90 bins in [−0.5, . . . , 0.5], normalized to unity. For the single-QD system A(1) figure 9(a) shows two examples, ψ 503 and ψ 1356 , which are fairly well described by equation (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 equation (21) (not shown), for instance the bouncing-ball state ψ 429 and the checkerboard-like states ψ 866 , ψ 1197 ; see figure 3. It is well known that such regular states are not described by a Gaussian random variable [54,55].
For the 4-QD system B we find that nearly all among twenty randomly chosen sample states between ψ 492 , . . . , ψ 1581 show a remarkably good agreement with equation (21). Two examples are shown in figure 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.

Localization
In this section we take up a question raised in section 4 about the WL states' tendency to avoid or seek the region of the QD(s). For this purpose, we define the localization of the ith eigenstate as where, again, B k is the region where the confinement potential V k of the kth QD is nonzero. These double integrals are summed over all QDs and, like equation (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 • ) and the area of the WL (A) which can be written as For system A(1) we compute the localization W i for the WL states with i = 4, . . . , 1000. We consider here a smaller amount of states than shown in figure 3 due to the high sensitivity of the W i on the convergence of the series expansion (equation (5)), which here again contains 3000 terms. For comparison, we calculate W i also for system A(1) under withdrawal of the QD having the plane wave-like energy eigenstates in equation (6). The resulting probability distributions P(W) are depicted in figure 10. In the system without QD the distribution essentially shows one sharp peak at W ≈ 0.03 ≈ A • /A (see equation (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 W which is therefore again very close to the area ratio in equation (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 γ i of the ith eigenstate would be in first-order perturbation theory for weak coupling proportional to the localization W i . Hence, γ i / γ = W i / W . It has been shown that for weakly open systems (see, e.g. [56]) the lifetime statistics is given by a χ 2 distribution with normalized decay rate z = γ i / γ and number of decay channels ν . Equation (24) is the many-channel generalization of the Porter-Thomas distribution [57]. The χ 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 [58] and references therein. We therefore assume P(W) ≈ χ 2 ν (W/ W ) and use the number of decay channels ν as a fitting parameter. An upper bound for ν 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 (λ/2) 2 with the wavelength λ corresponding to the considered energy. Second, we count how many of such areas fit into the area A • of the QD region. With the largest (unfolded) energy e max used in the statistics and equation (23) we find It can be seen in figure 10 that the numerically obtained P(W) of A(1) is rather well fitted by the χ 2 distribution for ν = 12. This is consistent with the upper bound ν max ≈ 39 using e max = 1000.
The system A(1) without QD, see figure 10, cannot be well fitted by a χ 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 B. Here, the agreement with the χ 2 distribution is even slightly better than for system A(1). The fitted value of ν = 21 is consistent with the upper bound ν max ≈ 30, using e max = 1000 in equation (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 equation (22) for the single-QD system A(1) to be a spatially shifted region of the same shape and size, still in Ω, 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 ν = 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.

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 [14], we consider the single-QD system 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 with C = (Mω/π ) 1/2 . In equation (26) the parameter d 1 is incorporated via the frequency ω = 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 ϕ dot : It should be noted that in a set {u i } containing u 1 = ϕ 1 the QD state ϕ dot could also constitute an additional state.
The approximated states are determined as  x and y are measured in nm.
where i = 2, . . . , N and N (OPW) j with j = 1, . . . , N is a normalization constant of the j th state, such that the state is normalized to unity in Ω. But unlike the exact eigenstates the OPWs ψ (OPW) i are just orthogonal to the ground state ψ but not necessarily to each other. Two examples of OPWs are shown in figures 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 figures 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 figures 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 [14], 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 ψ (OPW) 1 |u i in equation (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 [17].

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 shortand 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 χ 2 distribution.
The observed localization behavior can be important for capture time calculations, such as in [14]. 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 [59] 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. [21]. 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. [60], 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.