Valley-selective energy transfer between quantum dots in atomically thin semiconductors

In monolayers of transition metal dichalcogenides the nonlocal nature of the effective dielectric screening leads to large binding energies of excitons. Additional lateral confinement gives rise to exciton localization in quantum dots. By assuming parabolic confinement for both the electron and the hole, we derive model wave functions for the relative and the center-of-mass motions of electron–hole pairs, and investigate theoretically resonant energy transfer among excitons localized in two neighboring quantum dots. We quantify the probability of energy transfer for a direct-gap transition by assuming that the interaction between two quantum dots is described by a Coulomb potential, which allows us to include all relevant multipole terms of the interaction. We demonstrate the structural control of the valley-selective energy transfer between quantum dots.


Valley-selective energy transfer between quantum dots in atomically thin semiconductors
Anvar S. Baimuratov 1* & Alexander Högele 1,2 In monolayers of transition metal dichalcogenides the nonlocal nature of the effective dielectric screening leads to large binding energies of excitons. Additional lateral confinement gives rise to exciton localization in quantum dots. By assuming parabolic confinement for both the electron and the hole, we derive model wave functions for the relative and the center-of-mass motions of electronhole pairs, and investigate theoretically resonant energy transfer among excitons localized in two neighboring quantum dots. We quantify the probability of energy transfer for a direct-gap transition by assuming that the interaction between two quantum dots is described by a Coulomb potential, which allows us to include all relevant multipole terms of the interaction. We demonstrate the structural control of the valley-selective energy transfer between quantum dots.
The unique properties of two-dimensional (2D) materials provide versatile opportunities for nanomaterial physics 1 . Within this realm, monolayers of transition metal dichalcogenides (TMD) represent 2D crystalline semiconductors with unique spin and valley physics for opto-valleytronic applications [2][3][4] . Coulomb electron-hole attraction and the nonlocal nature of the effective dielectric screening lead to a large binding energy of excitons, which dominate both light absorption and emission [5][6][7][8] . The combination of exceptional brightness and spin-valley coupling opens up novel opportunities for tunable quantum light emitters for quantum information processing and sensing 9 realized on the basis of excitons confined in TMD-based systems.
There are various approaches to realize exciton confinement in TMD monolayers. Impurities, vacancies, or strain in monolayer TMD crystals, as well as local modulations of the immediate environment, modify the energy gap of 2D materials 10 . By providing additional lateral confinement, local disorder is known to confine excitons in a relatively small area of TMD monolayers and give rise to quantum dot (QD) excitons. Spectral signatures of quantum dot exciton localization with bright and stable single-photon emission were observed from unintentional defects in monolayer tungsten diselenide [11][12][13][14][15] . Subsequently, strain engineering has proven as a viable deterministic approach to obtain spatially and spectrally isolated quantum emitters in monolayer and bilayer TMDs [16][17][18][19] , and controlled positioning has been achieved by irradiating monolayer crystals with a sub-nm focused helium ion beam 20 . Alternatively, QDs have been realized by electrostatic confinement [21][22][23] , or by creating lateral TMD heterostructures forming a potential well 24 . In vertical TMD heterostructures, moiré superlattices give rise to periodic QD arrays hosting localized excitons 25,26 .
By preserving strong spin-valley coupling, TMD QDs inherit optovalleytronic properties from their 2D host system 27 , as the intervalley coupling is weak due to the vanishing amplitude of the electron wave function at the QD boundary and hence valley hybridization is quenched by the much stronger spin-valley coupling 28 . As in conventional QDs, the oscillator strength and radiative lifetime of confined excitons are strongly size-dependent, which results in oscillator strength enhancement and ultrafast radiative annihilation of excitons, varying from a few tens of femtoseconds to a few picoseconds 29 . In the presence of an external magnetic field bound states in TMD QDs can be considered as quantum bits for potential applications in quantum technologies [30][31][32] .
In this work, we study nonradiative resonance energy transfer between two adjacent QDs in TMD monolayers 33 . Building on theories initially developed for molecular systems by Förster in the framework of dipole-dipole interaction 34 and generalized by Dexter for quadrupole and exchange interactions 35 , we derive the theory of nonradiative resonance energy transfer for atomically thin QDs hosted by 2D crystals. For conventional QDs with sizes in order of tens of nm, the multipole nature of Coulomb interactions and energy transfer through dipole-forbidden states must be taken into account 36,37  www.nature.com/scientificreports/ here, we account for multipole terms of the transfer process and analyze the effect of the donor-acceptor system geometry on transfer efficiency.

