Quenching of exciton recombination in strained two-dimensional monochalcogenides

TL;DR

This paper predicts that strain-induced localization in two-dimensional monochalcogenides can produce long-lived excitons with high binding energies, significantly reducing recombination rates.

Contribution

It introduces a theoretical framework showing how local strain fields in 2D monochalcogenides create long-lived, highly bound excitons by spatially separating electrons and holes.

Findings

Localized excitons have high binding energies.

Strain fields hinder exciton recombination.

Long-lived excitons are predicted in 2D monochalcogenides.

Abstract

We predict that long-lived excitons with very large binding energies can also exist in a single or few layers of monochalcogenides such as GaSe. Our theoretical study shows that excitons confined by a radial local strain field are unable to recombine despite of electrons and holes co-existing in space. The localized single-particle states are calculated in the envelope function approximation based on a three-band $\boldsymbol{k}\cdot \boldsymbol{p}$ Hamiltonian obtained from DFT calculations. The binding energy and the decay rate of the exciton ground state are computed after including correlations in the basis of electron-hole pairs. The interplay between the localized strain and the caldera-type valence band, characteristic of few-layered monochalcogenides, creates localized electron and hole states with very different quantum numbers which hinders the recombination even for singlet excitons.

Figures

Equations

$$H^{\text{eff}}(\mathbf{k})=\left(\begin{array}{ccc}A(\mathbf{k})& 0& B(\mathbf{k})\\ 0& C(\mathbf{k})& 0\\ B^{\ast }(\mathbf{k})& 0& D(\mathbf{k})\end{array}\right),$$

$$\begin{array}{rl}A(\mathbf{k})& =-1.4+19.6k_x^2+20.0k_y^2,\\ B(\mathbf{k})& =55.5k_x^2+63.4k_y^2,\\ C(\mathbf{k})& =1.4+67.7k_x^2+164.8k_y^2,\\ D(\mathbf{k})& =2.5+74.3k_x^2+148.8k_y^2\end{array}$$

$$\sum _m{H}_{nm}^{\text{SQD}}(\mathbf{r},-i{\nabla }_{\mathbf{r}})F_m^{(\gamma )}(\mathbf{r})={\epsilon }_{\gamma }{F}_n^{(\gamma )}(\mathbf{r}).$$

$$H_{nn}^{\text{SQD}}=E_n+{\Delta }_n(\mathbf{r})+\sum _{\alpha }{M}_n^{\alpha }{\nabla }_{\alpha }^2,$$

$$F_n^{(\gamma )}=\sum _{n_x{n}_y}{c}_{n_x{n}_y}^{(n\gamma )}{\varphi }_{n_x{n}_y},$$

$$H=H^{\text{SQD}}+H^{\text{e-h}},$$

$$\begin{array}{rl}{H}^{\text{SQD}}& =\sum _i{ϵ}_i{c}_i^{†}{c}_i+\sum _j{ϵ}_j{d}_j^{†}{d}_j,\\ H^{\text{e-h}}=& -\sum _{\begin{array}{c}{i}_1,j_1,i_2,j_2\\ \sigma ,{\sigma }^{\prime }\end{array}}{c}_{i_1\sigma }^{†}{c}_{i_2\sigma }{d}_{j_2(-{\sigma }^{\prime })}^{†}{d}_{j_1(-{\sigma }^{\prime })}{V}_{i_1,j_1,j_2,i_2}\\ & +\sum _{\begin{array}{c}{i}_1,j_1,i_2,j_2\\ \sigma ,{\sigma }^{\prime }\end{array}}{S}_{\sigma {\sigma }^{\prime }}{c}_{i_1\sigma }^{†}{c}_{i_2{\sigma }^{\prime }}{d}_{j_2(-\sigma )}^{†}{d}_{j_1(-{\sigma }^{\prime })}{V}_{i_1,j_1,i_2,j_2},\end{array}$$

$$|X⟩=\sum _{i,j,\sigma }{A}_{ij}{c}_{i\sigma }^{†}{d}_{j(-\sigma )}^{†}|0⟩,$$

$$I\propto\frac{1}{E_{X}^{(0)}}\sum_{\bm{q}\lambda}\Big{|}\sum_{i,j}A_{ij}p_{ji}^{(\bm{q}\lambda)}\Big{|}^{2}\delta(E_{X}^{(0)}-\hbar\omega_{q}),$$

$${\psi }_{\gamma }=\sum _{n\in A}{F}_n^{(\gamma )}{u}_{n0}+\sum _{m\notin A}{\tilde{F}}_m^{(\gamma )}{u}_{m0}.$$

