Single-particle properties of topological Wannier excitons in bismuth chalcogenide nanosheets

We analyze the topology, dispersion, and optical selection rules of bulk Wannier excitons in nanosheets of Bi2Se3, a topological insulator in the family of the bismuth chalcogenides. Our main finding is that excitons also inherit the topology of the electronic bands, quantified by the skyrmion winding numbers of the constituent electron and hole pseudospins as a function of the total exciton momentum. The excitonic bands are found to be strongly indirect due to the band inversion of the underlying single-particle model. At zero total momentum, we predict that the s-wave and d-wave states of two exciton families are selectively bright under left- or right-circularly polarized light. We furthermore show that every s-wave exciton state consists of a quartet with a degenerate and quadratically dispersing nonchiral doublet, and a chiral doublet with one linearly dispersing mode as in transition metal dichalcogenides. Finally, we discuss the potential existence of topological edge states of chiral excitons arising from the bulk-boundary correspondence.


Introduction
Three-dimensional topological insulators, and all other topological materials for that matter, are presently receiving much attention because of their excellent prospects for energy-efficient electronics, (pseudo)spintronics devices, and quantum information processing [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].Prototypical examples of three-dimensional topological insulators are the bismuth chalcogenides Bi 2 Se 3 and Bi 2 Te 3 .Since in linear response these materials are ideally conducting only due to the presence of massless Dirac fermions on their surface, most experiments with topological insulators have focused on these unusual topologically protected surface states.However, the situation changes dramatically upon photoexcitation, as excitons and unbound charges may be produced in the bulk with a topologically nontrivial band structure.Consequently, it is important for the understanding of light-matter interactions to investigate also the bulk properties of topological insulators and in particular the precise topological nature of the excitons, whose presence or absence is crucial in optoelectronic devices such as lasers [19][20][21][22], light-emitting diodes [23][24][25], and photovoltaic cells [26][27][28][29].In the context of quantum information processing a particularly interesting question is if the exciton topology is transferred to the quantum state of the photons emitted via photo-or electro-luminescence [30].Apart from such applications, the many-body physics of topological excitons is thought to be very exciting and accessible experimentally by well-established pump-probe techniques.Important examples of interesting many-body states are the topological excitonic insulator [31][32][33][34][35] and a Bose-Einstein condensate of topological excitons or possibly even biexcitons [36].In recent years, therefore, there has been an increasing interest in the study of excitons formed in topological insulators.Semiclassically and within the effective-mass approximation, a pioneering approach has been to introduce Berry-phase corrections to the electron-hole interaction due to the topological band structure and determine the exciton energy spectrum [37][38][39][40][41][42].The approach presented in these works mostly focuses on the case of a total exciton momentum Q = 0, which is sufficient for studying optically active excitons, but not enough to obtain their global topological properties.The latter has been achieved for Frenkel excitons within a Hubbard-like model with on-site Coulomb interactions [43].However, because of the insulating nature of the bulk of bismuth chalcogenide nanosheets and their geometry, which leads to a lower-than-bulk dielectric constant due to the surrounding medium, we expect in these materials long-ranged electron-hole interactions allowing for the formation of so-called topological Wannier excitons [44,45].In this article we therefore study the concrete example of (quasi-)two-dimensional bulk excitons in Bi 2 Se 3 nanosheets.

Band structure of Bi 2 Se 3 nanosheets
Because we are only interested in the physics around the Γ point, we restrict ourselves to the bands closest to the Fermi surface.For all numerical purposes we choose a nanosheet thickness of 6 nm, which approximately corresponds to 6 quintuple layers (QLs).This is sufficiently large to preserve the nontrivial bulk topology, but also sufficiently small for the problem to be mainly regarded as two-dimensional.The effective Hamiltonian governing the physics around the Γ point in Bi 2 Se 3 nanosheets is equivalent to the BHZ Hamiltonian for the quantum spin Hall effect [46].As a consequence of spin-orbit coupling it features a band inversion around the Γ point as illustrated in Fig. 1.As a combined effect of the time-reversal and inversion symmetries of the system, the conduction and valence bands are each two-fold degenerate.We label each band by an index denoted as the spin-orbit parity and obtain four eigenstates |χ µ k , where µ ∈ {c, v} ∪ {+, −}. 2

Topological exciton eigenstates
The four different spin-orbit-parity combinations for an electron and a hole give four different exciton families |Q; s, t that generalize the singlet and triplet states in normal semiconductors.Here, s and t are the spinorbit parities of the bands where the electron and the hole are located, respectively, and Q ≡ (Q x , Q y ) is the total in-plane exciton momentum.In the absence of an exchange interaction, these states diagonalize the two-body Hamiltonian and have a well-defined Chern number equal to C s,t = s + t.Hence, in this idealized case there are two topologically nontrivial exciton states, characterized by a nontrivial winding of the momentum-space pseudospin texture Γ e,h (Q) (defined precisely in the Methods section).As illustrated in Fig. 2, the pseudospin texture gives an intuitive picture of the nontrivial topology of this basis of exciton states, as it allows us to visualize the total Chern number of each basis element as a combination of two winding numbers by looking at the path of the electron (hole) pseudospin from pointing upwards (downwards) at the origin Q = 0 to pointing downwards (upwards) at Q → ∞, represented as the circles in Fig. 2. Thus, the pseudospin texture of each constituent particle can be seen as a skyrmion, so that the topological exciton is represented by a double skyrmion texture.