Model for confined excitons
We start our analysis from delocalized excitons in TMD monolayer, by focusing on excitons formed by states at the bottom of the conduction band and the topmost valence band at K ± points of the first Brillouin zone. Using the two-band effective mass model the wave function of the exciton can be written as a factorization of the relative motion of charge carriers and their center-of-mass motion 38,39 : where α = ±1 is the valley index, r e(h) is the radius vector of the electron (hole), σ is the normalization area, R = (m e r e + m h r h )/M is the center-of-mass vector, m e(h) is the effective mass of the electron (hole), M = m e + m h , Q is the wave vector of the center-of-mass motion, ψ N (ρ) is the wave function of the relative motion with coordinate ρ = r e − r h , and u α (r e ) and v α (r h ) are the Bloch functions of the electron and hole in valley α . For very small QDs (nanoflakes) the effective mass approach is not applicable, but it is possible to use ab-initio calculations to study their properties 40 . We solve the Schrödinger equation for the relative motion of states with circular symmetry, namely S-states with zero angular momentum, where ρ = |ρ| , µ = m e m h /M is the reduced mass, ǫ N is the eigenenergy of the S-state with the principal quantum number N. The nonlocally screened electron-hole interaction is described by the Rytova-Keldysh potential [41][42][43] where e is the electron charge, ρ 0 is the screening length, ε is the effective dielectric constant, and H 0 (x) and Y 0 (x) are Struve and Neumann functions.
Without considering the details of lateral confinement we focus on TMD QDs with in-plane localization of charge carriers and make the approximation of a harmonic confinement. If the coordinates of the relative motion, ρ , and the center-of-mass motion, R , of the confined exciton are not separated, one can use a variational procedure without the separation of coordinates to find the wave function 44 . Here we assume for simplicity that both the electron and hole are confined by parabolic potentials of the form which are characterized by the confinement frequency ω . With this potential we separate the coordinates of the relative motion and the center-of-mass motion 45 . Therefore, the energies and wave functions of the excitons in QDs are written as: respectively. The total energy of the exciton confined in the QD takes discrete values and is dependent on the band gap of the TMD monolayer, E g , and on E nl and ǫ N , which are the energies of the confined center-of-mass and relative motions.
We solve the Schrödinger equation for the center-of-mass motion in polar coordinates R ≡ (R, θ) . The exact eigenenergies and eigenstates are those of the 2D harmonic oscillator 46 , namely A nl = 2 · n!/(|l| + n)!, www.nature.com/scientificreports/ where n = 0, 1, 2, . . . and l = 0, ±1, ±2, . . . are the principal and angular momentum quantum numbers, respectively, L = √ /(Mω) is the QD size, and L |l| n (x) are the associated Laguerre polynomials. The radial part of the relative motion with zero angular momentum is determined not only by the nonlocally screened potential from Eq. (3) as in the case of free excitons, but also by the parabolic potential µω 2 ρ 2 /2 . The Schrödinger equation for the relative motion is written as To find an approximate solution of this equation one can use the 2D hydrogen-like wave functions with the Bohr radius as variational parameter 5 . Using the material parameters from Table 1, we solve this equation numerically for the first three exciton S-states in four specific TMDs, namely MoS 2 , MoSe 2 , WSe 2 , and MoTe 2 .
The energies ǫ N and wave function overlaps of the electron and hole at the same spatial position |ψ N (0)| 2 are shown in Fig. 1 as functions of the QD size L. In the top panel of Fig. 1 we observe for all materials the same trends for the energies ǫ N , they decrease with the size of the QD. The 1S state (N = 1) is less dependent on the QD size as it is the most localized state. The wave function overlaps |ψ 1 (0)| 2 are shown in Figs. 1e-h. Due to the confinement effect these overlaps are larger for smaller QDs and decrease monotonically with size. QDs in WSe 2 exhibit larger size effect than in MoS 2 , MoSe 2 , and MoTe 2 due to smaller reduced mass. For large QDs all three states converge to those of delocalized excitons shown by the dashed lines in Fig. 1.