$${\tilde{F}}_m^{(\gamma )}=-i\frac{\hslash }{m_0}\sum _l^{N_A}\frac{1}{E_l-E_m}(p_{ml}^x{\partial }_x+p_{ml}^y{\partial }_y)F_l^{(\gamma )}.$$

$$p_{ji}^{(\mathbf{q}\lambda )}=\sum _{m,n\in A}⟨u_m|{\mathbf{ϵ}}_{\mathbf{q}\lambda }\cdot \mathbf{p}|u_n⟩⟨F_m^{(j)}⟩F_n^{(i)}$$

$$p_{ji}^{(\mathbf{q}\lambda )}\approx ⟨u_{v_0}|p_z\cos {\theta }_z|u_{c_0}⟩⟨F_{v_0}^{(j)}⟩F_{c_0}^{(i)},$$

$$I\propto \frac{|{\sum }_{ij}⟨F_{v_0}^{(j)}⟩F_{c_0}^{(i)}{|}^2}{E_X^{(0)}}\delta (E_X^{(0)}-\hslash {\omega }_q)$$

Full text

††preprint: APS/123-QED

Introduction.— Long-lived excitons have been long sought both for creating quantum Bose gases in semiconducting materials Lozovik and Yudson (1976); Butov et al. (2001) and for the development of novel excitonic devices. To date, these have been extensively studied in III-V and II-VI semiconductor heterostructures where the range of existence is limited to low temperatures. Very recently, they have also been observed at room temperature in van der Waals transition-metal dichalcogenide (TMD) two-dimensional (2D) heterostructures Unuchek et al. (2018); Calman et al. (2018). In both cases the long-lived nature stems from the fact that the electron and the hole are spatially separated, thereof the name of indirect excitons.

The spatial separation between electron and hole composing an indirect or interlayer exciton in TMD 2D heterostructures allows for long recombination life-times, but it is also a limitation. Leaving aside the technical difficulties of creating heterostructures in a controlled manner, the required spatial separation limits, in turn, the binding energy. Larger inter-layer or electron-hole distances amount to longer recombination times, but weaker binding energies. Here we show that both requirements, impossible to meet in TMD heterostructures can be found in strained few-layered monochalcogenides such as GaSe.

It is well known that the application of external strain on a semiconductor leads to a modification of the band gap. This is particularly important in 2D crystals where large strains can be applied before the crystal rupture Lee et al. (2008); Bertolazzi et al. (2011). The modulation of the band gap due to strain has been quantified for MoS2 Castellanos-Gomez et al. (2013), phosphorene Quereda et al. (2016), and GaSe Maeso et al. (2019), to name a few. In the multilayer form of the latter and in other 2D TMDs, a causality relation between localized strain deformations and single-photon emission has been demonstrated Kern et al. (2016); Tonndorf et al. (2017a, 2015); Kumar et al. (2015); Tonndorf et al. (2017b), although the ultimate reason for this relation remains unclear.

The proven strong band-gap modulation with strain, the possibility of single-photon emission, along with an anomalous shape of the GaSe valence band, motivates the present study. Below approximately 7 layers, the valence band turns from a common inverted parabola to a caldera-type shape or ring-type shape Rybkovskiy et al. (2014); Maeso et al. (2019). Free excitons in GaSe monolayer have been already studied with ab-initio methods Antonius et al. (2018), showing a variety of behaviours due to the peculiar valence band. We are interested here in localized radial deformations as those invoked to explain single-photon emission in multilayer GaSe Kern et al. (2016); Tonndorf et al. (2017a). Here we show that the combined peculiarity of the valence band along with local strain quenches the photoluminescence by making the localized excitons extremely long-lived, as shown below.

The effective Hamiltonian.— We start by constructing an accurate low-energy Hamiltonian for a monolayer of GaSe using $\bm{k}\cdot\bm{p}$ theory. Firstly, we find the DFT band structure of monolayer GaSe as shown in Fig. 1 (blue solid lines). We have used the CRYSTAL code implementation of the hybrid functional HSE06 Dovesi et al. (2014) to obtain a good approximation to the actual quasi-particle gap. Our DFT calculation will serve a two-fold purpose: to check the validity of our effective model Hamiltonian and to obtain the necessary information (Bloch states and band edges at the $\Gamma$-point) to compute the momentum matrix elements for a parameter-free description. The unperturbed Hamiltonian is $H^{(0)}_{nm}=E_{n}\delta_{nm}$, where the $n,m$ indices represent Bloch states at the $\Gamma$-point, which we will use as our basis, and $E_{n}$ is the energy of these Bloch states. Following the discussion in Ref. Li and Appelbaum (2015), we pursue a 3-band model to properly account for the caldera form of the valence band (see Fig. 1). In $\bm{k}\cdot\bm{p}$ theory, the first order correction is linearly dependent in $k$ with the form $H^{(1)}_{nm}=\frac{\hbar}{m_{0}}\sum_{\alpha}k_{\alpha}.{p}^{\alpha}_{nm}$, where $\alpha$ (and later $\beta$) denotes Cartesian coordinates and $p^{\alpha}_{nm}$ is the $\alpha$ component of the momentum matrix element between $\Gamma$-point Bloch states.