Hole Electron
Figure 2: Momentum-space skyrmion textures of an electron and a hole.This idealized illustration corresponds to the configuration in a weakly bound exciton.Due to rotational symmetry the two circles actually represent only a radial slice of the two spheres onto which the complete momentum plane is mapped with a unit winding number.The south (north) pole of the top (bottom) sphere corresponds to the origin Q = 0, whereas the opposite pole corresponds to Q → ∞.Note that in a normal semiconductor like CdSe the pseudospins always point in the same direction, independent of Q, and thus have a zero winding number.
Including the exchange interaction couples |Q; +, + with |Q; −, − , and this idealized picture complicates somewhat.Analyzing the effective two-dimensional interaction between electrons and holes, we find that the exciton eigenstates split into two doublets as with φ Q the polar angle of Q with respect to the x-axis.
and its behavior on the momentum plane is sketched in Fig. 3.We see that it has a winding number of 2 around Q = 0, analogously to the valley pseudospin in Ref. [47].We mention, in passing, that the chiral excitons considered here are of different nature from those observed in the experiment of Kung et al. [30] Indeed, the latter result from high-energy transitions between massive holes and massless Dirac electrons on the surface, whereas the ones in this work arise from long-wavelength transitions between bulk bands near the Fermi level.

Solution of the Bethe-Salpeter equation
define a spin-orbit-parity pseudospin operator σ ≡ (σ x , σ y , σ z ) that acts on the space spanned by |Q; +, + and |Q; −, − .Its expected value for the states |Q; 2 ± is The corresponding pseudospin texture is sketched in Fig. 3.1.We see that it has a winding number of 2 around Q = 0, analogously to the valley pseudospin in Ref. [100].
Note that the wave functions in |Q; 2 + are in general different from those in |Q; 2 − .This is to say that the states |Q; +, + and |Q; −, − are not actually well-defined on their own in the presence of an exchange interaction, which couples them and modifies them.However, the expressions above still make sense for small Q, and more generally if we take into account that the wave functions depend on the eigenstate itself.

Solution of the Bethe-Salpeter equation
In this section we discuss the solutions to the BSE (3.64).The numerical method we have employed is based on Ref.
[99]; a brief overview is given in App.3.B.We have only considered angular momenta |m| ≤ 3 after verifying that higher-m corrections do not play any significant role for the low-lying states.

Labeling of the eigenstates
Before discussing the results, let us introduce a notation to label the different individual states within the families |Q; 2 ± and |Q; 0 ± belonging to the subspaces with s = t and s = t, respectively.It will be enough to specify two additional quantum numbers in order to identify them unambiguously: a principal quantum number n ∈ {0, 1, 2, . ..}, and the azimuthal angular momentum m at Q = 0, which can be chosen to satisfy m ∈ {−n, . . ., 0, . . ., n}.In this way, the states are labeled in the same way as the orbitals of the 2D hydrogen problem [101].The number m is well-defined due to the rotational invariance of the problem at zero momentum.Away from Q = 0 the states are labeled in the same way by connecting them to their zeromomentum counterparts.This can be done unambiguously for the subspaces with s = t and Figure 3: Behavior of the spin-orbit-parity pseudospin around the origin of momentum.As the polar angle φ Q varies from 0 to 2π, the pseudospin σ of the states |Q; 2 + (left) and |Q; 2 − (right) winds around twice, as its orientation is coupled to that of the total exciton momentum.
We may further label the different individual states in |Q; 0 ± by a principal quantum number n ∈ {0, 1, 2, . . .} and by their at Q = 0 well-defined relative angular momentum m ∈ {−n, . . ., 0, . . ., n}, introducing the notation |Q; 0 ± ; n, m .On the other hand, the individual states in |Q; 2 ± are similarly labeled by the at Q = 0 well-defined angular momentum m of the |Q; +, + component in their linear combination.

Exciton dispersion relations and wave functions
We have solved the associated Bethe-Salpeter equation with both the Rytova-Keldysh potential and the twodimensional Coulomb potential after neglecting the effects of the quantum confinement in the z-direction.The former is often used to approximate the long-distance behavior of the microscopic Keldysh potential [48].We find that only the Rytova-Keldysh interaction is compatible with the neglect of quantum confinement, so unless otherwise specified all results correspond to those obtained with this potential.Figs. 4 shows the dispersion relations of all four excitonic ground states, and also of several excited states that have a nonzero angular momentum at Q = 0.The same is shown in Fig. 5 for the Coulomb interaction, whose weaker nature allows us to better visualize the features of the band structure that are qualitatively independent of the interaction details due to the robustness of the topology.Notice the difference in behavior of the energy between the |Q; 2 ± ; 0, 0 and |Q; 0 ± ; 0, 0 doublets around zero momentum, as the former linearly splits off.This effect is heavily influenced by the topology, as we have checked that the dispersions become strongly parabolic in the absence of a band inversion, with the nonanalytic mode being barely noticeable for all states.The linear dispersion may be analytically understood by expanding the effective exchange potential around this point, which is further analyzed in Ref. [49].Moreover, the effective 2 × 2 Hamiltonian for the |Q; 2 ± ; 0, 0 chiral doublet has been discussed previously [47,50,51] and is further studied below.The energy of the |Q; 0 ± ; 0, 0 doublet follows a quadratic dispersion both at small and large momenta, with a crossover taking place around the minimum of the electron-hole continuum, as the inversion of the bands after that has less of an effect.Interestingly, the states |Q; 2 ± with opposite angular momenta +m and −m are split in energy, with the former having higher energy than that of the latter, as seen in the figures for the solid and dashed lines.This is a consequence of the Berry curvature within the nanosheet [39], which provides an anomalous contribution to the single-particle velocity that is only present in the subspace spanned by the |Q; +, + and |Q; −, − states.Our results agree qualitatively with previous works that perturbatively incorporate the effects of the Berry curvature [37][38][39][40][41][42].However, numerical discrepancies are expected between these references and the present work, as the former all consider the effective-mass approximation.This is not appropriate in our system due to the band inversion, which prevents us from decoupling the relative and center-of-mass motions.Consequently, our exciton spectra at zero momentum deviate from the well-known Rydberg series.
Also, the energies of the higher excited states closely follow the electron-hole continuum, as the excitons are more delocalized in real space and their wave functions become those of a separate electron-hole pair.An important feature of the Rytova-Keldysh spectrum is that all excitons are strongly indirect, which leads to long lifetimes due to the strongly reduced radiative recombination rate.As revealed by the shape of the electron-hole continuum, this is a direct consequence of the band inversion of the underlying single-particle bands.Note, however, that in this work both the electron and the hole are taken to reside in the same layer and are thus still direct in real space.Fig. 6 shows the relative exciton wave function of the states |Q; 2 ± for the ground state and several excited states from Fig. 4 at Q = 0. Note the significantly different behavior of the ground-state wave function compared to that obtained from the hydrogen problem, which is proportional to ((a 0 k) 2 + 1) −3/2 .Furthermore, Figs. 7 and 8 show the relative real-space probability density of the low-lying states |Q; 0 ± and |Q; 2 ± , respectively.For the states |Q; 0 ± at Q = 0 we exploit their opposite-angular-momentum degeneracy to obtain linear superpositions resulting in hydrogen-like orbitals shown in the left column, and note how these become deformed in the direction of a nonzero exciton momentum in the right column due to the breaking of the rotational symmetry.In the case of the states |Q; 2 ± , the splitting between states with opposite m at Q = 0 prevents us from writing down such orbitals, so the first column shows a rotationally invariant probability density.Nevertheless, a nonzero exciton momentum breaks again this rotational symmetry, and the wave functions develop lobes in either the transversal or the longitudinal direction as seen in the second and third columns.