Resonant energy transfer
In the following, we calculate the nonradiative resonant energy transfer between two 2D QDs coupled by a Coulomb potential. Due to the finiteness of QD dimensions in the xy-plane, the point dipole model developed by Förster leads to large errors when the distance between the QDs is of the order of their sizes. Therefore, the multipole nature of Coulomb interaction must be taken into account as in the case of conventional 3D QDs 36 . For conventional colloidal QD and molecular systems the orientations of QDs and molecules are random, whereas in layered systems the positions of the QDs are fixed and the localized excitons with lowest energies have an in-plane  www.nature.com/scientificreports/ circular polarization with sign reversal for K + and K − . This in-plane arrangement of the dipole moments leads to characteristic valley-selective orientation effects, which are considered in detail below. Let us consider two 2D QDs, namely donor and acceptor with sizes L D and L A , located in the same or in two different layers of the same TMD material. We assume the centers of QDs to be separated by a 2D-vector d in the xy-plane as shown in Fig. 2a. In the z-direction they are separated by a distance h as illustrated in the xz-plane projection. In our analysis, it is useful to distinguish two limiting cases. The first case corresponds to the situation when the QDs lie in different layers on top of each other and d = 0 . The second case is realized when the QDs are located in one monolayer and h = 0 . These two limiting cases allows us to analyze the orientation effects, e.g., for h = 0 the system is a truly 2D-object, in which the QDs and their dipole moments are in one plane. For d = 0 the problem is quasi-3D, because in principle all dimensions matter.
Without loss of generality we assume that only the 1S states (N = 1) contribute to the energy transfer between the donor and acceptor. Then, energy transfer is related to the annihilation of an exciton F (α) 1,nl in the donor and the creation of an exciton F (β) 1,mk in the acceptor (see Fig. 2b), where α and nl (β and mk) are the valley index and two quantum numbers of the exciton in the donor (acceptor). Here we consider the range of distances between the donor and acceptor, s = (d 2 + h 2 ) 1/2 , much smaller that the de Broglie wavelength of the annihilated exciton and neglect effects of exchange and radiation transfer 47,48 , but take into the account the multipole nature of the Coulomb interaction 36 . With these assumptions, the energy transfer depends on the matrix elements of the Coulomb potential where r 1 and r 2 are the 2D center-of-mass vectors originating at the centers of the donor and acceptor. By using the Fourier expansion we find the matrix elements where If we use the long wave approximation qa ≪ 1 with the lattice constant of TMD a and express the integration as a sum of integrals over elementary cells, we simplify F D (q) = dre iqr � nl (r)u α (r)v α (r). and is the area of the unit cell. The corresponding expression for F A (q) can be found by replacements D → A , α → β , and nl → mk in Eq. (17). By carrying out the integration over the angular variable of q in Eq. (15) we obtain with two integrals I 1 and I 2 given by where J η (x) is the Bessel function of the first kind, η = 1, 2 , and n is the unit vector co-directional with d.
Within the two-band approximation of the band structure in TMD monolayer, the interband matrix elements of the dipole moment operator in the donor and acceptor are written in Cartesian coordinates as where D = eat/E g and t is the nearest-neighbor hopping integral 49,50 , and angles ϕ and ϑ determine the crystal coordination axes of the donor and acceptor. Using these formulas for dipole moments and choosing n = e x , we find where δ αβ is the Kronecker symbol and p = e i(βϑ−αϕ) is the phase factor, which is determined by the alignment of the crystal coordination axes. Further, for the sake of simplicity, we assume ϕ = ϑ = 0 and p = 1 . According to the result, the intravalley matrix element (α = β) corresponds to the annihilation and creation of excitons in the same valley and is proportional to the difference I 1 − I 2 /2 . The intervalley matrix element (α = β) is proportional to −I 2 /2 , because the dipole moments of the excitons are orthogonal to each other and the first scalar product in Eq. (19) is zero. It is noteworthy that if we consider the interlayer excitons in TMD homo-or heterobilayers instead of intralayer ones, we must take into account the z-component of dipole moments in Eqs. (21) and (22).
Finally, using the Fermi's golden rule we obtain the rate of the resonant energy transfer from the donor state α, nl to all final states of acceptor β, mk as: where Ŵ is the sum of the total dephasing rates of interband transitions in the donor and acceptor, and = E nl − E mk is the energy detuning between the exciton levels in the donor and acceptor involved in the energy transfer process. Notably, if the magnitudes of the Coulomb matrix elements in Eq. (24) are much larger than and Ŵ , the formation of the entangled states in QDs must be considered 51 . Assuming the simplified case when the energy transfer occurs for equal QDs, L D = L A = L , from the state nl = 00 in valley α to the states with mk = 00 in valleys β = ±α , we find where the rate is dependent on the intravalley and the intervalley matrix elements. To quantify these matrix elements for h > 0 and d > 0 we evaluate the integrals in Eq. (20) numerically. It should be noted that all other S-states of the relative motion and all other mk states of the center-of-mass motion in acceptor contribute to the energy transfer, but their contributions are negligibly small.
Before proceeding with the evaluation of the integrals, it is instructive to introduce the dipole-dipole approximation (DDA) of the Coulomb interaction developed by Förster 34 . This approximation of point dipoles is obtained by increasing the distance between the donor and acceptor s to infinity, in particular d → ∞ for h = 0 or h → ∞ for d = 0 . Then we find the DDA limit for the matrix element from Eq. (23) and the energy transfer rate from Eq. (25) www.nature.com/scientificreports/ where the matrix element indexes were omitted for simplicity.
As illustrative examples of our theory, we consider two limiting cases shown in Fig. 2a, as they allow us to analyze the matrix elements and rates, and find their asymptotics. For h = 0 , when both QDs are located in the same monolayer, we find exact expressions where C 0 = 2πD 2 |ψ 1 (0)| 2 /(εL) , ξ = d 2 /(8L 2 ) , and I 0 (ξ ) and I 1 (ξ ) are the modified Bessel functions of the first kind. In Fig. 3 we show the matrix elements for MoS 2 monolayer and ϕ = 0 starting from the close-contact distances between QDs, d = 2L . For the shown range of distances the intravalley matrix element C 1 first starts with a positive value and decreases with distance d. After crossing the zero at d ≈ 2.51L it reaches a minimum C 1 ≈ −0.03C 0 at d ≈ 3.41L and further increases monotonically. On the other hand, the intervalley matrix element C 2 starts from a negative value, decreases monotonically with distance, and after reaching a minimum C 2 ≈ −0.15C 0 at d ≈ 2.26L it increases monotonically.
The analysis shows that for ξ → ∞ , which is the DDA in Eq. (26) for h = 0 , we obtain For small distances in the monolayer limit ( ξ ≪ 1 ) exchange effects must be taken into account, because they substantially change the intravalley matrix element. By substitution of Eqs. (28) and (29) into (25) we find the rate of the energy transfer for h = 0 as where γ 0 = 2C 2 0 /( 2 Ŵ) . By using the material parameters of monolayer MoS 2 and assuming L = 3 nm, ε = 4.5 and Ŵ = 5 meV 20 , we estimate the absolute values C 0 = 17 meV and γ 0 = 176 ps −1 . Again for ξ → ∞ we obtain the DDA in Eq. (27) γ (dd) d = 10γ 0 (L/d) 6 . This result clearly shows that for this limit, the system can be considered as two interacting point dipoles in the Förster model. We plot the rate of the energy transfer γ d in Fig. 4a and show explicitly the intravalley and intervalley contributions to the transfer. Evidently, the intervalley transfer is larger than the intravalley, particularly for large distances it is nearly one order of magnitude larger. It should be noted, that for the distance d ≈ 2.51L only the intervalley transfer occurs, this position is marked by the red arrow in Fig. 4a.
To compare the exact result for the energy transfer in Eq. (31) with the DDA we calculate the ratio between the energy transfer rates γ /γ (dd) d in Fig. 4b. The gray dashed line corresponds to the DDA asymptotic. For the monolayer limit with h = 0 we use Eq. (31), whereas for QDs separated in the z-direction with distances www.nature.com/scientificreports/ h = 0.1L, 0.5L, L we substitute Eqs. (23) and (20) into (25). Evidently, for small distances between the QDs the multipole contribution becomes sizable and the distance dependence deviates from the DDA. It is a result of the finite sizes of the QDs in the xy-plane. With increasing distance d the ratios γ /γ (dd) d become larger than unity, and also have maxima, e.g. for the monolayer limit we observe the maximum γ d ≈ 2.35γ (dd) d at d ≈ 4.14L . With increasing separation in the z-direction, h, we observe a decrease in the energy transfer and a shift of the ratio maximum to larger distances d. Summarizing the above, the monolayer limit exhibits the largest energy transfer rate, that is up to 2.35 times larger than the energy transfer for the DDA.
Another limiting case with d = 0 is shown in Fig. 2a and corresponds to the geometry, where the centers of the two QDs lie on the z-axis. For this symmetry only the intravalley energy transfer occurs, and its matrix element is given by where ζ = h/(2L) and Erfc(ζ ) is the complementary error function. Using Eq. (25) we calculate the transfer rate as From the asymptotic behavior of these functions for ζ → ∞ , which corresponds the DDA in Eqs. (26) and (27) for d = 0 , we obtain 2C 0 (L/h) 3 for the intravalley matrix element and γ (dd) h = 4γ 0 (L/h) 6 for the energy transfer rate. In Fig. 4c we show that with increasing distance h the transfer rate γ h decreases subexponentially. For all distances this ratio is smaller than unity and the DDA overestimates the energy transfer (see Fig. 4d). This holds generally if all multipole contributions are taken into account, which is necessary when the distance s is in the order of the QD sizes and the energy transfer deviates from the Förster model.

Summary and conclusions
In summary, we developed a simple theory for excitons confined in QDs of atomically thin TMD semiconductors. Using the approximation of harmonic confinement we derived the energy spectrum and the wave functions of QD excitons. By calculating the intravalley and intervalley Coulomb matrix elements and taking into account not only dipolar contributions but also multipole corrections, we determined resonant energy transfer rates between two adjacent QDs. We derived exact expressions for two simplified cases of possible donor-acceptor geometries, and found that for small distances all multipoles must be taken into account. At large distances, the energy transfer was found to converge asymptotically towards the Förster DDA. The largest energy transfer rate was found in the monolayer limit, where the major contribution stems from the intervalley matrix element of the Coulomb interaction and the intravalley matrix element is small. Moreover, the intravalley matrix element can be drastically suppressed by geometry, rendering the energy transfer valley selective. This aspect is important when considering collective light-matter effects such as superradiance for QD ensembles, and might prove useful for spin-valley selective transfer of quantum information.