Electric field induced tuning of electronic correlation in weakly confining quantum dots

We conduct a combined experimental and theoretical study of the quantum-confined Stark effect in GaAs/AlGaAs quantum dots obtained with the local droplet etching method. In the experiment, we probe the permanent electric dipole and polarizability of neutral and positively charged excitons weakly confined in GaAs quantum dots by measuring their light emission under the influence of a variable electric field applied along the growth direction. Calculations based on the configuration-interaction method show excellent quantitative agreement with the experiment and allow us to elucidate the role of Coulomb interactions among the confined particles and -- even more importantly -- of electronic correlation effects on the Stark shifts. Moreover, we show how the electric field alters properties such as built-in dipole, binding energy, and heavy-light hole mixing of multiparticle complexes in weakly confining systems, underlining the deficiencies of commonly used models for the quantum-confined Stark effect.


I. INTRODUCTION
Quantum optoelectronic devices capable of deterministically generating single photons and entangled photonpairs on demand, are considered key components for quantum photonics. Of the different available systems, semiconductor quantum dots (QDs) are one of the most promising candidates, because they combine excellent optical properties with the compatibility with semiconductor processing and the potential for scalability. 1-8 A prominent example is represented by GaAs/AlGaAs QDs fabricated by the local droplet etching (LDE) method [9][10][11][12][13][14] via molecular beam epitaxy (MBE). These QDs can show ultra-small excitonic fine-structure-splitting (FSS), with average values of ≈ 4 µeV 14,15 , ultra-low multi-photon emission probabilities, with g (2) (0) below 10 −4 , 16 , stateof-the-art photon indistinguishabilities 17 and near-unity entanglement fidelities of 0.978 (5) 18 . Devices based on LDE GaAs QDs have recently achieved high performance as sources of polarization-entangled photon pairs 18,19 , which led to the demonstration of entanglement swapping 20,21 and quantum key distribution 22,23 .
In addition to their excellent optical properties, semiconductor QDs also provide a platform for photon-tospin conversion 24,25 , building up bridges between photonic and spin qubits 26 . In addition, the nuclear spins of the atoms building up a QD are emerging as longlived quantum storage and processing units that can be interfaced to photons via coupled electron spins 27,28 . To efficiently initialize and manipulate single spins confined in QDs, the QD layer is typically embedded in a diode structure, which allows the charge state to be deterministically controlled 29 . By tuning the diode bias, not only is the charge state modified, but the magnitude of the electric field (F d ) along the QD growth direction is as well. In turn, F d modifies the energy and spatial distribution of the confined single particle (SP) states as well as the Coulomb and exchange interactions among the charge carriers via the so-called quantum-confined Stark effect (QCSE), leading to deep changes in the electronic and optical properties of the QDs [30][31][32][33] . Therefore, a fundamental understanding of the effects of F d in this kind of quasi-zero dimensional structures is highly desirable. LDE GaAs QDs formed by filling Al (or Ga) dropletetched nanoholes (NHs) at high substrate temperature (∼ 600 − 650°C) present advantages over conventional strained QDs and QDs obtained by droplet epitaxy. These advantages include negligible strain, minimized intermixing of core and barrier material, a low QD density of ≈0.1 µm −2 , high ensemble homogeneity, and high crystal quality, [9][10][11]34 thus providing a particularly clean and favorable platform for both fundamental investigations and applications of QCSE. To the best of our knowledge, only a few works have been dealing with the physics of GaAs QDs in externally applied electric fields. 29,[35][36][37][38][39][40] As an example, Marcet et al. 39 and Ghali et al. 40 used vertical fields (perpendicular to the growth plane) to modify the FSS of neutral excitons confined in natural GaAs QDs (thickness or alloy fluctuations in thin quantum wells, with poorly defined density, shape and optical properties). Besides that, several simulation models based on SP assumption were also built up to explain the arXiv:2105.11244v4 [cond-mat.mes-hall] 14 Feb 2022 charge noise (emission line broadening caused by fluctuating electric field around the QDs produced by charge trapping/detrapping occurring at random places). 41 Nevertheless, those models neither fully explain the behavior of the charge carriers in the electric field, nor take into account correlation effects 35 completely. On the contrary, we note that correlation is of particular importance in the GaAs/AlGaAs QD system because of the generally large size of the studied QDs. [42][43][44] For example, without including the effects of correlation, the binding energy of X + with respect to X 0 shall be rather small and attain negative values (anti-binding state) rather than positive ones (binding state), 45 which is in contrast with the experimental observations. [46][47][48] Although positive binding energies have been theoretically calculated for GaAs QDs obtained by "hierarchical self-assembly", 43 quantitative agreement between theory and experiment has not been demonstrated so far. In addition, detailed studies of the electric field effects on the Coulomb interactions between electrons (e − ) and holes (h + ) in GaAs QDs are still lacking.
In this work, we conduct a combined experimental and theoretical study of the QCSE in individual GaAs QDs. Our experiments, based on micro-photoluminescence (µ-PL) spectroscopy, offer direct information on the permanent electric dipole moment (p) and polarizability (β) of the neutral exciton X 0 (X 0 ≡ 1e − + 1h + ) and X + (X + ≡ 1e − + 2h + ) states in GaAs QDs, which sensitively depend on carrier interactions in those nanostructures. In the experiment, we are able to tune the QD emission energy over a spectral range as large as 24 meV thanks to the large band offsets between QD material (GaAs) and surrounding Al 0.4 Ga 0.6 As barriers. Such "giant Stark effect" 31 allows us to observe a crossing of the X + emission line with that of the X 0 with increasing F d , see also Appendix I.. The evolution from a binding to an anti-binding X + state (relative to X 0 ) indicates substantial electric-field-induced changes in Coulomb interactions and possibly correlation. The calculations of the aforementioned complexes are performed using the configuration-interaction (CI) method, 49-52 see also Appendix II., with SP basis states obtained using the eightband k·p method computed with the inclusion of the full elastic strain tensor and piezoelectricity (up to second order 53,54 ) by Nextnano 55 software package. Our computational approach provides consistent results with all experimental data. These calculations not only extend the investigated F d 's to the range inaccessible in the experiments and explore different QD morphologies but also maps the behavior of the corresponding direct Coulomb integrals (electron-hole J eh , hole-hole J hh ) and valence band mixing as F d is varied. Interestingly, we find that the often overlooked correlation effects among e − and h + plays a central role for describing the QCSE and that the commonly assumed quadratic dependence of the emission energy shift on F d in QDs is questionable.