Optical properties
We have also derived selection rules for circularly polarized light in the xy-plane moving in the positive zdirection at Q = 0, where the exchange interaction vanishes and the exciton families |Q; +, + and |Q; −, − become uncoupled.For left-circularly polarized light (with angular momentum m γ = +1) we find that the excitons |0; +, +; n, 0 and |0; −, −; n, +2 are bright, whereas the rest are dark.On the other hand, for rightcircularly polarized light (m γ = −1) the only optically active excitons are |0; −, −; n, 0 and |0; +, +; n, −2 .Note that these results combined are in accord with the time-reversal symmetry.By contrast, the excitons |0; +, − and |0; −, + are all dark irrespective of their angular momentum.These results greatly differ from the situation in ordinary semiconductors, where only the s-wave singlet is bright.

Topological excitonic edge states
For small momenta, the behavior of each pair of |Q; +, + and |Q; −, − excitons can be understood by means of an effective 2 × 2 model.Here we are interested in the effects of their topological properties, and in particular we focus on the s-wave ground state.For our purposes it is enough to study the minimal Hamiltonian Here, J is the exchange coupling   4 and the first available n.The wave functions themselves are obtained by multiplying the magnitude by the corresponding phase e imφ k .The wave functions have a maximum around the momentum at which the energy gap presents a minimum.In particular, the wave function in the s-wave case (red line) significantly differs from that obtained from the hydrogen problem, which is proportional to ((a 0 k) 2 + 1) −3/2 (grey dashed line).Here we have set a 0 = 10/ √ 8π nm 1.99 nm, which is suitable for comparison.The values E m shown in the figure correspond to the excitonic eigenenergies for the same states.The quantities ∆ m are the binding energies of each state, that is, the difference between the electron-hole continuum and E m .and σ the spin-orbit-parity pseudospin, and we have added a possible Zeeman-like term ∆σ z that breaks the time-reversal symmetry.The topological properties become apparent after calculating the winding number of the pseudo-magnetic field b(Q), which is w = sgn ∆.This unit winding when ∆ = 0 can be understood from the fact that the normalized vector b(Q) represents a meron configuration with chirality two.Even though we have neglected higher-order terms in the effective Hamiltonian, exciton edge states can in principle be found with these terms included as well, as argued in the Methods section.The approximation performed here admits a simple analytic solution, whereas the expressions become very cumbersome when quadratic terms are kept.
We consider the setting of two semi-infinite systems defined by y < 0 and y > 0, with Zeeman-like couplings ∆ 1 and ∆ 2 , respectively.We have solved the Hamiltonian for Q x = 0, in which case there are two zero-energy eigenstates that read where η ≡ − sgn ∆ 1 = sgn ∆ 2 .These solutions describe states that are localized around the interface at y = 0 and decay exponentially away from this point.We stress that they only exist when sgn ∆ 1 = − sgn ∆ 2 , as there is no continuous normalizable solution when the signs of these Zeeman-like couplings are equal.Furthermore, these states are robust against perturbations of the gap parameter at either side, and in particular the solution exists even in the limit |∆ 1 | → ∞, even though in this case the continuity condition must be relaxed.In Fig. 9 we have plotted a typical example of the corresponding probability density We briefly comment on the term ∆σ z , as it has been introduced by hand in the effective Hamiltonian.This is necessary because the underlying BHZ model is time-reversal symmetric and thus by itself will not give rise to such a term.However, the topological properties only depend on its sign, not on its magnitude, (n, m) = (0, 0) (0, 0)  and emerge no matter how small ∆ may be as long as it is nonzero.We are therefore allowed to add it to the effective model as a perturbation.Experimentally, this term may be realized by a time-reversal-breaking perturbation such as a small magnetic field, via contact to a thin magnetized layer, or with the injection of a small concentration of magnetic impurities.Furthermore, even though the topological properties imbued by ∆ may seem independent of the state of the underlying BHZ Hamiltonian, one must keep in mind that the effective model for the chiral excitons arises only in the topological regime.Consequently, both ingredients are required for the emergence of topological excitons.
For small but nonzero Q x , perturbation theory yields an equal energy shift for both edge states, namely (n, m) = (0, 0) (0, 0)  The values of n and m are indicated in each plot.We have set Q along the x-direction with Q = 0.7 nm −1 .Thus, the system hosts two chiral modes of excitons at the interface, with a Fermi velocity of magnitude 2|J| and with both states moving in the same direction.Due to the factor of η in the dispersion relation, this direction is directly specified by the bulk winding number at either side of the boundary.Hence, reversing the Zeeman-like coupling at both sides also reverses the direction of the chiral excitons, as expected from the analogy of ∆ with a magnetic field.These results can be directly understood from the bulk-boundary correspondence, which predicts a total number of states equal to the difference in winding numbers at both sides, |w y<0 − w y>0 | = 2, with their direction being determined by the sign of this difference [52].