Based on group theory considerations we can anticipate the momentum matrix elements that are non-zero and, thereby, the general form of the Hamiltonian. We will follow the notation in Ref. Li and Appelbaum (2015) for the $D_{3h}$ point group, to which monolayer GaSe belongs. The valence and the two first conduction bands, $v_{0}$, $c_{0}$ and $c_{1}$, belong to the $\Gamma^{+}_{1}$, $\Gamma^{-}_{2}$, and $\Gamma^{+}_{1}$ irreducible representations (IRs) respectively, and will form our working subspace $A$ from now on. In order for the momentum matrix elements to be non-zero, the decomposition of the direct product of the eigenstates IRs and the IR of the momentum operators, which in the $D_{3h}$ group transform as $\Gamma^{+}_{3}$ ($p_{x}$, $p_{y}$) and $\Gamma^{-}_{2}$ ($p_{z}$), must contain the fully symmetrical representation $\Gamma^{+}_{1}$. One can see that only $\braket{\Gamma_{3}^{+}}{p_{x,y}}{\Gamma^{+}_{1}}$, $\braket{\Gamma^{+}_{3}}{p_{x,y}}{\Gamma^{+}_{2}}$, $\braket{\Gamma_{3}^{-}}{p_{x,y}}{\Gamma^{-}_{1}}$, $\braket{\Gamma^{-}_{3}}{p_{x,y}}{\Gamma^{-}_{2}}$, $\braket{\Gamma_{3}^{+}}{p_{x,y}}{\Gamma_{3}^{+}}$, $\braket{\Gamma_{3}^{-}}{p_{x,y}}{\Gamma_{3}^{-}}$ and their complex conjugates are non-zero. Therefore, the resulting Hamiltonian will not contain terms linear in $k$. The absence of linear terms in the Hamiltonian forces us to consider the coupling of the subspace $A$ with the rest of the bands, which now becomes the dominant contribution. This is typically done through the Löwdin partitioning technique Glutsch (2004). Including the coupling with all bands and evaluating numerically the momentum matrix elements, the effective Hamiltonian becomes

[TABLE]

with

[TABLE]

in units of eV (and atomic units for wave vectors). As discussed earlier, $v_{0}$ and $c_{1}$ bands transform as $\Gamma_{1}^{+}$ which only couples with $\Gamma_{3}^{+}$ states. On the other hand, $\Gamma_{3}^{+}$ does not couple with $\Gamma^{-}_{2}$, the IR of $c_{0}$. As a result, one could have expected $v_{0}$ and $c_{1}$ to remain decoupled from $c_{0}$, as can actually be seen in Eq. (1). The three bands resulting from the effective Hamiltonian are shown in Fig. 1 (red dashed lines). The agreement with the original DFT bands is very good near the $\Gamma$-point, included the caldera shape of the valence band. Since we have used all the bands in the evaluation of the second order perturbation terms, no further improvement can be envision to second order within our 3-band model.

Local strain model.— Assuming that strain changes occur on length scales much larger than the lattice spacing we can work in the envelope function framework where one writes the eigenfuntions in the form $\psi_{\gamma}=\sum_{n}F^{(\gamma)}_{n}(\bm{r})u_{n\bm{0}}(\bm{r})$. Here $u_{n\bm{0}}$ are the periodic Bloch functions at $\bm{k}=0$ and $F_{n}^{(\gamma)}$ are 2D envelope functions 111Note that in the expansion of the SQD states, the envelope functions are defined in 2D while the Bloch periodic functions are ab-initio 3D functions. calculated from a set of coupled differential equations:

[TABLE]

$H^{\text{SQD}}$ is the hamiltonian of the strain-induced quantum dot (SQD) created by the local perturbation. We assume the strain to be locally biaxial which can be modeled by adding a local “potential” to the $k$-independent diagonal terms in $H^{\text{eff}}$,

[TABLE]

