Diversity of trion states and substrate effects in the optical properties of an MoS2 monolayer

Almost all experiments and future applications of transition metal dichalcogenide monolayers rely on a substrate for mechanical stability, which can significantly modify the optical spectra of the monolayer. Doping from the substrate might lead to the domination of the spectra by trions. Here we show by ab initio many-body theory that the negative trion (A−) splits into three excitations, with both inter- and intra-valley character, while the positive counterpart (A+) consists of only one inter-valley excitation. Furthermore, the substrate enhances the screening, which renormalizes both band gap and exciton as well as the trion-binding energies. We verify that these two effects do not perfectly cancel each other, but lead to red-shifts of the excitation energies for three different substrates ranging from a wide-bandgap semiconductor up to a metal. Our results explain recently found experimental splittings of the lowest trion line as well as excitation red-shifts on substrates.

Left side: Isolated MoS 2 monolayer with an excited electron (red) -hole (blue) pair. Also given is a sketch of the energy positions involved (not to scale). Right side: When the supercell size is reduced, the CBM is shifted to lower energies, but the binding energy is also reduced, leaving the excitation energy (A) unchanged.

Supplementary Note 1: The LDA+GdW approximation
All results presented in the main text are obtained within the framework of many-body perturbation theory (MBPT), which allows a reliable, ab initio description of excited electronic states [1,2]. We calculate the quasiparticle band structure starting from a density functional theory (DFT) calculation within the local density approximation (LDA) followed by the GW approach within the LDA+GdW approximation. The latter is numerically efficient and allows the inclusion of substrates. To handle trions, we use an extension of MBPT for excitons which includes a third excited particle. Both approaches are elaborated further followed by a careful analysis of the numerical stability of our results.
The LDA+GdW method is a first principles theoretical approach to the many-body physics of a given system. It states a further approximation in calculating the GW self-energy operator. Ref. [3] introduces it in detail showing its biggest advantage: numerical efficiency.
Furthermore, it allows the straightforward inclusion of screening effects of substrates.
Starting point for most GW calculations are density functional theory energies and wave functions. They are used as input to calculate the self-energy Σ = iGW , which replaces the DFT-exchange correlation potential V xc in the quasiparticle Hamiltonian The observation that motivates the LDA+GdW approach is that the usage of a hypothetical metallic screening approximately reproduces the DFT-LDA exchange correlation energy, i.e., iGW metal ≈ V xc (provided that iGW is a good approximation to Σ). This way, the construction ofĤ yields a good approximation for the quasiparticle Hamiltonian. Using Supplementary Eq. (2) to calculate the quasiparticle corrections states a huge advantage in efficiency, since now ∆Σ := iGdW (of the order of 1 eV) has to be calculated and converged instead of calculating V xc and iGW separately (both often of the order of 10 eV) and then subtracting them. Within LDA+GdW we employ a second approximation: the usage of a model dielectric function (ε) for W and W metal , which replaces the numerically demanding calculation of ε within, e.g., random-phase approximation. Due to the comparably small size of iGdW (compared to iGW ) the resulting error is tolerable, while such approximation would be more troublesome in full GW .
In the employed model, the dielectric function of the system is constructed by the contributions of the individual atoms. To simulate the screening of a substrate we only need to add the contributions of the substrate atoms to ε. Treated this way, the substrate only enters the calculation through the dielectric function, making it possible to get a clear picture of the changes induced by its screening only. We note that covalent bonding between the substrate and the monolayer is neglected. In principle, quantum-mechanical interaction of the electronic orbitals across the interface could add exciplex (or charge-transfer) contributions to the MoS 2 excitons, leading to red or blue shifts depending on orbital overlap (see Ref. [4] for similar effects in interacting carbon nanotubes). In practice, however, such effects would require perfect matching of the complex phase of the electronic states on both sides of the interface [4], which is already destroyed by the lattice mismatch between substrate and monolayer and can thus be considered as irrelevant.
To conduct a calculation, the unit cell of the substrate is stretched accordingly to fit the unit cell of the monolayer. The atomic polarizabilities were modified such that the substrate still has its natural dielectric constant. It was carefully checked that this distortion has a neglectable impact on the induced shifts due to screening. As explained by Rohlfing et al.
[4], the shifts due to additional screening are not expected to depend strongly on details of the lattice structure. As mentioned in the main text, we use values for the dielectric constant ∞ of 3.9 (SiO 2 [5]), 4.95 / 4.10 (ε ∞ /ε ⊥ ∞ h-BN [6]), and ∞ (Au). We also note in passing that our current LDA+GdW approach, although being an additional approximation to the GW approach on an absolute energy scale, is ideally suited to address the red-shifts (i.e., differences) with sufficient accuracy [4].