Discussion
Our results show that, in principle, the topology of the conduction and valence bands is indeed inherited by the Wannier excitons, as the exciton wave function contains electron and hole pseudospin textures with a topologically nontrivial winding.This is in particular true for the excitonic basis states that diagonalize the Wannier problem with only the direct interaction included.As the exchange interaction couples these states the physical picture complicates somewhat, but the nontrivial pseudospin winding remains.Furthermore, the bare two-dimensional interaction is heavily modified by the topological band structure.For the excitonic ground state, this ultimately results in a nonchiral doublet with quadratic dispersion relation at low momenta, and a chiral doublet with one linearly dispersing mode and one quadratic mode.Diagonalization of the resulting low-momentum effective model including a Zeeman-like coupling leads to topological edge states of chiral excitons.Since our work does not impose many restrictions on the model parameters, we expect that it remains valid for other similar topological materials such as Bi 2 Te 3 and Sb 2 Se 3 by using the appropriate parameter values.
A future direction of research is to consider the motion of the excitons under the effects of lattice strain, which yields a confining potential in real space.In this scenario we expect that the nontrivial Berry curvature, which constitutes a momentum-space magnetic field, will give rise to an anomalous Hall effect [53,54].Experimentally, we want to resolve the linear dispersion in the chiral doublet, which would be a first step towards the observation of the topological properties of excitons.This is in principle possible by means of terahertz spectroscopy, as the details of the dispersion affect the chemical equilibrium between excitons and free charges in such pump-probe experiments [55][56][57].We also expect the polarizability of the excitons to be strongly affected by their topology, which can again be observed in terahertz conductivity measurements.In particular, the polarizability of the topological s-wave excitons should be significantly reduced with respect to that of their trivial counterparts due to their modified relative wave function.Another interesting feature of the obtained excitonic band structure is that it is indirect as a consequence of the band inversion of the single-particle Hamiltonian, which leads to a greatly reduced rate of radiative recombination processes and thus to long-lived excitons.
The selection rules we have derived can be understood by noting that the eigenstates |χ µ k are also eigenstates of the total angular momentum operator j z = l z + 1 2 s z , where s z is the Pauli matrix acting on the spin part of the single-particle basis and l z = −i∂ φ k .Consequently, the single-particle states |χ c,s k and |χ v,t k have angular momenta s 2 and − t 2 , respectively.By noting that the angular momentum of a hole is opposite to that of the destroyed valence electron, one finds that an exciton can couple to a left-circular photon if its relative angular momentum is 0 and the underlying electron and hole lie in the bands with s = t = +1, or if its relative angular momentum is +2 and the particles lie in the bands with s = t = −1, as we have seen.The argument follows analogously for right-circular photons.
Finally, the (top and bottom) surface electronic states, which have been neglected in this work, may contribute to several aspects that will be analyzed in a follow-up article [58].First, by screening the interaction between bulk electrons and holes, which may be treated as static screening of the Coulomb or Rytova-Keldysh potentials.We have in first instance neglected this effect because of the reduced density of states around the Dirac cone, to which the screening length will be inversely proportional, at least in the Thomas-Fermi regime.Note that we are always considering single-exciton excitations in the material, so that electronelectron correlation effects are effectively incorporated in the parameters of the band structure, for instance via GW corrections.Second, by causing a nonzero transition probability from bulk states to surface states mediated by surface plasmons, which thus in principle provide a mechanism for surface electron-hole decay of bulk excitons.The effects of both the screening and exciton decay will be analyzed perturbatively to verify the correctness of our results, as it depends on the magnitude of the energy shifts and spectral broadenings induced by these phenomena.Away from zero chemical potential, we expect quasi-2D nanosheets of Bi 2 Se 3 with around 6 QLs to provide a promising platform for exciton-plasmonics, which will constitute the main topic of the sequel.Third, by a contribution to the complex conductivity measured in pump-probe terahertz conductivity experiments, which may be modelled as a separate contribution from that of the bulk states, as similarly done for unbound charges in Ref. [55].However, we expect the contribution from the surface states to the absorption to be small compared to that of the bulk states, as the oscillator strength of the latter is much larger than that of the former.