II. QUANTUM-CONFINED STARK EFFECT IN A SINGLE GaAs QD
We start by measuring the Stark shifts of X 0 and X + states of GaAs QDs by µ-PL spectroscopy. The shape of the QD is defined by the Al-droplet-etched NH [see Fig. 1 (a)], with a depth of ∼ 7 nm, a full width at half maximum depth of ∼ 33 nm), and ∼ 1 − 2 nm thick "wetting layer" (WL) above the NHs formed by the GaAs filling. 14,15 To apply an electric field F d along the growth direction, the QDs were embedded in the intrinsic region of a p-i-n diode structure (see the details in Appendix I) as sketched in Fig. 1 (b). The direction of F d and the corresponding movement of the e − (h + ) wavefunction is marked in Fig. 1 (c). F d is calculated as Figure 1 (d) shows typical µ-PL spectra obtained from a QD (marked as QD1) as a function of F d . Near F d = 0, an isolated X 0 transition is found at 1.611407(2) eV, accompanied by multiexciton states at lower energies (1.60843 − 1.60381 eV). This configuration agrees qualitatively with other reports on GaAs QDs grown by LDE, 29,45,48 droplet epitaxy 56 and hierarchical, self-assembly, 42,43 and it is different from that observed in InGaAs QDs, for which X + usually attains higher energy, and X − attains lower energy compared to X 0 . 32,57,58 The X 0 state was identified by the polarization and power mapping. The other charged complexes can be calibrated by combining power mapping and temperaturedependent µ-PL measurement, as shown in our previous work ? . Here we would like to focus only on the most intensive X + , as it has minimal interaction (mixing) with other charged states. That X + was paired to the X 0 by the position check. Our sample has an ultra-low QD density (0.3 − 0.4 QD/µm 2 ), allowing single QD excitation. Energy shifts for F d 30 kV/cm are not observed in our experiments because of the current injection in the diode. Investigations on the electroluminescence (EL) of this type of device have been reported previously in Ref. 59. The XX transition is usually not recognizable under above-band excitation (except for some values of F d ) due to the fact that it competes with other charged states. At large F d (F d 240 kV/cm) the µ−PL signal becomes faint and cannot be tracked because of the field ionization of excitons. 58 Overall, the emission energy is red-shifted by almost 24 meV upon increasing F d . We extract the energy of X 0 and X + by performing Gaussian fitting of their µ−PL spectra for the corresponding F d , and we plot those for QD1 in Fig. 2 (a) along with the data for another QD (marked as QD2). In both cases we observe a smaller energy shift for X + compared to X 0 , leading to a crossing for sufficiently large values of F d .
In the simulation we have modeled the NH as a   cone with the basal diameter of 40 nm, a height (h) of 4 − 9.5 nm and a wetting layer thickness of 2 nm. Note, that later on we also provide the theory result for lensshaped dots with the same basal diameter as reference cone-shaped dots. The lens shape, although it does not reproduce the real NH shape, has an increasing lateral space for taller QDs. In the experiments, the taller (larger) QDs will also be "wider" than the short (smaller) one. The simulated Stark shifts of the QDs are plotted together with the experimental data from 5 dots in Fig. 2 (b). Calculation results are also shown for F d < 0, which is however not experimentally accessible with the present diode structure. It is interesting to note that the parabolic shifts are not symmetric around F d = 0, as already predicted in Ref. 35. Concomitantly, the maximum of the emission energy appears at F d > 0. Both effects are the result of the asymmetric shape of the QDs along the F d direction, i.e., the z-axis combined with the different behaviors of e − and h + as their wave functions move along the z-axis, thus, experiencing different lateral confinements. On the other hand, the maximum of emission energy at non-zero F d can be interpreted with the existence of a permanent electric dipole, which we will discuss in the following section.

