Topologically Protected Photovoltaics in Bi Nanoribbons

Photovoltaic efficiency in solar cells is hindered by many unwanted effects. Radiative channels (emission of photons) sometimes mediated by nonradiative ones (emission of phonons) are principally responsible for the decrease in exciton population before charge separation can take place. One such mechanism is electron–hole recombination at surfaces or defects where the in-gap edge states serve as the nonradiative channels. In topological insulators (TIs), which are rarely explored from an optoelectronics standpoint, we show that their characteristic surface states constitute a nonradiative decay channel that can be exploited to generate a protected photovoltaic current. Focusing on two-dimensional TIs, and specifically for illustration purposes on a Bi(111) monolayer, we obtain the transition rates from the bulk excitons to the edge states. By breaking the appropriate symmetries of the system, one can induce an edge charge accumulation and edge currents under illumination, demonstrating the potential of TI nanoribbons for photovoltaics.


I. METHODS
The material of choice here is a monolayer of Bi(111), also known as beta-bismuthene, which was predicted [1,2] and observed to be a 2D TI [3,4].We use a Slater-Koster tightbinding model which properly describes the electronic structure of Bi(111) [5].Monolayer Bi(111) appears as a puckered honeycomb lattice and presents two possible standard edge terminations: zigzag and armchair [2].We consider ribbons with zigzag terminations, defining N as the number of Bi dimers across the unit cell.In particular, we choose N to be even so that the space group of the ribbon is symmorphic [6] to restrict the symmetry discussion to the point group.Nevertheless we do not expect this to be relevant for the calculations forward.
With the material model established, we now turn our attention to the description of excitons.As they are bound electron-hole pairs, one must consider an electrostatic interaction between the electrons of the solid, i.e.: ε nk denotes the bands of the ribbon, and V ijlm are the matrix elements of the electrostatic interaction.The indices i, j, l, m are short-hand notation for pairs of quantum numbers (n, k), with n the band index and k the crystal momentum, which is an scalar since the system is one-dimensional.The exciton states are then given by a superposition of electronhole pairs between valence and conduction bands: where |Ω⟩ denotes the Fermi sea, which we assume is the ground state of the system, i.e. we work in the Tamm-Dancoff approximation [7].Q denotes the center-of-mass or total momentum of the exciton.For the bulk excitons we discuss throughout the article, we restrict v and c to bulk bands, i.e. all bands except for the edge ones (those closing the gap).Then, the problem of determining the exciton coefficients A Q vc (k) amounts to diagonalizing the total Hamiltonian (1) projected over the sector of single electron-hole pairs, i.e.P HP .Therefore, we need to obtain the matrix elements of the interacting Hamiltonian represented in the basis of electron-hole pairs.Denoting the electron-hole pairs as |v, c, k, Q⟩ = c † ck+Q c vk |Ω⟩, the matrix elements ⟨v, c, k, Q|H|v ′ , c ′ , k ′ , Q⟩ ≡ H cc ′ vv ′ (k, k ′ , Q) are given by: The D, X terms represent the direct and exchange interaction between the electron and the hole, and are given by some specific matrix elements of the two-body interaction V : This approach to excitons has been used previously in several works [8][9][10][11], since it is a simpler and faster alternative to many-body perturbation theory.To compute the excitons, we need to evaluate explicitly the direct and exchange terms (4) in terms of the tight-binding Bloch states.If {C nk iα } iα are the coefficients corresponding to the eigenstate |nk⟩ in the lattice gauge [12], then the interaction matrix elements D can be written as: where R are Bravais vectors, R = |R| and t i are the atomic positions of the motif.V ij denotes a lattice Fourier transform of the potential centered at t j − t i .This expression is obtained evaluating the four-body integral of the interaction in real-space, using the tight-binding approximation (i.e.point-like orbitals) [13].It is also possible to derive an alternative expression for the interaction matrix elements with the interaction in reciprocal space through its Fourier series [14,15].Given that the platform for our excitons is a semi-infinite ribbon, the real-space approach is better suited since it takes into account the finite boundaries.As for the exchange term X, we set X = 0 assuming that its contribution is negligible.
For the electrostatic interaction, we use the Rytova-Keldysh potential [16,17].This model has gained interest to describe excitons in two-dimensional materials since it includes screening, as opposed to the bare Coulomb interaction.This is particularly relevant since the exact diagonalization approach does not screen the bare interaction, while manybody perturbation theory does, usually through the random-phase approximation [18].The Rytova-Keldysh potential is given by: where ε is an environmental dielectric constant, given by ε = (ϵ s + ϵ v )/2, with ϵ s , ϵ v the dielectric constants of the substrate and vacuum respectively.r 0 is a screening length, such that r 0 = dϵ/(ϵ s + ϵ v ), with d the height of the layer and ϵ the dielectric constant of the material of interest [19].H 0 , Y 0 are Struve and Bessel functions of second kind respectively.
We follow the prescription of [8] to renormalize the divergence of the potential at r = 0, by setting V (0) = V (a), a being the lattice parameter of Bi(111).For our monolayer Bi(111), we set the dielectric constant to ϵ = 40, and the environmental one to ε = 2.45, corresponding to a SiO 2 substrate.
The edge electron-hole pairs are defined as conventional pairs, but instead we specify the boundary rather than the band number: where c † sk+Q creates a conduction electron such that it is located at side s with the specified momentum.The analogue is done with c s ′ k for the valence hole, with the overall restriction that the energy of the pair is the same as that of the exciton.Note that for complex edge bands, there could be several electron-hole pairs verifying the same conditions, and so all of them should be taken into account in the transition rates.The calculation of the transition rates is done with Fermi's Golden rule, which involves the matrix elements of the electrostatic interaction between the exciton and the final edge electron-hole pair.Expanding the exciton into electron-hole pairs, the transition rates can be cast in terms of the direct and exchange terms of eq. ( 4).Since we neglect the exchange term, the final expression for the rates is: The calculation of the velocity is based on the second quantized version of the velocity operator, namely v = ij v ij c † i c j where the indices i, j correspond to pairs (n, k).Taking the expected value with the edge electron-hole pair yields the following expression: this is, the total electronic velocity is given by the velocity of the electron minus the velocity of the hole [20].For the computation of the single-particle velocity elements we refer to [12].
Note that as with the electron-hole pair, the velocity of the exciton can be separated into contributions from the electron and the hole separately.In this context, the velocity operator can be rewritten as: where we have defined two new operators.This leads to defining the center-of-mass velocity of the exciton as: which is the relevant quantity for determining the propagation direction of the excitons.
All the quantities shown in the article, from the exciton spectrum to the probability densities and the transition rates have been calculated using the Xatu code [13].For the calculation of excitons we have used ) is the number of bulk valence (conduction) bands taken into account, and a minimum of N k = 800 points in the Brillouin zone to ensure convergence of all transition rates.
All the excitons have been computed using N v = N c = 4 and a minimum of N k = 800, while in some cases the number of k points was even higher to ensure that the lowest magnitude transition rates were properly converged, up to N k = 1200.Here we show additional results relative to the convergence of both the excitons and the transition rates with the number of k points, as well as the effect of the number of bands on the results of the main text.
The exciton energies converges relatively quickly with N k , requiring around N k ∼ 200 to be converged.On the other hand, the number of bands included in the BSE is particularly relevant since we are working with a ribbon, i.e. a semi-infinite system where the previous degree of freedom k ⊥ of the periodic system now corresponds to additional bands in the system, which are more closely packed.Therefore, in principle one should include a number of bands in the ribbon equivalent to the number of k points one would use for the 2D periodic system.However, as we will see, the convergence of transition rates with N k imposes a computational restriction on the number of bands we can consider, given that the BSE matrix size scales as When computing the exciton spectrum with Q = 0 and without spin-orbit coupling (SOC), we observe a four-fold degeneracy of all energy levels.This is to be expected, since the point group of the ribbon is C 2h , whose irreducible representations (irreps) are all of dimension 1.For a spinless system, all states would be non-degenerate, while for a spinful system with neither SOC nor exchange, four-fold degenerate.As we turn on the SOC, the exciton spectrum acquires a fine structure.The character table of the C 2h double group shows one-dimensional irreps [21].However, we observe two-fold degeneracy for some states, in particular for the ground state.This is produced by the presence of inversion symmetry which results in degenerate bands, together with time-reversal symmetry even in presence of spin-orbit coupling.The introduction of the sublattice potential, which breaks inversion symmetry, lifts the quasiparticle band degeneracy and also results in a splitting of the degenerate exciton bands.
The real-space electronic density probability of the excitons, is computed from the following expression: which is evaluated using the tight-binding approximation.In Fig. (S3) we show the wavefunctions with and without SOC.Without SOC we observe the standard hydrogenic spectrum, starting from s states, then p and so on.Instead, for the topologically insulating state the ground state is now p-like along the finite direction, instead of the periodic one.This appears to be an effect of the confinement of the wavefunctions, as the same computation for bulk Bi(111) shows that the ground state for both the trivial and topological phases are s-like.As it is shown in the section for Q ̸ = 0 in the main text, we also characterize the total spin projection of the ground state exciton as well as the center of mass velocity of the exciton, which are both calculated the same way.The total spin projection is simply given by the the second quantized version of the S z operator: where i, j denote pairs of quantum numbers i ≡ (n, k), and σ ij = ⟨i|S z |j⟩.Taking the expectation value with an exciton state, one obtains: which can be regarded as the expectation value of the spin taken over the conduction electron plus the expectation value taken over the valence hole.We have defined here new spin operators associated with the electron and hole, respectively: into which S z can be split.Analogously, to compute the total electronic velocity of the exciton, one generically defines the velocity in second quantization: transition rate, the higher is the number of k points required to converge the plot.It can be seen that the rates oscillate as they convergence; if these oscillations are smaller than the final rate, it converges quickly.Otherwise, for smaller rates, these oscillations can be of the same order of magnitude or higher, preventing the rate from converging until they are sufficiently dampened.For N = 14, we see that around N k ∼ 800 the rates are properly converged, and thus is the baseline we have taken in the calculations shown in the main text.If the width of the ribbon, or the different parameters of the model change, the rates might diminish even further, in which case more k points are added to the calculations.
Finally, we show the plots of the main text but for N v = N c = 2, for both Q = 0 and Q ̸ = 0.While the results are slightly different quantitatively, the same discussion holds.
We observe that in presence of the edge offset potential, a charge imbalance can be created at the edges of the material from the dissociation of excitons into edge electron-hole pairs.
Analogously, for finite Q a net current can form at the edges of the material as seen from all the transition rates.In the same way, it is possible to engineer the different parameters of the system to tune the rates.For completeness, we also study the effect of the edge termination of the ribbon on the transition rates.The main effect stems from the different band structure, which for the armchair termination results in more complex edge bands.As shown in Fig. S12(a), for the armchair ribbon there are 4 edge bands, which in presence of the edge offset potential results in 8 different bands.Therefore, for one given rate Γ ± ss ′ with fixed sides s, s ′ there can be multiple electron-hole pairs available.For simplicity, we only compute the transition rates to the electron-hole pair that gives the dominant contribution to each rates, namely the pairs closer to the Γ point.These rates are shown in Fig. S12(b) as a function of the width of the armchair nanoribbon for Q = 0 and N vc = 4.We observe the same behaviour as with the zigzag, except for the fact that the domining rate now is Γ LR instead of Γ RL .This is attributed precisely to the more complex band structure, since for the exciton energies present there are two available pairs for the Γ LR rate: one closer to K and another one closer to Γ in the BZ (the latter shown in Fig. S12(a)), while in the zigzag ribbon there was only one close to K. Thus, even though the dissociation direction is reversed, the charge separation mechanism still takes place.Likewise, if we were to consider the rates as a function of Q, we would expect the same behaviour as in the zigzag, and more importantly, the same direction for the currents given the slope of the bands for the relevant edge electron-hole pairs.
We assume that the ground state exciton consists of a degenerate subspace of N states {|X n ⟩} n .Then, since those excitons have Q = 0, the time-reversal operator T will map those excitons onto themselves, i.e.: where U is the unitary matrix of the representation.Then, with the transition rate defined as the sum over all degenerate states, we may use that the interaction is time-reveral invariant, [V, T ] = 0 to prove the identity: • Estimation of the photocurrent produced We can estimate the current with j = ηn(Q)ev, where n(Q) is the linear density of excitons available for dissociation, η is the efficiency of the effect, e the charge and v the velocity of the final electron once the exciton has dissociated.First, we assume that the linear density of excitons available for dissociation is steady and that the laser is intense enough so that there is always at least one exciton available in the bulk for dissociation.Then, the efficiency η is defined as η = Γ − RL / ss ′ Γ ± ss ′ , where Γ − RL is the dominating rate relative to the edge current generation.Thus, ηn(Q) is the fraction of edge electron-hole pairs relevant to the effect.Since only one exciton can dissociate at a time (since the edge bands are saturated if there is already one electron-hole pair), we set n(Q) = 1.Doing the calculation for N = 14, Q = 0.1 as in Fig. 5, we estimate a current of j ≈ 0.5 µA.Summing over Q would yield a higher current, but again in the order of µA.
FIG. S1.(a) k probability density of the first four exciton levels for N v = N c = 2. (b) Convergence of the ground state exciton energy with N k and the number of bands included for both valence N v and conduction N c , i.e.N c = N v = 2, 4 or 6

FIG
FIG.S3.Real-space electronic probability density of the ground state exciton (left) and first excited state (right).The top row densities correspond to the system with SOC, while the bottom ones are computed without SOC.

FIG
FIG. S5. (Left) Total expected spin projection S T z of the exciton as a function of Q for N v = N c = 2, N k = 200, for different values of the staggered potential V st .(Right) Total exciton spin S T z for Q = 0.1 as a function of the staggered potential V st , for N v = N c = 2, N k = 200.
FIG. S7.Total velocity of the final electron-hole pair, and velocity of the individual components for N v = N c = 2, N k = 300.
FIG. S8.Convergence of the transition rates Γ ± ss ′ with N k .(Left) Rates for a ribbon of width N = 12, with N v = N c = 2 and Q = 0. (Middle) Rates for a ribbon of width N = 14, for N v = N c = 2 and Q = 0.1.(Right) Transition rates for a ribbon of width N = 14, for N v = N c = 4 and Q = 0.1.