Band structure of Bi 2 Se 3 nanosheets
Our starting point is the k • p Hamiltonian derived in [59] and [60] for modeling three-dimensional Bi 2 Se 3 and Bi 2 Te 3 around the Γ point in the Brillouin zone.To account for the quantization in the z-direction in a nanosheet geometry of thickness L z , we solve the model at the two-dimensional Γ point (k x = k y = 0) with the substitution k z → −i∂ z and hard-wall boundary conditions as done by [61] and [62].Everywhere in our numerics we consider L z = 6 nm, which in the case of Bi 2 Se 3 is sufficiently large for the single-particle Dirac states on the opposite surfaces not to be gapped out by tunneling processes and the nontrivial topology of the bulk to survive [63][64][65], but also sufficiently small for bulk electrons and holes to still behave as in two dimensions.We then project the 3D Hamiltonian onto the energetically highest-lying valence and lowestlying conduction subbands in order to integrate out the z-dependence.Hence, we assume that the relevant low-energy physics of individual particles is confined to this subspace, which has been verified a posteriori by comparing the obtained exciton binding energies to the energy splitting of the bulk subbands due to the confinement in the z-direction.
Ultimately, after projecting the three-dimensional k • p Hamiltonian onto the single-particle ground states of the 3D model, we obtain our desired nanosheet Hamiltonian.Written in terms of Pauli matrices in spin and orbital space, s and τ respectively, it is given by the 4 × 4 matrix where k . Note that products of matrices in different spaces, e.g., s x τ x , are Kronecker products and not matrix products, and that identity matrices are implied.Our effective two-dimensional Hamiltonian is equivalent to that of the BHZ model after a suitable unitary transformation.Furthermore, as expressed in Eq. ( 6) it has the same form as the three-dimensional Bi 2 Se 3 Hamiltonian with k z = 0 and some renormalized values of the parameters with respect to those given by [59] However, the basis in which it is expressed differs from that of the 3D model, since the original Bi + and Se − orbitals have hybridized in the corresponding eigenstates.We denote these hybridized orbitals by Bi + and Se − .The renormalized values of the parameters, with which all of our numerical results have been obtained, are A 2 = 0.41 eV nm, M = 0.28 eV, B 2 = 0.473 eV nm 2 , E = −0.0012eV, and D 2 = 0.202 eV nm 2 .
Diagonalization of H 0 (k) leads to the conduction and valence band energies as well as to the topologically nontrivial single-particle states that we use throughout in our description of excitons.As a consequence of time-reversal symmetry in combination with inversion symmetry, the Hamiltonian does not couple the two subspaces {|Bi + ; ↑ , |Se − ; ↓ } and {|Bi + ; ↓ , |Se − ; ↑ }, and the conduction and valence bands are each two-fold degenerate.The corresponding eigenstates are labeled by their so-called spin-orbit parity, defined as the eigenvalue of the operator s z τ z that commutes with H 0 (k).
We next introduce the vector of pseudospin operators Γ ≡ (s x τ x , s y τ x , τ z ), which are used to rewrite Eq. ( 6) as . The nontrivial topology of the Hamiltonian is now very explicit in this form, since d(k) is a skyrmion texture in the momentum plane (k x , k y ) and the operator Γ reduces exactly to the three Pauli matrices in the uncoupled even and odd spinorbit-parity subspaces.Therefore, the expectation value of the pseudospin as a function of k follows the nontrivial winding of d(k) and the winding number of the latter is up to a possible sign equal to the Chern number of the conduction and valence bands.

Exciton basis and wave functions
The conduction and valence states are x; a|k; µ = (e ik•x / √ V ) a|χ µ k , where V = L x L y , the states |a denote our four-dimensional combined spin and orbital basis states, and |χ µ k with µ ∈ {c, ±; v, ±} are the eigenstates of the Hamiltonian in Eq. ( 6).Note that the latter are not periodic due to the fact that our Hamiltonian is derived from k • p theory and only accurately models the band structure around the Γ point.An exciton state X ≡ {Q, s, t} is labeled by the total momentum Q and the spin-orbit-parities s and t of the electron and hole bands, respectively, as well as by an additional set of ro-vibrational quantum numbers that we introduce later.Such an exciton state is a bound state in the polarization [66], which in second quantization is determined by the pair correlation function where s, t ∈ {+, −} and R ≡ (x 1 + x 2 )/2 is the position of the exciton.Here Φ s,t Q,k is the relative wave function of the exciton in momentum space.One may also perform a particle-hole transformation so that holes become positive-energy excitations.In terms of electrons and holes, the conduction states remain unchanged, a|χ e,s q = a|χ c,s q .By contrast, the hole states satisfy b|χ h,t −q = χ v,t q |b .In this picture the pair correlation function thus reads If desired we can then also obtain a first-quantized wave function, which is given by This is normalized provided that 1 Note that interchanging both particles, thus x 1 ↔ x 2 and (a, 1) ↔ (b, 2), results in an overall minus sign after we perform the transformation k → −k.
The single-particle states, together with the relative momentum states |k rel , describe the exciton states |Q; s, t fully in Dirac notation as where R|Q = e iQ•R / √ V .The four combinations of the spin-orbit-parity labels s and t lead to four distinct exciton basis states, similarly to the singlet and triplet excitons in regular semiconductors.As explained in the next section, the states |Q; s, t in Eq. ( 11) are exact eigenstates of the Wannier problem with only the direct interaction included.This is in fact a much used approximation in the literature [43], but in our case not sufficiently accurate as the exchange interaction in principle couples these states for Q = 0. Nevertheless, the above set of states can still be considered as the most appropriate basis for the full excitonic problem.Further note that, as advanced, |Q; s, t more precisely stands for an entire family of ro-vibrational states which must be labeled by additional quantum numbers describing the different relative wave functions that solve the exciton problem.
The topology of these exciton basis states can now be intuitively understood, as in the single-electron case, from the dependence of the expectation value of the pseudospin Γ e (Q) on the momentum Q, obtained from Eq. (11) as and similarly for the pseudospin of the hole.
To determine the topology of the states |Q; s, t mathematically more rigorously, we need to compute the winding number of the pseudospin.Note that for the moment we particularize to the case of a globally vanishing exchange interaction, in which case all states |Q; s, t are uncoupled and therefore have a welldefined Chern number.To access the topological properties we must compute the Berry connection [8][9][10][11][12], which in the excitonic case is given by a sum of three distinct terms, Each of these terms contributes to the local Berry curvature or momentum-space magnetic field through s,t (Q), where i = e, h, Φ.Finally, the Berry curvature is connected to the desired winding or Chern number by where the first surface integral is performed over the first Brillouin zone, which in our continuum model amounts to integration over the infinite momentum space, and the equivalent line integral is performed along a contour with Q → ∞.It is clear that the terms B s,t (Q) not only contain the direct weighted Berry curvature of the underlying single particles, but also an interference term between the single-particle Berry connection and the excitonic envelope wave function.However, the total contribution of these terms to the exciton Chern number is given only by the single-particle electron or hole Chern number.This can be shown, e.g., for the electron case, by first expanding χ e,s . In our model one then observes that all terms in the line integral decay faster with Q than the zeroth-order term χ e,s Q/2 |∇ Q |χ e,s Q/2 and thus vanish on the contour at infinity.Therefore, and the normalization condition 1 may now be used.We are then left precisely with the expression for the free electron Chern number.
There can in principle be an additional contribution to the Chern number due to the envelope wave function itself, given by s,t (Q).However, this can be shown to vanish here by analyzing the winding of the direct interaction around the origin of Q.This winding is trivial, meaning that the potential does not pick up a phase as one circles around this point.That is, , where φ q is the angle between the momentum q and the x-axis.A straightforward analysis of the exciton eigenproblem then shows that we can choose a gauge such that Φ Observing that the system enjoys reflection symmetry with respect to the x-axis when Q = Qx, this may now be used to verify that C Φ = 0.
We thus conclude that all global topological properties of the excitons are introduced by the electron and hole states |χ e,s Q/2+k and |χ h,t Q/2−k .Hence, the total Chern number of the above exciton basis states is C = C e + C h , which for our BHZ model becomes C = s + t by explicit calculation.Physically, this result can be most easily understood by the fact that the even and odd spin-orbit-parity subspaces are related by time-reversal symmetry and that the wave function for the hole is the complex conjugate of the electronic valence-band wave function, as we have seen.Note, however, that the local properties that influence for instance the electron and hole transport are quantified by the Berry curvature, and are thus still dependent on the precise shape of the relative wave functions and the interference terms.
We stress that the intuitive picture given above is only valid in the case of zero exchange interaction.As seen in the following section, when this is included the two subspaces with s = t become coupled together and cannot be treated individually.The Chern number is then technically not well-defined since the time-reversal symmetry protects the degeneracy at Q = 0.However, the nontrivial winding caused by the underlying Chern numbers remains, and as a result the full excitonic eigenstates still possess a chirality of ±2.This is crucially dependent on the fact that the exchange interaction that couples the two subspaces winds nontrivially around the origin of Q, transforming instead as