III. PERMANENT ELECTRIC DIPOLE MOMENTS AND POLARIZABILITY OF NEUTRAL AND POSITIVELY CHARGED EXCITONS
The shifts of the X 0 and X + energy induced by F d are commonly described by the following quadratic equation: where E 0 is the emission energy for F d = 0, and p z and β can be intuitively interpreted as the permanent electric dipole moment and polarizability of the corresponding complexes, respectively. 33,58,60,61 The quantity p z /e can be seen as the distance between the electron and hole probability densities along the z-axis. The results for QD1 and QD2 fitted by Eq. (1) for F d in the range 30 kV/cm < F d < 240 kV/cm are shown in Fig. 2 (a) and Table I. Data for X + at F d < 120 kV/cm were excluded as we could not unequivocally identify the X + band in that region. The same was done for data obtained from the other three QDs (marked as QD3-QD5 in Figures 2-4 and Table I) and the fit is performed in the F d range of ≈ 100 − 250 kV/cm. Figure 3 (a) summarizes the fitted values of p z /e for X 0 and X + for five QDs. The negative values of p z /e for X 0 (see Table I) indicate that the e − wavefunction is shifted closer to the bottom of the NH (tip of the dot) compared with h + for F d = 0, as sketched in the bottom inset of Fig. 3 (a). The corresponding positions of the e − /h + wavefunction and the p z /e value (p z /e=−0.39 and −0.31 nm for QD1 and QD2) are close to the experimental data reported in Ref. 62  We start evaluating our theoretical results for X 0 or X + given in Fig. 2 (b) by performing the same fitting procedure using Eq. (1) as for experiment. However, we find that the values of p z /e obtained using that pro-  , for a similar F d range as for experiment, we find p z /e < 0, in agreement with experimental data (for comparison of fits see Fig. 6 in Appendix IV). Thus, the aforementioned way to obtain the value of permanent electric dipole moments is unsatisfactory. It actually points to the fact that the evolution of energy of QD multi-particle complexes does not follow equation (1) faithfully. In order to access the intrinsic distance p z /e in GaAs QDs, we can use directly the SP h + and e − states, similarly to Refs. 33 and 63. However, this approach is reasonable only when the e −h + distance is evaluated between the SP ground states of those quasiparticles. Thus, this option is available only for X 0 (not X + or any complex consisting of more than two particles) and for systems that can be reasonably well described in the single-particle picture, which is not the case for GaAs/AlGaAs QDs where already X 0 is sizeably influenced by correlation. 44 Hence, instead we develop a method of obtaining p/e directly during our CI calculations 51 as