while keeping unchanged the $k$-dependent diagonal (here generically represented by $M$) and non-diagonal terms $B^{(*)}(\bm{k})$. This approximation captures the main effect, which is the confinement of electrons and holes, while maintaining the complexity of the implementation to a minimum. It neglects changes in the effective mass or band curvature and the possible appearance of new Hamiltonian terms due to the breaking of the $C_{3h}$ symmetry, but it has been shown to be quantitatively accurate in a somewhat similar context Tadić et al. (2002). In practice we consider a tensile gaussian strain field of the form $S_{0}\exp{(-r^{2}/\sigma_{s}^{2})}$ with a maximum strain at its center of 5%. The local lattice deformation is mapped onto $\Delta_{n}(r)$ thanks to a series of DFT calculations performed with biaxial tensile strain. In the inset of Fig. 1 we show the $\Gamma$-point energies as obtained for several values of the tensile strain. We have also implicitly assumed, as seen in the inset of Fig. 1, that the chemical potential (set at zero) stays in the mid-gap at all points in space. A detailed calculation is beyond the scope of this work, but this choice complies with local charge neutrality and overall thermodynamic equilibrium at low temperatures.

Single-particle results.—

In order to solve Eq. (3) we further expand the envelope functions as

[TABLE]

where $\phi_{n_{x}n_{y}}$ are the cartesian 2D harmonic oscillator (HO) basis functions. The kinetic part of the matrix elements in the HO basis can be easily evaluated using the ladder operators, while the strain contribution coming from the term $\Delta_{n}(r)$ is evaluated numerically. As an example we show in Fig. 2 the density of states of electrons and holes, obtained from the numerical diagonalization of Eq. (3) for two representative cases. In both cases we find bound states within the band gap of the unstrained monolayer. These states distribute in a more dense spectrum in the valence band than in the conduction band, a consequence of the caldera shape of the former, which also gives rise to the expected van Hove singularity Antonius et al. (2018). The number of bound states in each band increases with $\sigma_{s}$ and $S_{0}$, as could intuitively be expected. The transition from bound to delocalized states (shaded regions in Fig. 2) occurs near the band edges of the unstrained conduction and valence bands. A minimum of $\sim 30$ and $\sim 100$ bound states have been found in the valence and conduction band respectively in all the studied cases, which is sufficient for our purposes.

A plot of the envelope functions reveals the very different nature of valence and conduction bound states. In Fig. 3(a) and 3(c) we show a chart of the HO basis coefficients (absolute value) of the $v_{0}$ and $c_{1}$ envelope functions corresponding to the highest-energy hole state and in Fig. 3(e) the $c_{1}$ envelope function of the lowest-energy electron state. Both hole envelope functions show a very similar behaviour, although $v_{0}$ gives the main contribution to the state. The other envelope function is zero due to the absence of coupling in Eq. (1). The electron state has mainly $(n_{x}=0,n_{y}=0)$ character (as could be anticipated for a parabolic band), but the hole state shows a dramatic shift to high $n_{y}$ quantum numbers. This behaviour can be traced back to the caldera ring away from the $\Gamma$-point in the unstrained band structure. We also show the 20th electron and hole states (moving away from the band edges) in Fig. 3(b, d and f). The hole state spreads now also to higher values of $n_{x}$ while the electron state still presents the largest weight at low HO quantum numbers.

Excitons.— Interactions are now added to the single-particle Hamiltonian expanded by electron and hole bound states:

[TABLE]

with

[TABLE]

where $c_{i}$ ($d_{j}$) operators act on electrons (holes) and the phase factor $S_{\sigma\sigma^{\prime}}=(-1)^{(1/2-\sigma-\sigma^{\prime})}$ preserves essential symmetries in the electron-hole language Fetter and Walecka (2012). Note that the renormalization of the single-particle energies into quasi-particle excitations Rohlfing and Louie (2000) is already approximately accounted for in our hybrid DFT band structure. No other interaction terms are needed since we only consider single electron-hole pairs. At this level, excitonic states can be constructed as

[TABLE]

where $\ket{0}$ denotes the vacuum state with no electrons and holes. The excitonic states are then calculated by projecting Eq. (7) onto the electron-hole pair basis and converged in the number of these pairs (we restrict ourselves to 200 electron-hole pairs which is enough for our purposes). Note that, to leading order, the exchange part in Eq. (7) will not contribute to the matrix elements Azizi and Machnikowski (2015); Sheng et al. (2005). This is due to the specific form of our Hamiltonian in Eq. (1), which results in conduction and valence states not sharing envelope functions. As a consequence, singlet and triplet states of the exciton are degenerate in energy. States with well defined total spin can be constructed a posteriori, if desired.