Electron-hole interaction and Bethe-Salpeter equation
Having introduced the free part of the full Hamiltonian, we now proceed to discuss the electron-hole interaction potential that binds these particles together to form an exciton state.The interaction potential is , where V D and V X denote the direct and exchange interactions, respectively.These are given by where V (q) is the bare electrostatic potential within the nanosheet.A variational approach with the trial wave functions in the previous section leads to the Bethe-Salpeter equation [67,68] k ,s ,t where we have restored the spin-orbit-parity indices in the wave functions.The matrix elements of the total Hamiltonian, including the electron-hole interaction, read We have solved this equation with both the bare two-dimensional Coulomb potential V C (q) with the dielectric constant of the surrounding medium ε s , and with the Rytova-Keldysh potential V RK (q), which also takes into account the dielectric constant ε d of Bi 2 Se 3 and is typically used in a nanosheet geometry.The explicit forms of these potentials in momentum space read where r 0 = (ε d /2ε s )L z is a screening length that depends on the dielectric constants of Bi 2 Se 3 and the surrounding environment.We have chosen the relative permittivity ε s = 6, which is a typical low value for the environment [69].The bulk dielectric constant has been set to ε d = 28 in accord with recent first-principles studies of rhombohedral Bi 2 Se 3 in the near-infrared region at 6 QLs [70,71].
Our approach implies that we neglect the effects of quantum confinement on the classical electrostatic interaction.This is acceptable if the in-plane separation of the bound electron and hole is larger than the nanosheet thickness, as in this case the electric field lines will mostly lie in the surrounding environment.We find that the long-wavelength Coulomb interaction alone cannot accurately model excitons in bismuth selenide nanosheets because the resulting exciton diameters of the low-lying states are in fact smaller than the slab thickness, as seen in Table 1.Coincidentally, the Rytova-Keldysh interaction closely resembles the quantum-confined Keldysh potential at all momenta [72], which is in any case guaranteed to accurately model the electrostatic interaction in this geometry.For this reason, it is not important which of these potentials enters Eqs. ( 18) and ( 19), and we have additionally checked that the difference in binding energies obtained with the Rytova-Keldysh and the quantum-confined Keldysh potentials are no larger than 1 meV, leading to errors lower than 10%.Note, however, that fully incorporating the quantum confinement would in principle lead to electrostatic interactions between subspaces with different spin-orbit parity due to the small overlap of the wave functions in the z-direction.This may have a small effect on the Rytova-Keldysh ground-state excitons, which are slightly smaller than the nanosheet thickness, but we note that the dielectric constant we have employed for the surrounding environment is relatively low and using a higher value would also lead to bigger excitons and thus mitigate this issue.
Additionally, we must compare the binding energies of the excitons with the splitting of the first two bulk subbands due to confinement, as we have neglected the subspaces of higher energy.Our approximation of projecting the 3D Hamiltonian onto the bulk subbands closest to the Fermi surface will be justified if the former is smaller than the latter.A binding energy larger than the bulk-subband splitting would imply the need to include the wave functions of higher excited states and thus effectively hinder the two-dimensional treatment of the problem.The splittings of the conduction and valence subbands for L z = 6 nm are found to be about 61 meV and 24 meV, respectively, and the binding energies are given in Fig. 6.Except for the ground state, all states have a binding energy that is smaller than the relevant subband splittings.In the case of the ground state, the binding energy is only around 1 meV larger than the valence-band splitting, so we do not expect that including the second subband will lead to significant modifications.We conclude that the Rytova-Keldysh potential is apt for our study of excitons in Bi 2 Se 3 nanosheets, especially if we use a higher relative permittivity for the environment.Hence, unless otherwise specified, all figures in this article correspond to results obtained via the Rytova-Keldysh potential.Note, however, that the topology is neither affected by the values of the dielectric constants, nor by the precise form of the interaction potential.1: Mean exciton diameters of the zero-momentum states |0; +, +; 0, m and |0; −, −; 0, −m as obtained with the Coulomb and Rytova-Keldysh potentials with ε s = 6 and ε d = 28.In the case of the Coulomb potential the radii must be compared to the film thickness, which is L z = 6 nm.
Due to the orthogonality of the states in the different spin-orbit-parity subspaces the exchange interaction V X only contributes when the spin-orbit parities of the electron and the hole are equal in both the initial and final states and when Q = 0, as the inner products in Eq. ( 19) tend linearly to zero as a function of Q.Thus, specifically, the states |Q; +, − and |Q; −, + are not coupled by the exchange interaction, whereas the states |Q; +, + and |Q; −, − are.We can therefore solve the problem in the subspaces spanned by |Q; +, + and |Q; −, − on the one hand, and by |Q; ±, ∓ on the other hand.Since the interaction potentials are, up to complex conjugation, the same for the states |Q; +, − and |Q; −, + , the corresponding wave functions also only differ by complex conjugation, and both states are degenerate in energy for all Q.By contrast, the eigenstates spanned by |Q; +, + and |Q; −, − are degenerate in energy only for Q = 0.