and the simulated result
where η l m is an m-th element of the l-th CI matrix eigenvector |M l = η l 1 , . . . , η l nSD T corresponding to m-th wherer h (r e ) marks the position operator of h + (e − ) SP eigenstate |Ψ h k ( Ψ ej ), the indices j and k mark the SP states included in SD m , and the bra-ket integrals are evaluated over the whole simulation space. Note, that in Eq. (2) the CI eigenstates η l m are used as "weights" of the expectation values computed from SP states. Thus, it provides a rather general way of including the effect of correlation to the "classical" properties related to SP states. Note that the method is partly motivated by our previous results in Ref. 44.
We show the p z /e component of Eq.
(2) in Fig. 3 (a) for X 0 and X + . The small computed values of p z /ethat can be expected also from the probability density plots in Fig. 3 (c)) (see also Appendix III.) -are plotted together with the values (also negative) extracted by fitting the experimental data with Eq. (1). The calculations indicate that the permanent electric dipole of excitons confined in GaAs QDs is very small. This is very different from the situation typically encountered in strained QDs, where the dipole is mostly determined by opposite effects, namely the alloy gradient and the strain inhomogeneities combined with piezoelectricity. 33,60,[64][65][66][67][68][69][70] In view of the minuscule values of p z /e that we find in both experiment and theory it is reasonable to discard the p z /e term in fitting using Eq. (1) in the case of our data; see also the comparison of the fitting with/without a linear term in Eq. (1) in fig. 2 (a), as |p z /e| is in atomic scale In contrast to p z /e, we find for β/e of X 0 (X + ) a more consistent agreement of fits by Eq. (1) between theory and experiment, see Fig. 3 (b). The results of the fits for different intervals of F d are again given in Appendix II. Furthermore, β/e of X 0 (X + ) shows a clear dependence on E 0 . The larger QDs, with smaller E 0 , tend to have a larger magnitude of β X 0 (β X + ) for X 0 (X + ), consistent with the results reported in Ref. 62. The theoretical prediction in Ref. 41 and 65 also pointed out that with a fixed shape and chemical composition profile, β is mostly sensitive to the QD height. A taller QD provides in fact more room along the z-direction for the confined e − -h + pairs to move away from each other when pulled apart by F d , resulting in a stronger red-shift in spite of the reduced e − -h + binding energy.
We will discuss the detailed role of e − − h + Coulomb interaction and correlation in the Stark shift with the help of simulation in the following section.