For the interaction kernel, we have considered the bare Coulomb interaction and the Keldysh interaction model, more appropriate for 2D systems Keldysh (1979). Both types of interactions give qualitatively similar results with unimportant differences in the actual exciton binding energies which are not relevant for our discussion here. Since the interactions preserve all the symmetries, excitons can be classified according to how the transform under in-plane reflections in both electron and hole coordinates. In the following we focus in the excitons whose electron and hole wavefunctions transform similarly, as these are, in principle, the active ones in photoluminescence. In Fig. 4(a) we show the quasiparticle gap and binding energy of the optically active (symmetry-allowed) lowest-energy excitonic state $E_{X}^{(0)}$ for different parameters characterizing our SQD. Depending on the parameters and interaction details this can be the ground state, but this is unimportant for our main conclusion (see below). In particular we have used the Keldysh potentialPrada et al. (2015) in an environment of vacuum. The binding energies decrease with $\sigma_{s}$ since the single-particle states become more extended in space. On the other hand, the binding energy increases with the strength of the strain as they become more localized.

Quenching of the exciton recombination.— The most relevant result is found in the photoluminescence. The recombination rate to the ground state mediated by the emission of a photon can be computed using Fermi’s golden rule (ignoring proportionality constants) Hawrylak (1999); Sheng et al. (2005); Spataru et al. (2005)

[TABLE]

Here $p_{ji}^{(\bm{q}\lambda)}$ is the scalar product of the momentum matrix element and the photon polarization vector $\bm{\epsilon}_{\bm{q}\lambda}$. The squared quantity is known as the oscillator strength. Thus, the recombination is encoded in those electron-hole pairs which have a relevant contribution to the exciton wave function. For completeness and to strengthen our conclusion, we invoke now the first order correction to the bound states Glutsch (2004),

[TABLE]

The first term represents the main contribution to the bound states, as explained previously, while the second one now includes Bloch functions out of the subspace $A$ of the effective Hamiltonian. The new envelope functions can be calculated as a linear combination of the derivatives of those in $A$,

[TABLE]

With this the momentum matrix element in Eq. (9) now becomes (up to first order in $\bm{\nabla_{r}}$)

[TABLE]

The intra-band terms (not shown), which are usually neglected Sheng et al. (2005), here they are strictly zero because of the structure of our Hamiltonian. For in-plane polarization, we have already seen that the momentum matrix element between states in $A$ vanish. If the polarization vector has a $z$-component, then the only coupling allowed is between $v_{0}$ and $c_{0}$ Antonius et al. (2018). Thus in the leading order

[TABLE]

being $\cos\theta_{z}$ the $z$-direction cosine of the polarization vector. For a single polarization Eq. (9) reduces to

[TABLE]

Since we are considering symmetry-allowed active excitons, the scalar product can be non-zero. However, in our case the overlap is always vanishing small. Note again that the exciton is composed of electron-hole pairs with respective wavefunctions of very different HO quantum numbers, at least in the relevant range of energies. Therefore, at this point we can conjecture that even the dipole-allowed transitions will vanish due to the nearly orthogonal nature of the hole and electron envelope functions. In Fig. 4(b) we show the rate of Eq. (14) in logarithmic scale. We obtain, as expected, very low values for all the cases studied. Note that in a common semiconductor with parabolic bands, one would expect the overlap between electron and hole envelope functions to be of the order of one and therefore, for similar momentum matrix elements, a decay rate many orders of magnitude larger. Higher order corrections to the envelope functions are not expected to give a relevant contribution to Eq. (12). They are obtained by differentiation of the envelope functions, which translates in a spreading of the order of unity of the quantum numbers (see, e.g., Fig. 3), thus maintaining its quasi-orthogonal character. Finally note that this result does not depend on the singlet or triplet nature of the exciton or the magnitude of the dipolar matrix elements because it is encoded in the envelope part of the wave functions. We thus conclude that photoemission is generically quenched.

Conclusions.— We have shown that excitons localized by a radial strain perturbation in GaSe monolayer present extremely long recombination times. This is a direct consequence of the very different nature of the wave function of electrons and holes. While the parabolic conduction band presents bound states with low HO quantum numbers, the caldera-shaped valence band makes the hole states to have a more disperse structure of HO quantum numbers. This translates in an orthogonality of electron and hole envelope functions which, in the end, quenches the recombination. Although we have focused on a specific 2D monochalcogenide, everything seems to indicate that similar results would be obtained in other monochalcogenides such as GaS and InSe Rybkovskiy et al. (2014) and even in different families of materials presenting a caldera-shaped valence or conduction band.