Derivation of optical properties
The oscillator strength of an exciton |0; s, t; n, m reads [73,74] where ê is the Jones vector of the outgoing beam, ε m s,t;n is the energy of the zero-momentum exciton, and v(k) is the velocity operator.Since only interatomic transitions are expected to play a role in Bi 2 Se 3 , the velocity operator is well approximated by v(k) ≈ ∇ k H 0 (k) [75].The Jones vectors of left-and right-circularly polarized light are ê+ and ê− , respectively, with ê± = 1 √ 2 (x ± iŷ).Computing the matrix elements that enter the above equation readily yields the reported results.

Effective exciton Hamiltonian
The small-Q behavior of excitons can be understood by constructing an effective model for the pairs of states |Q; +, +; n, m and |Q; −, −; n, −m that are degenerate at zero momentum, as already done before [47,50,51].By computing the relevant matrix elements with the wave functions for Q → 0 we obtain the effective Hamiltonian where X = (n, m) labels the particular exciton doublet.Here, ω X is the energy of the doublet at zero momentum, M X is an effective mass, and We note that these parameters depend on which particular doublet the model attempts to describe.According to Eq. ( 26), the states with odd angular momentum at Q = 0 should remain degenerate for small Q, which is indeed the case in Figs. 4 and 5.However, this degeneracy is broken to third order in Q, but this effect is not captured by our perturbative scheme and requires including corrections to the zero-momentum wave functions as well.On the other hand, the states with even angular momentum may split into a linear mode and a quadratic mode, as is the case for the lowest s-wave and d-wave states in the figures.This effect is expected to be most important for states with m = 0, as the lowest-order contribution to the splitting parameter J X is proportional to |Φ Q (r = 0)| 2 , which is nonzero only for s-wave excitons.
To solve this model we restrict ourselves to s-wave excitons and neglect the terms proportional to the identity matrix, which do not contribute to the topological properties.Also, we drop the remaining terms proportional to Q 2 , as we are interested in the low-energy physics.We introduce a small Zeeman-like term ∆σ z that breaks the time-reversal symmetry and makes the model gapped, resulting in H eff (Q) = b(Q) • σ with b(Q) = (JQ cos 2φ Q , JQ sin 2φ Q , ∆).The winding number of the Hamiltonian is computed via the formula and gives w = sgn ∆.The ∆ term also introduces Berry curvatures Ω ± xy (Q) = ±J 2 ∆/ ∆ 2 + J 2 Q 2 3/2 for the positive-and negative-energy solutions, which give Chern numbers C = ± sgn ∆.
To solve the model with a position-dependent Zeeman term we must perform the substitution Q y → −i∂ y in the effective Hamiltonian due to the breaking of translational invariance in the y-direction.This introduces some complications, as the object Qe 2iφ Q does not have a unique operator generalization.However, for Q x = 0 we unambiguously have e 2iφ Q = −1, whereas Q reduces to √ −∂ 2 y .We interpret this as an operator whose square is −∂ 2 y and whose eigenfunctions are plane waves.However, there are still two possibilities for the action of this operator on a plane wave, namely √ −∂ 2 y e Λy = ±iΛe Λy .These two choices precisely give the two modes we are looking for.For Q x = 0, the model to solve is thus and we look for solutions with E = 0 consistent with the sought-after chiral excitons.For each choice of operator action, the ansatz Ψ 1 (y) = e Λy and Ψ 2 (y) = χΨ 1 (y), with χ ∈ C, gives a single solution of the Jackiw-Rebbi type [8,76].In this way we find the two linearly independent eigenfunctions of Eq. ( 4).For small but nonzero Q x we can use perturbation theory to find the first-order correction to the spectrum by using Qe , where in analogy to the above we define 1 ∂y e Λy = 1 Λ e Λy as the action of the inverse derivative operator on plane-wave solutions.
This analytic solution has been found by neglecting the terms quadratic in Q in the effective exciton Hamiltonian.However, we expect the presence of exciton edge states even when including these higher-order terms, although an explicit calculation is much more involved.Firstly, we note that the inclusion of the Zeeman-like term will produces a globally gapped exciton spectrum.This is true not only for the effective model with or without quadratic corrections, but also for the full excitonic spectra of Figs. 4 and 5.The presence of edge states is then inferred directly from a nonzero Chern number.Including the quadratic terms yields the following Berry curvatures for the valence and conduction bands: Integration of this expression yields C ± = ± sgn ∆, the same as before, so edge modes are expected in this full model as well.
As an additional check, we further note that in this more general case the Berry curvature is nonzero even for J = 0. Hence the topology remains even in this case, for which the Hamiltonian is parabolic and only contains integer powers of Q x and Q y .Thus, we may use the argument of Ref. [52] to infer the existence of topological edge states.The Hamiltonian can be written as H(Q) = b 0 (Q) + b(Q) • σ, where for each fixed Q x the vector b(Q) is a parabolic function of the perpendicular momentum Q y .For all momenta Q y = 0 the origin lies on the concave side of the projection of this parabola onto the plane, and so there should be edge states, which will be chiral due to the nonvanishing topological invariant.