IV. TRION BINDING ENERGY AND THE ROLE OF COULOMB INTEGRALS IN ELECTRIC FIELD
To describe the evolution of the relative binding energy E b = E(X 0 ) -E(X + ) with F d we assume a quadratic dependence as in Eq. (1) with an omitted linear term (see above discussion) where E b,0 marks E b for F d = 0. Thereafter, using Eq. (4) we fit the difference between E(X 0 ) and E(X + ) taken from corresponding dependencies in Fig. 2 (b) and we obtain the parameters E b,0 and β * E b , which we show alongside the calculated values in Fig. 4 (a) and (b), respectively. From Fig. 4 (a) we see that the calculated E b,0 is satisfyingly close to the experimental data for both the cone-and the lens-shaped dots, in contrast to former CI calculations. 43 Remarkably, a positive trion binding energy as large as large as 5 meV is obtained from realistic calculations. The E b,0 values are also close to those reported in Ref. 71 (= E b,0 linearly increasing from ∼ 2.4 to ∼ 2.9 meV for emission energies increasing from ∼ 1.56 to ∼ 1.61 eV). We ascribe the agreement between our theory and experiment to an almost full inclusion of the correlation effects, which will also be discussed and tested in the following.
However, we first show that the physical reason for the disagreement of Eq. (1) with theory is due to the omission of the effect of correlation in Eq. (1) as well. We start by writing the energies of the final photon states after recombination of X 0 and X + as 44,72 where E X + is the energy of X + before recombination, and J eh,X 0 , J eh,X + , and J hh are the Coulomb interactions of e − -h + pairs in X 0 and X + , and of the h + -h + pair, respectively; ε e (ε h ) is the single particle e − (h + ) energy, and δ(X 0 ) (δ(X + )) marks the energy change due to the effect of correlation for X 0 (X + ). Consequently, the E b can be written as: where δ = δ(X 0 ) − δ(X + ). Note that we have completely neglected the exchange interaction for elaborating the simplified model in Eq. (7) since we found that to be ≈ 100 times smaller than direct Coulomb interaction in our CI calculations (for which the exchange interaction was of course not neglected). In Fig. 4 (d) we plot β * E b /e for J eh , J eh − J hh , J hh , and E b,sim from simulation on E 0 . Note, that β * E b /e values were obtained by fits using Eq. (4) of the theory dependencies of J eh , J eh − J hh , J hh , and E b,sim on F d computed by CI with a 12×12 SP basis, for the fits see Appendix V. Clearly, we find that β * E b /e depends on the QD size. For bigger QDs (smaller E 0 ), with steeper side facets and larger height, |β * E b /e| of J eh is more pronounced compared to that in flatter QDs. The reason is that taller QDs facilitate the e − -h + separation (polarization) under the influence of vertical F d . On the other hand, |β * E b | for J hh is smaller in larger QDs. The reason is that larger QDs allow the separation between h + to be larger, thus reducing the Coulomb repulsion. Since the value of |β * E b | for J hh is smaller than that of J eh for every QD, β * E b for E b,sim has a larger contribution of that corresponding to J eh . However, we notice that |β * E b | for J eh − J hh is still smaller than that of E b,sim (see the corresponding curves in Fig. 4 (c)). That means, besides J eh and J hh there must be another important variable in Eq. (7) changing with F d . Therefore, the last component in Eq. (7), i.e., the correlation effect δ, must also vary with F d , i.e., δ = δ(F d ).
To prove the importance of the correlation effect in our system, we calculated E b based on the CI model for the simulation with increasing SP basis from two e − and two h + (2×2) states to twenty-four e − and twenty-four h + (24×24) states. The result is plotted in Fig 4 (c). Clearly, in the absence of correlation, i.e., using 2×2 and 4×4 basis, X + is anti-binding with respect to X 0 , in contradiction with the experiment. However, with increasing basis size, the effect of correlation gains importance and X + becomes binding with respect to X 0 . The increase  Fig. 2 where blue corresponds to h = 4 nm and red to h = 9.5 nm. The Bloch state contents for both X 0 and X + were calculated using the CI model with the basis consisting of 12 SP e − and 12 SP h + state, with the effects of the direct and the exchange Coulomb interaction, and the correlation effect being included, see also Appendix III.
of E b is steep up to 12×12 basis, where it almost saturates. Note, that the dependence was computed for the largest considered QD, i.e., h = 9.5 nm, where the effect of correlation was expected to be the most significant.