Supplementary Note 2: Electron-hole excitations
Excitons are bound states of an excited electron and a hole. A well established framework for their theoretical description is given by the GW /BSE approach [7,8]. Here, excitations are described through the Hamilton-matrix elements Where |vc :=â † câ v |0 denotes electron-hole pair excitations from the ground state of an N -electron system |0 . The index v (c) refers to occupied valence (empty conduction) states. For a periodic system v = (v, k v ) and c = (c, k c ) combine the band index (v, c) with a wave number (k v , k c ) for holes and electrons. The short-hand notation v (c) is chosen for brevity. The Hamilton matrix has two contributions: the first is the bandstructure termĤ BS containing the quasiparticle energies c , v . The second is the electron-hole interactionĤ eh which mixes the excitations and includes the bare Coulomb interaction Note that the discussion is restricted to the Tamm-Dancoff approximation.

Supplementary Note 3: Trion states
Trion states (also known as charged excitons) consist of three excited particles. Negatively charged states are formed from two electrons and one hole, while positively charged states are formed from one electron and two holes. In this section we focus on negatively charged states; the corresponding expressions for positive trions are analogous. A negative trion can be constructed as To avoid double counting, we require that c 1 < c 2 , which assumes that the empty states are ordered. The definition of the order is arbitrary but has to be kept. Setting up similar matrix elements as in Supplementary Eq. (3) yields A detailed derivation of Supplementary Eq. (6) can be found in Ref. [9] and its Supplemen- where the summation is restricted to c 1 < c 2 . The total momentum of the trion is given by K, which is constructed by the momenta of the electrons and the hole as which states a restriction to the sum in Supplementary Eq. (7). Note that the corresponding eigenvalues of Supplementary Eq. (7) denoted as E (T,K) are the energies of the trions. Optical transition energies found in an experimental spectrum would correspond to transitions from the trion state |T, K to a single excited electron state having the same momentum as the trion |cK . Transition energies are thus given by

Supplementary Note 4: Numerical results -band-structure calculations
Starting point for our many-body calculations is a DFT calculation carried out in the local density approximation (LDA), in the parametrization of Perdew and Zunger [10]. Normconserving pseudopotentials [11] in Kleinman-Bylander form [12] are used that also include spin-orbit interaction. We employ a basis of three shells of Gaussian orbitals with s, p, d . We find a theoretical lattice constant of 3.16Å that coincides with the experimental value of bulk MoS 2 3.16Å [13]. The position of the S atoms are structurally relaxed until forces are smaller than 10 −4 Ry a B . In reciprocal space a k-mesh of 10 × 10 × 1 is employed.
All calculations are carried out using a code written by ourselves [14]. The resulting bandstructure is depicted in dashed red in Supplementary Fig. 1 (a) showing a direct band gap at K of 1.82 eV. The DFT calculation provides the single particle wave functions and energies that are used as input for the subsequent quasiparticle calculations.
As mentioned above, our MBPT calculations are carried out within the LDA+GdW approximation. The resulting band structure in direct comparison with the DFT bands is shown in Supplementary Fig. 1 (a). In MBPT the appearance of two-point functions requires a second, auxiliary basis. We use plane waves with a cutoff of 2.5 Ry ≈ 205 plane waves, which shows a convergence of the bands better than 0.05 eV as can be seen in Supplementary Fig. 1 (b). To account for interlayer interaction between different super cells, that lead to a closing of the electronic band gap, we calculate the electronic structure for multiple interlayer distances L (up to 80Å), and extrapolate as L → ∞ (see Supplementary Fig. 1 (c)). We note that this extrapolation scheme is only needed for the electronic properties, i.e., the gap and band structure, the absorption spectra are almost independent on the interlayer distance (see below). Supplementary Fig. 1 (c) also shows the k-mesh convergence of the quasiparticle gap, which is converged better than 0.05 eV for meshs greater than 14 × 14 × 1.