Computational details
The Bethe-Salpeter equation has been solved with an independently developed MATLAB code that implements the discretization procedure of Ref. [77] after rewriting all equations in terms of the dimensionless momentum u = kL z and keeping only contributions from angular momenta |m| ≤ 3. The integrals in the matrix elements have been performed up to a momentum cutoff U = 10 and with a discretization ∆u = 0.05, and we have verified that the results do not depend on the cutoff by performing additional calculations up to U = 40.The long-wavelength divergence of the Coulomb potential has been regularized via an infrared cutoff ∆ V u chosen such that ∆ V u/∆u ≈ 0.2262, for which we have checked that the energies do not depend on ∆u.If ∆ V u/∆u is chosen differently, a well-defined extrapolation of the exciton energy levels for ∆u → 0 can be done by using different values of the discretization.
In order to identify the angular momentum quantum numbers unambiguously after solving the Bethe-Salpeter equation, we choose a nonsingular gauge for the single-particle eigenstates that enter the direct and exchange potentials [39,73].We have checked that with this choice the level ordering in the trivial regime reduces to that of the 2D hydrogen atom.

1− 1 Figure 1 :
Figure 1: Band structure of two-dimensional Bi 2 Se 3 around the Γ point.Both the conduction bands and the valence bands are shown in solid lines.Each is two-fold degenerate, and due to spin-orbit coupling features a band inversion and an avoided crossing of the corresponding uncoupled states depicted with dashed lines.The dispersions are isotropic in the k • p approximation we are using.

Figure 3 . 1 .
Figure 3.1.Orientation of the spin-orbit pseudospin vector σ for the states |Q; 2 + (left) and |Q; 2 − (right) as the polar angle φ Q varies from 0 to 2π.The pseudospin winds twice around itself in this process.

Figure 5 :
Figure 5: Excitonic dispersion relation for the Coulomb potential.All colors and line types are equivalent to their counterparts in Fig. 4.

Figure 6 :
Figure 6: Magnitudes of momentum-space exciton wave functions.These correspond to the zeromomentum states |0; +, +; n, m and |0; −, −; n, −m and are shown for several values of m as ordered in Fig.4and the first available n.The wave functions themselves are obtained by multiplying the magnitude by the corresponding phase e imφ k .The wave functions have a maximum around the momentum at which the energy gap presents a minimum.In particular, the wave function in the s-wave case (red line) significantly differs from that obtained from the hydrogen problem, which is proportional to ((a 0 k) 2 + 1) −3/2 (grey dashed line).Here we have set a 0 = 10/ √ 8π nm 1.99 nm, which is suitable for comparison.The values E m shown in the figure correspond to the excitonic eigenenergies for the same states.The quantities ∆ m are the binding energies of each state, that is, the difference between the electron-hole continuum and E m .

Figure 7 :
Figure 7: Relative real-space probability densities of the exciton eigenstates |Q; 0 ± ; n, m .The rows correspond, from top to bottom, to the first five dotted-line states of Fig. 4. In the left column, Q = 0 and we take appropriate linear combinations resulting in hydrogen-like orbitals.The right column corresponds to Q along the x-direction, with Q = 0.7 nm −1 .The values of n and m are indicated in each plot, as well as the orbital name in the standard notation.

Figure 8 :
Figure 8: Relative real-space probability densities of the exciton eigenstates |Q; 2 ± ; n, m .The rows correspond, from top to bottom, to the first five solid-and dashed-line states of Fig. 4.More precisely, the states shown are |0; 2 ± ; n, m (left column), |Q; 2 − ; n, m (middle column), and |Q; + ; n, m (right column).The values of n and m are indicated in each plot.We have set Q along the x-direction with Q = 0.7 nm −1 .

Figure 9 :
Figure9: Probability density of chiral exciton edge states.Here, Q x = 0, and we have used the parameters ∆ 1 = 3∆ 2 and |∆ 2 /J| = 1.The states are localized at the interface y = 0 and their wave functions decay exponentially with the distance from this point.
The wave functions contained in |Q; +, − and |Q; −, + are related by complex conjugation, and so are those contained in |Q; +, + and |Q; −, − .Due to the combined phase factors e ±iφ Q , the members of the doublet |Q; 2 ± have chirality two, a signature of the nonzero Berry curvature within the nanosheet.This chirality can be understood via the introduction of a spin-orbit-parity pseudospin operator σ ≡ (σ x , σ y , σ z ), where σ i are the Pauli matrices acting on the space spanned by |Q; +, + and |Q; −, − .Its expected value for the states |Q; 2 ± is Excitonic dispersion relation for the Rytova-Keldysh potential.Shown are the exciton eigenstates |Q; 0 ± (dotted lines), |Q; 2 + (solid lines), and |Q; 2 − (dashed lines), with Q the total exciton momentum.Different colors correspond to different values of m at Q = 0 as indicated in the legend.The 1s state (not shown) has higher energy than the 1p, 2d and 3f states.The excitons |Q, 2 ± with opposite angular momenta are split in energy, with the higher state having angular momentum m > 0 and the lower state m < 0. The magenta region at the top represents the electron-hole continuum, delimited by a solid purple line resulting from analytically minimizing the energy gap with respect to the relative electron-hole momentum in the absence of interactions.The purple dotted line represents the continuum threshold as calculated numerically including only angular momenta |m| ≤ 3 and tends to the solid line upon inclusion of higher values of |m|.