V. VALENCE BAND MIXING OF THE NEUTRAL EXCITON AND THE POSITIVE TRION
In this section we study the effect of F d on heavy-(|HH ), light-(|LH ), and spin-orbit (|SO ) hole Bloch state mixing for X 0 and X + ground states. The corresponding contents divided by the sum of those components, i.e., κ(HH) + κ(LH) + κ(SO) where κ marks the respective content, is shown in Fig. 5.
Note that the method of extracting the Bloch band content of CI states we show in Appendix III. (see also Ref. 44) and the conversion between {|HH , |LH , |SO } and {|p x , |p y , |p z } bases is provided in Appendix VI.
We observe asymmetric dependencies around F d = 0. The content of |HH increases with F d with a concomitant decrease in the contribution of |LH states. Since the holes are pushed towards the bottom of the QD by positive F d (Fig. 3 (c)), the h + SP state barely feels the broken translation symmetry along z-axis, since the lateral confinement is weaker at the bottom of the QD. Without broken symmetry the hole states tend not to mix, which causes an increase of the amount of |HH Bloch states. On the other hand, negative F d (F d applied along the opposite direction) pushes the holes towards the top of QD, thus, increasing the valence-band mixing (increase of the content of |LH and |SO Bloch states). According to Appendix VI while |HH Bloch states are purely |p x and |p y -like, |LH and |SO Bloch states consist also of a non-negligible amount of |p z states. However, for the |SO states, the same amount of |p x , |p y , and |p z Bloch states is involved, which leads to a more symmetric trend than in the case of |HH and |LH states.
Interestingly, for negative F d , the content of |HH states changes the trend after an initial decrease for F d values close to zero and starts to grow again for F d < F d,crit , which is dependent on the QD height. Note, that this change is more pronounced for X 0 . Since the contents of |HH , |LH , and |SO are normalized to the total sum of all valence band components, we can directly compare X 0 and X + . In the case of X + the direct and exchange Coulomb interaction between e − and h + is twice as large as that for X 0 . Also the direct and exchange Coulomb interaction between two holes is included and the correlation affects the complexes in a different way; see Eqs. (5) and Eq. (6). As one can see, the aforementioned effects influence valence-band mixing rather strongly.
Now we focus on the dot size dependence of the contents of |HH and |LH Bloch states. For F d < 50 kV/cm (F d < 125 kV/cm) for X 0 (X + ), the amount of |HH (|LH ) Bloch states decreases (increases) with increasing height of the dot, as smaller QDs display larger energy separation between confined |HH and |LH SP states. Since the variation of valence band mixing is observed to be more pronounced in larger QDs (increased height), we observe the crossing of the HH curves for F d = 50 kV/cm (F d = 125 kV/cm) in case of X 0 (X + ). Thereafter, for F d > 50 kV/cm (F d > 125 kV/cm) for X 0 (X + ), the trend of the size dependence is reversed, i.e., bigger QDs have a larger amount of |HH states than QDs with smaller height. For such large fields the dominant part of the SP hole wavefunction leaks into the wetting layer and laterally delocalizes, leading to the a faster increase of the content of |HH states. We assume that for the same F d (F d > 50 kV/cm for X 0 ), all wavefunctions leak into the wetting layer with the same amount of probability density. Hence, the wavefunctions, with larger volume, i.e., for bigger QDs, consist of more |p x and |p y Bloch states and so also the larger contribution of |HH .

VI. CONCLUSIONS
In summary, by conducting detailed µ-PL spectroscopy measurements of the emission from LDE-grown GaAs/AlGaAs QDs modulated by an externally applied electric field and in conjunction with conscientious calculations of multiparticle states, we reveal the influence of the electric field on the Coulomb interaction among charge carriers in GaAs QD. The experimental data and the configuration interaction calculation clearly show the dot size dependence of the polarizability of X 0 and X + . Thorough analysis of configuration interaction calculations sheds light on the deficiencies of the commonly used analysis of the quantum confined Stark effect by highlighting the striking effect of correlation and the direct Coulomb interaction energy between holes, which change with applied field and which are also significantly influenced by the asymmetry of the QD along the field direction, especially in large quantum dots. Moreover, we analyzed the Bloch state composition of exciton and trion complexes as a function of applied electric field, and we emphasize the influence of QD height as well. Finally, we note that our multiparticle simulation model based on the full configuration-interaction approach with large number of single-particle basis states provides excellent quantitative agreement with the experiment, and proves the non-negligible role of the correlation effect on the Stark shift for the nanosystems.
In the experiments, the QDs were embedded in the intrinsic region of a p-i-n diode structure (thickness of 95 nm-124 nm-170 nm). The thickness of the diode and the location of the QDs were chosen to obtain a simple Au-semiconductor-air planar cavity after transfer on an Au-coated substrate to enhance the out-coupling efficiency (see the details in Ref. 59). Note that, minor bi-axial strain can be introduced during processing.
The FSS of X 0 s from this sample is ∼ 12 − 15 µeV near zero-field and increases slightly to ∼ 20 µeV at the maximally available field due to a slight in-plane asymmetry. The linewidth of one single component of X 0 is ∼ 40 µeV . The X 0 energy is chosen to be the average of the two components. We've tested the consequence of choosing different polarization components. The result showed that this ±10 µeV tuning has a negligible effect (<< than the uncertainty) on the fitting results of β and E 0 , since ±10 µeV is a quarter of X 0 linewidth and two magnitudes less than the energy difference between different dots.