Supplementary Note 5: Numerical results -BSE calculations
The excitonic properties are derived by solving Supplementary Eq. (7). Linearly polarized light (within the plane spanned by the monolayer) is used for all spectra. Supplementary   Fig. 2 (a) shows the calculated BSE spectrum for different supercell sizes and thus different interlayer distances L. The spectrum shows the prominent A, B, and C excitons. Note that a k-mesh of 15 × 15 × 1 is used to analyze the dependence on L. Most strikingly, the spectrum seems not to depend much on the interlayer distance (if L is large enough, i.e., 27Å). Supplementary Fig. 2 (c)  orange dots), leaving the excitation energies almost constant (black dots, labeled A and B). We note that the same holds for the trion states as well.
For a comparison with trion spectra two valence and four conduction bands are used.
Thereby the lower energy part of the spectrum, i.e., up to 2.6 eV is converged. Supplementary   Fig. 3 shows the A exciton energy for different k-grids on different substrates. Note that the A excitons are already converged better than 0.05 eV at meshs as small as 15 × 15 × 1, which we find equivalently for the B exciton. The excited A 2s state converges a bit slower, due to its more complicated spatial structure. The convergence of all three states is shown in Supplementary Fig. 4. This extremely fast convergence results from a combination of several aspects to our approach: In conventional GW calculations for d-electron systems, the s and p semicore electrons of the same shell (here: Mo 4s+p) have to be included as valence electrons to incorporate their bare exchange interaction. Within LDA+GdW , on the other hand, the bare exchange interaction is not considered explicitly, and the semicore states can be treated as core states. This reduces the number of valence electrons (per Mo atom) from 14 to only 6. Moreover, we use identical k-meshs for the quasiparticle corrections and the electronhole interaction, omitting the need for an interpolation scheme. The k-meshs for internal summation in the self-energy operator are also chosen exactly equivalent to the grids for the quasiparticle corrections. With such exactly matching meshs, the GW quasiparticle gap and the electron-hole interaction show the same asymptotic behavior with increasing mesh density, and the exciton energy (i.e., the combination of both) converges rapidly with the number of k-points. Additionally, the Coulomb interaction between adjacent supercells is not truncated, but exrapolated for an infinite interlayer distance for the electronic properties, which also reduces the convergence requirements [15].
Enhanced screening due to the substrate also affects the wave functions of the excited states. Due to stronger screening, the excitonic binding is weakened and the real space wave function gets further extended. This is illustrated in Supplementary Fig. 7 for the A exciton. The distribution of the electron in the plane spanned by the Mo atoms is shown in Supplementary Fig. 7 (a), assuming the hole is fixed close to the Mo atom in the center. In reciprocal space this trend is inverted. Supplementary Fig. 7 (b) show how the contribution in k-space of the A exciton gets more concentrated around K. The full width at half maximum for the vacuum situation (0.116 2π a ) is reduced by 4.1% or 8.7% for the SiO 2 or Au(111) substrate, respectively.

Supplementary Note 6: Numerical results -Trion calculations
In the last step, one additional electron is included and bound negative trion states are calculated (the treatment of positive trions is fully analogous). This requires a diagonalization of the trion Hamilton-operator described by Supplementary Eq. (6). For our presented trion spectra the dimension is ∼ N 2 k N v N 2 c ∼ 10 7 , using a grid of N k = 27 × 27 × 1 k-points and two valence and four conduction bands (N v , N c ), while the corresponding dimension of the BSE-Hamilton matrix would be merely N k N v N c ≈ 6000. The size of the Hamilton operator puts a direct diagonalization out of reach. Fortunately, about 99.8 % of the matrix elements are zero, enabling an efficient route through recursive schemes like the Haydock recursion [16] when focusing on spectra, or the Lanczos algorithm or similar [17] to access eigenstates. It was carefully checked that both approaches results in almost identical spectra (not distinguishable by eye on the presented energy scales). Sufficient processors and memory were supplied by the HPC system PALMA of the University of Münster and the super-computer JURECA at Jülich Supercomputing Centre (JSC).
Independently from the BSE, the convergence of the trion spectrum with respect to the used k-grid has to be carefully evaluated. The (Haydock-)spectra are shown in Supplementary Fig. 5. Not all trions converge at a similar rate, especially the first resonance shows a fast convergence, which becomes especially evident in Supplementary Fig. 6, which shows the energetic position of the A exciton and A − (1) trion in comparison, both extrapolated to N → ∞, which yields a trion energy of 2.07 eV and thus a trion binding energy of 58 meV.