APPENDIX II.
For better readability we reproduce in Appendix II and III the description of our CI method 51 , given previously also in 44 . Let us consider the excitonic complex |M consisting of N e electrons and N h holes. The CI method uses as a basis the Slater determinants (SDs) consisting of n e SP electron and n h SP hole states which we compute using the envelope function method based on k · p approximation using the Nextnano++ simulation suite 55 . SP states obtained from that read Ψ ai (r) = ν∈{s,x,y,z}⊗{↑,↓} where u Γ ν is the Bloch wave-function of an s-like conduction band or a p-like valence band at the center of the Brillouin zone, ↑/↓ mark the spin, and χ ai,ν is the envelope function, where a i ∈ {e i , h i }.
The trial function of the excitonic complex then reads where n SD is the number of SDs |D M m , and η m is the constant that is looked for using the variational method. The m-th SD can be found as In Eq. (12) q i and q j label the elementary charge |e| of either electron (−e), or hole (e), and (r, r ) is the spatially dependent dielectric function. Note, that the Coulomb interaction is treated as a perturbation. The evaluation of the sixfold integral in Eq. (12) is performed using the Green's function method [50][51][52]73 ∇ (r)∇Û ajl (r) = 4πe 2 0 Ψ * aj (r)Ψ al (r), where a, b ∈ {e, h} and ∇ := ∂ ∂x , ∂ ∂y , ∂ ∂z T . Finally, note that (r, r ) was set to bulk values 44,52 for the CI calculations presented here.

APPENDIX III.
To visualize the contents of SP states computed in multi-particle complexes calculated by CI, we need to transform the results of CI calculations to the basis of SP states instead of that of SDs. 44 During the set-up of SDs within our CI algorithm, we create the matrixÂ with rank n SD × N , where m-th row consists of SP states used in the corresponding SD A m = Ψ ej , . . . , Ψ e j+Ne −1 ; Ψ h k , . . . , Ψ h k+N h −1 . (14) Further, resulting from diagonalization of the CI matrix, we get n SD eigenvectors with n SD components |M l = η l 1 , . . . , η l nSD T , where the index l identifies the eigenvector. We choose those values of η l m that correspond to the A m consisting of a particular SP state Ψ ej {Ψ h k }, we sum the squares of the absolute values and we obtain the vector The values c ej and c h k are then normalized by imposing that j c l ej + k c l h k = 1. Since |η l m | 2 describes the weight of the corresponding SD in the CI eigenvector, we look for the weights of individual SP electron or hole states.
The procedure described thus far allows us to study also other excitonic properties, such as the influence of multi-particle effects on band mixing or visualizing the probability density of the studied excitonic complexes.
For visualizing the probability density of an eigenstate of the complex |M l with wave-function Φ l M (r) as in Fig. 3 (c), we calculate Finally, the probability density is finally normalized, i.e., M l |M l = 1.
In the case of band mixing we multiply the contents of {|S , |HH , |LH , |SO } of the particular SP state by the corresponding coefficient from Eq. (18). Hence, we get the matrix with rank (n e + n h ) × 4 for each l and we sum separately all |S , |HH , |LH and |SO contents in that matrix to get the four corresponding values for each CI state. Again, we normalize the contents in the same fashion as for Eq. (18). The aforementioned procedure was used to obtain the results shown in Fig. 5. In Fig. VII, we present the permanent electric dipole moments (p z ) and the polarizability (β) plotted as a function of the zero field energy E 0 of the corresponding complex X 0 or X + .