Excitonic properties of hBN from a time-dependent hartree-fock mean-field theory

Francisco Ricardo Lobo Ribeiro Excitonic properties of hBN from a time-dependent Hartree-Fock mean-field theory august 2023 UMinho | 2023 Francisco Lobo Excitonic properties of hBN from a time-dependent Hartree-Fock mean-field theory University of Minho School of Sciences University of Minho School of Sciences Francisco Ricardo Lobo Ribeiro Excitonic properties of hBN from a time-dependent Hartree-Fock mean-field theory Masters Dissertation Master’s in Physics Dissertation supervised by Doctor Bruno Amorim Doctor Nuno Peres august 2023 Copyright and Terms of Use for Third Party Work This dissertation reports on academic work that can be used by third parties as long as the internationally accepted standards and good practices are respected concerning copyright and related rights. This work can thereafter be used under the terms established in the license below. Readers needing authorization conditions not provided for in the indicated licensing should contact the author through the RepositóriUM of the University of Minho. License granted to users of this work: CC BY-NC-SA https://creativecommons.org/licenses/by-nc-sa/4.0/ ii Acknowledgements First and foremost, I show my appreciation to University of Minho, in particular to the Department of Physics, for the institutional support and, most importantly, for its amazing and kind professors who guided me not only through my Bachelor’s and Master’s degree, but also through life. Moreover, I knowledge the support of the Center of Physics (CF-UM-UP) and the Portuguese Foundation for Science and Technology (FCT) for the research grants funded by the 43/ECUM/CFUM/2022 - QNEEX2D project. I would also like to express my gratitude towards my advisor Doctor Bruno Amorim and my co-advisor Doctor Nuno Peres. Both of them are, without a doubt, one of the best professors and mentors one could have. I thank them for their pedagogical, kind and friendly attitude and for the patience and support shown throughout the realization of this dissertation and also in the courses they taught me. I attribute much of my academic background, specifically in condensed matter physics, to them. Lastly but not least, I’m grateful for my family and girlfriend Sofia who always have shown me their unconditional love. I thank my six month old nephew Tomás for accompanying the writing of the thesis and for useful discussions on exciton physics, i.e sleeping on my bed while I externalize thoughts whilst writing the dissertation. I specially thank my good friends and colleagues Marcelo Barreiro and Tiago Antão for all the great memories we lived together and for all the help and fruitful debates in our academic years. iii Statement of Integrity I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the University of Minho. University of Minho, Braga, august 2023 Francisco Ricardo Lobo Ribeiro iv Assinado por: Francisco Ricardo Lobo Ribeiro Num. de Identificação: 30286406 Data: 2023.08.02 17:31:13+01'00' Abstract In this work we perform a generic derivation on how collective excitations emerge from a many-body system of interacting particles within a time-dependent Hartree-Fock mean-field theory at zero-temperature. To this end, we study the linear response of the system’s reduced density matrix in a many-body perturbation theory and demonstrate that it can be expressed in terms of a generalized eigen-problem of the effective two-particle Hamiltonian of the electron-hole interaction. We then specify this formalism for the case of a crystal system and an atomistic electron-electron interaction, structuring the generalized eigen-problem in terms of the Bloch momentum and spin degrees of freedom. At last, we apply this theory to the case of hexagon boron nitride structures in a nearest-neighbor tight-binding model for the electronic Bloch states. We then solve the generalized eigen-problem numerically and obtain the excitonic states energies and wavefunctions. Also, we comment on the role of screening in the Hartree and Fock interaction, on the numerical details of the generalized eigen-problem and on the reliability of the Tamm-Dancoff approximation. Keywords time-dependent, mean-field approximation, Hartree-Fock, screening, reduced density matrix, zero temperature, linear response, effective two-particle Hamiltonian, generalized eigen-problem, crystal, tight-binding model, exciton, hexagonal boron nitride v Resumo Neste trabalho realizamos uma derivação genérica sobre como excitações colectivas emergem de um sistema de muitos-corpos de partículas interactuantes numa teoria de Hartree-Fock de campo médio dependente do tempo a temperatura zero. Para tal, estudamos a resposta linear da matriz de densidade reduzida do sistema numa teoria de perturbação de muitos-corpos e demonstramos que esta pode ser expressa em termos de um problema generalizado aos valores próprios do Hamiltoniano efetivo de duas partículas da interação eletrão-buraco. Em seguida, especificamos este formalismo para o caso de um sistema cristalino e de uma interação eletrão-eletrão atomística, estruturando o problema generalizado aos valores próprios em termos dos graus de liberdade do momento de Bloch e do spin. Por fim, aplicamos esta teoria ao caso de estruturas de nitreto de boro hexagonal num modelo de tight-binding ao vizinho mais próximo para os estados de Bloch electrónicos. Em seguida, resolvemos numericamente o problema generalizado aos valores próprios e obtemos as energias e as funções de onda dos estados excitónicos. Para além disso, comentamos também no papel da blindagem na interação de Hartree e de Fock, nos pormenores numéricos do problema aos valores próprios e na fiabilidade da aproximação de TammDancoff. Palavras-chave dependente do tempo, aproximação de campo médio, Hartree-Fock, blindagem, matriz de densidade reduzida, temperatura zero, resposta linear, Hamiltonian efetivo de dois-partículas, model de tight-binding, problema generalizado aos valores próprios, cristais excitão, nitreto de boro hexagonal vi Contents 1 Introduction 1 1.1 Exciton basics .................................... 1 1.1.1 Definition and characterization ........................ 1 1.1.2 Excitons in 2D materials ........................... 2 1.1.3 Exciton coupling to light ........................... 2 1.2 Collective excitations in a many-body system ...................... 3 1.3 Thesis Structure ................................... 6 2 Excitonic generalized eigen-problem in a crystal 8 2.1 Introduction to Linear Response Theory ........................ 8 2.2 Time-dependent Hartree-Fock mean-field theory .................... 9 2.3 Linear response theory ................................ 13 2.3.1 Zero temperature regime ........................... 15 3 Excitonic generalized eigen-problem in a crystal 19 3.1 Electron-electron atomistic interaction in the Bloch basis ................ 19 3.1.1 General electron-electron interaction ..................... 21 3.2 Structure of eigen-problem in the Bloch momentum degree of freedom ......... 22 3.3 Structure of eigen-problem in the spin degree of freedom ................ 24 3.4 Screening in the Hartree and Fock terms ....................... 26 4 Numerical Implementation 30 4.1 Discretization of the eigenvalue problem ....................... 30 4.1.1 k-point sampling ............................... 30 4.1.2 Regularization of interaction for small momentum ............... 32 4.1.3 Cutoff of the interaction for large momentum ................. 34 vii malism see [(23)]. In particular, if one is interested in pedagogical discussion of the BSE for a two-particle Green’s function see also [(24)]. 1.3 Thesis Structure We perform a generic derivation on how the collective excitation of a many-body electron system are captured within a time-dependent Hartree-Fock mean-field theory for the special case of insulators at zerotemperature. We then apply this formalism to the case of a hBN isolated monolayer and a hBN-metal hetero-structure. Specifically, in Sec.II, we introduce a general theoretical description of excitons. We start by considering a many-body system of electron described by an Hamiltonian expressed in second quantization containing a single-particle term H0and an interaction term Hint with a generic potential Vαβ γδ of symmetry V(r− r′) = V(r′−r). We then perturb this system with a external applied force such that the equilibrium Hamiltonian is driven out-of-equilibrium. The perturbation Hamiltonian is described as a single-particle term coupled to the external force through some linear operator as usually done in linear response theory. From here, we study the time evolution of the system’s reduced density matrix within a time-dependent Hartree-Fock mean-field scheme Then we expand the reduced density matrix in a power series and retain only first order term in order to study the system’s linear response. This is done in the special case of an insulator at zero temperature such that the single-particle states can be described as being either occupied or empty. Finally, we show that the linear response can be obtained by solving instead a generalized eigenvalue problem for the effective two-particle (electron-hole) Hamiltonian He-h. Secondly, in Sec.III, we derive the specific form of He-h. We start by specifying the focus on crystals systems and derive the explicit form of the potential Vαβ γδ for an atomistic e-e interaction written in a Bloch basis assuming ultra-localized Wannier functions. Using those results we then show the structure of the eigen-problem (still for an insulator at zero temperature) in the Bloch momentum degree of freedom and in the spin degree of freedom, commenting on the role of screening on the Hartree and Fock terms. In Sec.IV we give an example on how to solve the generalized eigen-problem numerically, giving the broad strokes of said numerical implementation while discussing on some important details. Finally, in Sec.V, we solve the generalized eigen-problem to the specific case of an isolated hBN monolayer and obtain its complete excitonic band structure and respective wave-functions. For this, we first show how to obtain the electronic single-particle Bloch states in a nearest-neighbor tight-binding model. We then consider the case of a hBN-metal hetero-structure in order to study the effects of (bulk) metal 6 screening on the excitonic energy levels. 7 Chapter 2 Excitonic generalized eigen-problem in a crystal 2.1 Introduction to Linear Response Theory All local physical measurements of a many-body system amount, in practice, to a process in which a perturbation is created by an applied external force, in the neighborhood of some point r′at some time t′, and then the response of the system is measured at some other point rat some later time t>t′. Consider the case of a many-body system Hwhich is weakly coupled to a one-particle external perturbation Hext through some operator Bsuch that the total system is described by H=H+Hext. In first order perturbation theory, the expectation value of some local observable A(r, t)perturbed by the external force/source F(r, t)can be written as ⟨A(r, t)⟩=⟨A(r, t)⟩0−Zdr′Zdt′χAB (r, t;r′, t′)F(r′, t′),(2.1) where ⟨A(r, t)⟩0is the initial equilibrium expectation value and the “coefficient” of proportionality between the change in the expectation value and the force is the so called generalized susceptibility. On the other hand, the expectation value of A(r, t)in a given initial state |α⟩of H, under the action of the weak perturbation Hext is modified as ⟨A⟩ ≡ ⟨α|A|α⟩ → ⟨α|U−1AU |α⟩where U(t)is the evolution operator in the interaction representation of Htot. By expanding U(t)to first (linear) order in the perturbation Hext, the change in the expectation value yields δ⟨A(r, t)⟩=i ℏZt −∞ dt′⟨[Hext(r, t′), A(r, t)]⟩.(2.2) This expression is known as the generalized Kubo formula . From here we find that the general expression for the generalized susceptibility reads as χAB (r, t;r′, t′) = −i ℏΘ(t−t′)⟨[A(r, t), B(r′, t′)]⟩,(2.3) Fourier −→ χAB (p, ω) = ⟨A(p, ω)⟩ F(p, ω)=”response” ”force” .(2.4) 8 Therefore, if one can express and calculate the commutator [A(r, t), B(r′, t′)], the susceptibility in know and consequently the expectation value of A. However, this is not always such an easy task. In an alternative approach, one can always write a one-particle observable as A=X ab Aabcac† b(2.5) and consequently its expectation value as ⟨A⟩=Pab Aab Dca(t)c† b(t)E. The last term is nothing more that the system’s time-dependent one-particle reduced density matrix, ρba(t) = Dca(t)c† b(t)E,(2.6) which can be obtained by solving its equation of motion in the Heisenberg picture of quantum mechanics. This is exactly the topic of this work. We show how formulate a time-dependent Hartree-Fock theory with the reduced density matrix as the central object. 2.2 Time-dependent Hartree-Fock mean-field theory Consider a many-body system of electrons whose equilibrium state is described by the Hamiltonian H, written in an arbitrary electronic basis {ϕα}, as H=H0+Hint =X αβ hαβc† αcβ+1 2X αβγδ Vαβ γδ c† αc† βcγcδ,(2.7) where the operators c† α(cα) create (annihilate) an electron in the state described by the wavefunction ϕα(r)with rthe electron’s positions vectors, and the greek indices are generic degrees of freedom that the system might have (momentum k, band λ, spin σ, etc...). Also hαβ describes the single-particle matrix elements, hαβ =Zd3rϕ∗ α(r)p2 2m+U(r)ϕβ(r),(2.8) with U(r)a generic static potential (for example a crystal lattice potential) and Vαβ γδ describes the twoparticle interaction matrix elements, Vαβ γδ =Zd3rd3r′ϕ∗ α(r)ϕ∗ β(r′)V(r−r′)ϕγ(r′)ϕδ(r),(2.9) with V(r)a generic electrostatic potential. As usual in electrostatic potentials, we assume that V(r) has the symmetry V(r−r′) = V(r′−r)and consequently the interaction matrix elements obey the symmetries Vαβ γδ =Vβα δγ ,(2.10) Vαβ γδ ∗=Vγδ αβ ,(2.11) 9 as shown in Appendix A.1. Consider now that we add a time-dependent perturbation to the system’s equilibrium Hamiltonian, switched on at t=t0, which is generally described by a one-particle term, Hext =X αβ Bαβc† αcβF(t),(2.12) where Bαβ are matrix elements of a one-particle operator which couples to the time-dependent external (force) field F(t). This term will drive the system out-of-equilibrium by inducing transitions from a given state denoted by the degree of freedom βto another state α. Examples of perturbations can be the charge density coupled to a scalar potential, the current density coupled to a vector potential or the dipole moment coupled to an electric field. The Hamiltonian of the perturbed system reads H(t) = X αβ hαβ +Bi αβFi(t)c† αcβ+1 2X αβγδ Vαβ γδ c† αc† βcγcδ.(2.13) As discussed previously, in order to determine a specific property of a material, we can study its linear response to a given external perturbation by monitoring the time-evolution of the expectation value of a related one-particle observable, A=X ab Aabcac† b,(2.14) which in turn can be written in terms of the system’s time-dependent reduced density matrix (rDM), ρba(t) = Dca(t)c† b(t)E.(2.15) Our work is thus to study the time-evolution of the rDM, specifically in a time-dependent Hartree-Fock scheme. We start by its equation of motions written as d dtρab(t) = i ℏ*dc† b(t) dt ca(t)++i ℏc† b(t)dca(t) dt ,(2.16) where, in the Heisenberg picture of quantum mechanics, the fermionic operators evolve accordingly to the Heisenberg equation d dt ˆ O(t) = i ℏ[H, ˆ O(t)].(2.17) Therefore the rDM equation of motion reads d dtρab(t) = i ℏD[H, c† b(t)]ca(t)E+i ℏDc† b(t)[H, ca(t)]E.(2.18) As shown in Appendix A.2, given the equal-time fermionic anti-commutator properties nca, c† bo=nc† b, cao= δab and {ca, cb}=nc† b, c† ao= 0, each of the commutators are evaluated to [H, c† b] = hαb +Bi αbFi(t)c† α+1 2Vαβ γb c† αc† βcγ−1 2Vαβ bδ c† αc† βcδ,(2.19) [H, ca] = −haβ +Bi aβFi(t)cβ+1 2Vαa γδ c† αcγcδ−1 2Vaβ γδ c† βcγcδ,(2.20) 10 where we have hidden the summations on repeated indices using Einstein’s notation and omitted the time dependency in the fermionic operators for compactness. This practice is done throughout the remaining work with no further mention. Substituting Eqs.(2.19) and (2.20) back into Eq.(2.18) yields −iℏd dtρab(t) = ρaα(t)hαb +Bi αbFi(t)−haβ +Bi aβFi(t)ρβb(t) +1 2Vαβ γb Dc† αc† βcγcaE−1 2Vαβ bδ Dc† αc† βcδcaE +1 2Vαa γδ Dc† bc† αcγcδE−1 2Vaβ γδ Dc† bc† βcγcδE,(2.21) where we came back to the definition of the rDM in Eq.(2.15) for the first two terms. Notice, however, that the remaining interaction terms correspond to the expectation values of four-operators. Dealing with such terms is no easy task and thus, to simplify our model, we introduce a mean-field approximation where we assume the two-particle expectation value to simply behave as a product of two one-particle expectation values. From this mean-field decoupling we can then make use of Wick’s theorem . This is done for the first term of Eq.(2.21) and then we just show the results for the other three. We start by making all the possible two operator contractions from the four operators: Dc† αc† βcγcaE=Dc† αc† βcγcaE+Dc† αc† βcγcaE+Dc† αc† βcγcaE(2.22) with Dc† αc† βcγcaE=Dc† αc† βE⟨cγca⟩,(2.23) Dc† αc† βcγcaE=−c† αcγDc† βcaE,(2.24) Dc† αc† βcγcaE=c† αcaDc† βcγE,(2.25) where the minus sign in the second contraction appears directly from the anti-commutation of the fermionic operators (one could also check for the a minus sign by counting the ntimes the Wick’s contraction line intersects: if nis odd number then there will be a negative sign). Notice that the first term is actually zero since, for non-superconducting systems, we must have that c†c†=⟨cc⟩= 0. Coming back once more to the definition of the rDM in Eq.(2.15), the first four-operator expectation value in Eq.(2.21) within the mean-field approximation read Dc† αc† βcγcaE=−ργα(t)ρaβ(t) + ρaα(t)ργβ(t).(2.26) 11 Following the same treatment, once can easily see that the remaining terms read Dc† αc† βcδcaE=−ρδα(t)ρaβ(t) + ρaα(t)ρδβ(t),(2.27) Dc† bc† αcγcδE=−ργb(t)ρδα(t) + ρδb(t)ργα(t),(2.28) Dc† bc† βcγcδE=−ργb(t)ρδβ(t) + ρδb(t)ργβ(t).(2.29) Putting these terms back together into Eq.(2.21), playing around with mute indices and using the symmetry properties of the potential in Eqs.(2.10) and (2.11), we can rewrite the expression more compactly as iℏd dtρab(t) = haβ +Bi aβFi(t)ρβb(t)−ρaα(t)hαb +Bi αbFi(t) +Vaα δγ −Vaα γδ ρδα(t)ργb(t)−ρaβ(t)Vβα δb −Vαβ δb ρδα(t).(2.30) We now identify the Hartree and Fock self-energy terms respectively as ΣH aγ[ρ(t)] = Vaα δγ ρδα(t),(2.31) ΣF aγ[ρ(t)] = −Waα γδ ρδα(t),(2.32) such that the rDM equation of motion can be written in a time-dependent Hartree-Fock scheme as iℏd dtρab(t) = haβ +Bi aβFi(t)ρβb(t)−ρaα(t)hαb +Bi αbFi(t) +ΣH aγ[ρ(t)] + ΣF aγ[ρ(t)]ργb(t)−ρaβ(t)ΣH βb[ρ(t)] + ΣF βb[ρ(t)].(2.33) This can be written in a more compact and elegant manner as iℏd dtρ(t) = [HHF,ρ(t)] (2.34) with HHF =h+B·F(t) + ΣH[ρ(t)] + ΣF[ρ(t)] (2.35) where HHF is defined as the Hartree-Fock Hamiltonian. Notice that we purposely defined the Fock self-energy term in Eq.(2.32) as an attractive interaction (by inserting the negative sign) such that it lower the total energy of the system enabling the formation of bound electron-hole states. Conversely, the Hartree self-energy term will control details of the excitation spectrum, such as spin-splitting (which will be discussed further in the text). Furthermore, while there is no apparent distinction between this two terms in the time-dependent density matrix Hartree-Fock formalism, in many-body perturbation theory the Hartree and Fock interactions are actually screened in different manners. Thus, forecasting this result, we write Eq.(2.32) instead with an W, symbolizing a screened interaction instead of the bare interaction denoted by V. Further ahead, when discussing two-particle 12 interactions in a many-body perturbation theory, we explain in more detail as to why this is the case. For now, and to put it vaguely, we do not also screen the Hartree interaction (as one would naively expect) because we would be accounting for the screening twice, in the sense that the Hartree interaction already “naturally” screens itself. On another subject, in general, the Hamiltonian halready contains interaction at some kind of meanfield level. Therefore, it is reasonable to subtract this contribution for the self-energy terms, ΣH/F[ρ(t)] →ΣH/F[ρ(t)−ρ(0)].(2.36) This would be equivalent to start from the Hamiltonian with the interacting part subtracted by the equilibrium mean field contribution from the interaction, Vαβ γδ c† αc† βcγcδ→Vαβ γδ c† αc† βcγcδ−(ΣH βγ ρ(0)c† βcγ−ΣF βδ ρ(0)c† βcδ.(2.37) 2.3 Linear response theory The rDM equation of motion as written in Eq.(2.33) would gives us the full response of the system to the external perturbation F(t). However, very frequently, one is only interested in the case where the perturbation is weak and the system responds linearly to it. This is the so-called linear response regime . In linear response theory, as the name indicates, we only account for first (linear) order perturbations. Therefore, we start by expanding the rDM in a power series, ρ(t) = ρ(0) +ρ(1)(t) + ρ(2)(t) + ..., and retain terms only up to the first order, ρ(t)≈ρ(0) +ρ(1)(t).(2.38) Substituting Eq.(2.38) into Eq.(2.33), we see that in the linear regime most terms actually vanish: the terms B·F(t),ρ(1)(t)and ΣH[ρ(1)(t)] + ΣF[ρ(1)(t)],ρ(1)(t)are neglected because they are of second order and the term h,ρ(0)is right away zero since the equilibrium rDM ρ(0) is a function of occupation of the single-particle Hamiltonian h. We are left with iℏd dtρ(1) ab (t)−hacρ(1) cb (t)−ρ(1) ac t)hcb=Bi acFi(t)ρ(0) cb −ρ(0) ac Bi cbFi(t) + ΣH ac[ρ(1)(t)]ρ(0) cb −ρ(0) ac ΣH cb[ρ(1)(t)] + ΣF ac[ρ(1)(t)]ρ(0) cb −ρ(0) ac ΣF cb[ρ(1)(t)].(2.39) 13 Next, we write the time-dependent force and the rDM as their Fourier counterpart F(t) = Zdω 2πe−iωtF(ω),(2.40) ρ(1) ab (t) = Zdω 2πe−iωtρ(1) ab (ω),(2.41) yielding the rDM equation of motion in terms of the frequency of the field ω, ℏωρ(1) ab (ω)−haγρ(1) γb (ω)−ρ(1) aγ (ω)hγb=Bi aγFi(ω)ρ(0) γb −ρ(0) aγ Bi γbFi(ω) + ΣH aγ[ρ(1)(ω)]ρ(0) γb −ρ(0) aγ ΣH γb[ρ(1)(ω)] + ΣF aγ[ρ(1)(ω)]ρ(0) γb −ρ(0) aγ ΣF γb[ρ(1)(ω)].(2.42) Since, at linear order, the system response more aggressively at frequencies ωwhich are in resonance with the characteristic frequencies of the system, we try to rewrite the linear response as ℏωδγδ ab −Hγδ ab ρ(1) ab (ω) = Jab(ω).(2.43) with Hγδ ab the effective two-particle Hamiltonian and Jab(ω)the source term. Specifically, we expect to obtain a free-particles term (where the two particle to not actually interact with each other), a Hartree term and a Fock term. Comparing Eq.(2.42) and Eq.(2.43), we can readily identify the non-interaction term as being haγρ(1) γb (ω)−ρ(1) aγ (ω)hγb =X γδ (haγδδb −δaγhδb)ρ(1) γδ (ω),(2.44) and the source term as being Jab(ω) = Bi aγρ(0) γb −ρ(0) aγ Bi γbFi(ω).(2.45) Although not so straight forward, after some minor play with mute indices while making use of the symmetry properties of the potential in Eqs.(2.10) and (2.11), the two-particle Hamiltonian can be identified as Hγδ ab = (haγδδb −δaγ hδb) + Vδa ϵγ ρ(0) ϵb −ρ(0) aϵ Vϵδ γb −Waδ ϵγ ρ(0) ϵb −ρ(0) aϵ Wϵδ bγ .(2.46) Furthermore, working in the basis that diagonalizes the single-particle Hamiltonian, we have that hac =ϵaδac and ρ(0) ab =faδab,(2.47) where ϵais the occupational energy of the state aand fa=f(ϵa)is the occupation of the state given by the Fermi-Dirac distribution. In this basis, the two-particle Hamiltonian in Eq.(2.46) simply reads Hγδ ab = (ϵa−ϵb)δaγδδb + (fb−fa)Vδa bγ −Waδ bγ .(2.48) Furthermore, the source term in Eq.(2.45) also simplifies to Jab(ω) = (fb−fa)Bi abFi(ω).(2.49) 14 2.3.1 Zero temperature regime Consider the case of an insulator at absolute zero. Since at T= 0Kthe Fermi-Dirac distribution behaves as a Heaviside theta function (centered at the Fermi energy), in the base that diagonalizes the singleparticle Hamiltonian, the degrees of freedom can only be classified as either occupied |o⟩or empty |e⟩ such that fo= 1 or fe= 0. In this regime, the Hamiltonian Hγδ ab has 2×2×2×2 = 16 blocks arising from all the different combination of a→o1or e1,b→o2or e2,γ→o3or e3and δ→o4or e4. Therefore, we can write Eq.(2.43) explicitly as         ℏω1−         He3e4 e1e2He3o4 e1e2Ho3e4 e1e2Ho3o4 e1e2 He3e4 e1o2He3o4 e1o2Ho3e4 e1o2Ho3o4 e1o2 He3e4 o1e2He3o4 o1e2Ho3e4 o1e2Ho3o4 o1e2 He3e4 o1o2He3o4 o1o2Ho3e4 o1o2Ho3o4 o1o2                         ρ(1) e3e4(ω) ρ(1) e3o4(ω) ρ(1) o3e4(ω) ρ(1) o3o4(ω)         =         Je1e2(ω) Je1o2(ω) Jo1e2(ω) Jo1o2(ω)         ,(2.50) where, in order to highlight the absolute zero regime, we introduced the notation J=J(T= 0). It is important to have in in mind that, although we are specifying a degree of freedom, there are other degrees of freedom still hidden away in the indices. That is why we formally referred to the Hamiltonian “entries” as being block matrices (and some an entry). Also, the dimensions of this blocks define the identity matrix 1dimensions which was not yet specified (it is not just a 4×4matrix). In this context, although the occupation function in this regime does really only have the two values, fo= 1 or fe= 0, the same is not true for the occupational energy ϵ, as it depends not only on the empty or occupied definition of the state but also on the other possible degree of freedom. From Eq.(2.48), we see that if a term has a=γ or/and δ=bthen the free-particles term is zero and if it has a=bthen its the interaction term that is zero. In addition, from Eq.(2.49), we see that the source terms with indices only in the occupied or the empty states vanish. We are left with         ℏω1−         He3e4 e1e20 0 0 He3e4 e1o2He3o4 e1o2Ho3e4 e1o2Ho3o4 e1o2 He3e4 o1e2He3o4 o1e2Ho3e4 o1e2Ho3o4 o1e2 0 0 0 Ho3o4 o1o2                         ρ(1) e3e4(ω) ρ(1) e3o4(ω) ρ(1) o3e4(ω) ρ(1) o3o4(ω)         =         0 Je1o2(ω) Jo1e2(ω) 0         .(2.51) Now, the first and last rows of the system of equations in Eq.(2.50) must have the trivial solution ρ(1) e3e4(ω) = ρ(1) o3o4(ω) = 0, and therefore, the linear response system reduces to  ℏω1−  He3o4 e1o2Ho4e3 e1o2 He3o4 o2e1Ho4e3 o2e1     ρ(1) e3o4(ω) ρ(1) o4e3(ω) = Je1o2(ω) Jo2e1(ω) ,(2.52) 15 with rthe generic position of the electron, i.e r=R+sς+rWS with rWS the position within the WignerSeitz unit cell. As a matter of fact, the atomistic interaction approximation used in our model can be obtained directly from the general expression if one expresses the Bloch states in a Wannier basis {wς}, ψkλ(r) = 1 √NX R,ς eik·(R+sς)ϕς kλwς(r−R−sς),(3.11) and assumes ultra-localized orbitals by forcefully introducing the Kronecker delta terms δRR′and δςς′. Substituting Eq.(3.11) directly into Eq.(3.10), within this ultra-localized regime, we obtain ϱkλ k′λ′(q+G)≈X ς e−i(k′−k+(q+G))·sς(ϕς kλ)∗ϕς k′λ′ϱ00 ςς′(q+G)(3.12) with ϱRR′ ςς′(q+G) = ZdDrei(q+G)·rw∗ ς(r−R−sς)wς(r−R′−sς′).(3.13) Setting ϱ00 ςς′(q+G) = 1 yields the same result in Eq.(3.8) of our atomistic model. 3.2 Structure of eigen-problem in the Bloch momentum degree of freedom Consider an electron of momentum kin a fully occupied valence band vwith energy ϵk,v that, by the perturbation of an externally applied electric field, transitions to the empty conduction band cwith momentum k+Q, with Qbeing the center of mass momentum of the exciton, having now energy ϵ{k+Q}c. Note that, in this two-band model, the occupied state corresponds to the valence band index, o→v, and the empty state corresponds to valence band index, e→c. As it turns out, this two-band model is sufficient to describe excitons in hBN, due to the simplicity of its electronic band structure. However, this is not always the case, since in most TMD at least a three band model is necessary [(26)]. The source term that described such perturbation is J{k+Q}c1,kv2(ω, Q+G).(3.14) Accounting for the Bloch momentum degree of freedom, the linear response in Eq.(2.53) would read X k3k4 ℏωS−  Hk3c3,k4v4 {k+Q}c1,kv2Hk3v3,k4c4 {k+Q}c1,kv2 Hk3v3,k4c4 {k+Q}c1,kv2†Hk3c3,k4v4 {k+Q}c1,kv2∗    ρ(1) k3c3,k4v4(ω) ρ(1) k3v3,k4c4(ω) =S J{k+Q}c1,kv2(ω) Jkv1,{k+Q}c2(ω) , (3.15) where the summation over the momentum k3and k4accounts for all the possible different interactions (although, technically, it should have been an integral in the thermodynamic limit). This interaction must, 22 of course, conserve momentum as imposed by the Kronecker deltas in Eq.(3.9) or by the diagrammatics in Fig.(4). For instance, looking at the Hartree term in Fig.(4)(c) we see that the virtual transferred momentum qmust read q=Q+Gwhile for the Fock term in Fig.(4)(b) must read instead q={k−k′}+G, where we rewrote k4→k′. Therefore, the linear response yields X k′ ℏωS−  R(Q)C(Q) C(Q)†R∗(Q)    ρ(1) {k′+Q}c3k′v4(ω) ρ(1) {k′+Q}v3k′c4(ω) =S J{k+Q}c1,kv2(ω, q+G) Jkv1,{k+Q}c2(ω, q+G))  , (3.16) where the resonant block is given by R(Q) = H{k′+Q}c3,k′v4 {k+Q}c1,kv2=ϵ{k+Q}c1−ϵkv2δ{k+Q}c1,{k′+Q}c3δkv2,k′v4+Vk′v4,{k+Q}c1 kv2,{k′+Q}c3−W{k+Q}c1,k′v4 kv2,{k′+Q}c3, (3.17) and the coupling block by C(Q) = H{k′+Q}v3,k′c4 {k+Q}c1,kv2=Vk′c4,{k+Q}c1 kv2,{k′+Q}v3−W{k+Q}c1,k′c4 kv2,{k′+Q}v3.(3.18) We note that formally we should have made explicit composite indices kc1v1;k′c3v4in the resonant and coupling blocks but abstained to do so for the sake of space and tidiness. Finally, using the general expressions in Eqs.(3.9) and (3.8), we can explicitly show the interaction terms in either the resonant and coupling block. The resonant Hartree term yields Vk′v4,{k+Q}c1 kv2,{k′+Q}c3=1 VX G V(Q+G)ϱ{k+Q}c1 kv2 (Q+G)ϱk′v4 {k′+Q}c3 (−(Q+G))(3.19) with, ϱ{k+Q}c1 kv2 (Q+G) = X ς ei(k−{k+Q}+(Q+G))·sςϕς {k+Q}c1∗ϕς kv2(3.20) and, ϱk′v4 {k′+Q}c3 (−(Q+G)) = X ς ei({k′+Q}−k′−(Q+G))·sςϕς k′v4∗ϕς {k′+Q}c3,(3.21) while the resonant Fock term yields W{k+Q}c1,k′v4 kv2,{k′+Q}c3=1 VX G W({k−k′}+G)"ϱ{k+Q}c1 {k′+Q}c3 ({k−k′}+G)#ϱk′v4 kv2 (−({k−k′}+G)) (3.22) with, ϱ{k+Q}c1 {k′+Q}c3 ({k−k′}+G) = X ς ei({k′+Q}−{k+Q}+({k−k′}+G))·sςϕς {k+Q}c1∗ϕς {k′+Q}c3 (3.23) and, ϱk′v4 kv2 (−({k−k′}+G)) = X ς ei(k−k′−({k−k′}+G))·sςϕς k′v4∗ϕς kv2.(3.24) 23 Comparing the resonant and coupling indices in Eqs.(3.17) and (3.18) respectively, we see that the momentum indices stay exactly the same while the band indices swap as v4→c4and c3→v3. Therefore, the expressions for the coupling block are the same as for the resonant block in Eqs.(3.19) through (3.24) but instead with the band indices altered as v4→c4and c3→v3. 3.3 Structure of eigen-problem in the spin degree of freedom Until now we have omitted the spin degrees of freedom, effective dealing with spineless fermions. We will now study the effect of spin in the structure of the effective two-particle Hamiltonian. We consider, however, the spin-orbit interaction to be negligible such that the single-particle states can be simply classified as either spin-up |↑⟩or spin-down |↓⟩states, this way, we just have to introduce the spin degree of freedom σ=↑,↓and sum over all the possible spin combinations of σ3and σ4in Eq.(3.16). The Hilbert space of the electron-hole pairs consists now of eight subspaces: |c3↑3, v4↑4⟩,|c3↑3, v4↓4⟩,|c3↓3, v4↑4⟩, |c3↓3v4↓4⟩and the other four instead with c3→v3and v4→c4. Therefore, the linear response equation (shown within the TDA for the sake of space) yields         ℏω1−         R↑3↑4 ↑1↑2R↑3↓4 ↑1↑2R↓3↑4 ↑1↑2R↓3↓4 ↑1↑2 R↑3↑4 ↑1↓2R↑3↓4 ↑1↓2R↓3↑4 ↑1↓2R↓3↓4 ↑1↓2 R↑3↑4 ↓1↑2R↑3↓4 ↓1↑2R↓3↑4 ↓1↑2R↓3↓4 ↓1↑2 R↑3↑4 ↓1↓2R↑3↓4 ↓1↓2R↓3↑4 ↓1↓2R↓3↓4 ↓1↓2                         ρ(1) c3↑3,v4↑4 ρ(1) c3↑3,v4↓4 ρ(1) c3↓3,v4↑4 ρ(1) c3↓3,v4↓4         =         Jc1↑1,v2↑2 Jc1↑1,v2↓2 Jc1↓1,v2↑2 Jc1↓1,v2↓2         .(3.25) Notice that most of the matrix elements are actually zero due to spin conservation, since the virtual transferred particle does not carry spin i.e the interaction does not induce spin-flips. To see why, we turn once more to a diagrammatic approach and study the diagrams in Fig.(4). For example, in the freeparticles term and in the Fock term represented in Fig.(4)(a) and (b) respectively, we must force the spin σ3to be the same as σ1and the spin σ2to be the same as σ4. One the other hand, in the Hartree term in Fig.(4)(c), we must force the spin σ1to be the same as σ2and the spin σ4to be the same as σ3. Therefore, we are left with R=         ϵ+V−W0 0 V 0ϵ−W0 0 0 0 ϵ−W0 V0 0 ϵ+V−W         ↑↑ ↑↓ ↓↑ ↓↓ ,(3.26) 24 where, for compactness, we wrote ϵ:= ϵ{k+Q}c1−ϵkv2,W:= W{k+Q}c1,k′v4 kv2,{k′+Q}c3and V:= Vk′v4,{k+Q}c1 kv2,{k′+Q}c3. If we now change the basis to U=         1 √2 1 √20 0 0 0 1 0 0 0 0 1 1 √2−1 √20 0         ,(3.27) we effectively decouple the problem into a spin-singlet class of solutions with symmetric subspace (1/√2) (|c↑, v ↑⟩+|c↓, v ↓⟩)for which the Hamiltonian becomes Rs=ϵ−W+ 2V, (3.28) and a spin-triplet class of solutions, consisting of the subspaces |c↑, v ↓⟩ and |c↓, v ↑⟩, and the antisymmetric (1/√2) (|c↑, v ↑⟩−|c↓, v ↓⟩), for which the Hamiltonian becomes Rt=ϵ−W. (3.29) The linear response can thus be solved for singlet and triplet configuration separately, we no further regards to the spin degrees of freedom. This still holds true outside the TDA, considering the coupling block C. Notice that, if the spin-orbit interaction was not negligible, the singlet and triplet configurations would have mixed and the two-particle Hamiltonian must be discussed including its full spin structure. This increases the number of basis states by a factor of 4 and the evaluation of the linear response becomes more difficult. For a discussion on the effect of spin-orbit interaction on the optical spectra (not particularly in hBN but inMoS2) see [(27)]. See that, in full analogy to the electric charge case, an electron of spin σremoved from a fully occupied valence band can be interpreted as the valence band having an overall deficiency of spin σwhich can then be regarded as a hole with spin −σ. In this sense the singlet state can be understood as a zero spin electron-hole bound state. 25 Figure 4: Feynman diagram of the electron-electron interaction including spins for the (a) freeparticles/non-interacting term [corresponding to the first term in Eq.(3.17)], (b) for the resonant Fock interaction [in Eq.(3.19)] and (c) for the resonant Hartree interaction [in Eq.(3.22)]. 3.4 Screening in the Hartree and Fock terms While we prematurely introduced the screened potential W into the Fock self-energy term in Eq.(2.32), we still haven’t discuss the reason why the Fock term is screened while the Hartree term is not. Since we already discussed the many-body perturbation problem and obtained a description for the effective twoparticle collective excitations, we are now more then capable to clearly see why this is the case. Foremost to this discussion, we note that the Hartree term is also commonly referred as being the exchange term and the Fock term as the direct term although, when talking about two electrons interaction outside the context of a excitonic problem, the Fock term is actually the one being called the exchange term and the Hartree as the direct (this is just an unfortunate mismatch). For this discussion, it is useful to take a diagrammatic approach and write the Bethe-Salpeter equation (BSE) as its analogous Dyson equation G=G0+G0ΣGwhere Gis the two-particle propagator, G0is the non-interacting single-particle propagator and Σis the proper/irreducible self-energy of the interaction, i.e that no component of Σcan be written in terms of two self-energy connected by G0. As usual, an 26 iterative solution can be constructed by making the initial ansatz G=G0, obtaining G=G0+G0ΣG0 and then repeatedly substituting the new equality back into Guntil we arrive at G=G0+G0ΣG0+ G0ΣG0ΣG0+.... Diagrammatically, this expansion amounts to an infinite series of diagrams containing all possible combinations of the interaction vertex. Consider, initially, the interaction to be the bare Coulomb V. Setting Σ = ΣH, the Dyson equation for the Hartree interaction is as shown in Fig.(5)(a). Notice that all diagrams from 2rd order onward contain a series of “bubbles” diagrams formed from the combinations VG0V. Analogously, setting Σ=ΣFfor the Fock interaction, we obtain a corresponding series of “ladders” as shown in Fig.(5)(b). Consider now that we substitute our bare Coulomb potential for the random phase approximation (RPA) screened potential WRPA. In the RPA, electrons are assumed to respond to an effective potential WRPA which already accounts for an averaging of screening effects. We can write the renormalized interaction as WRPA =V+VLV+VLVLV+... where Lis the so called Lindhard function corresponding to the proper bubble diagram as shown in Fig.(5)(c). Now see what happens if we substitute the potential WRPA in the place of V in both the Hartree and Fock terms. For the Fock terms it’s easy to see that each “rung” term has now infinite terms, corresponding to each of the bubbles of the interaction line. However, if we try to do the same for the Hartree term, we just obtain diagrams that are already accounted elsewhere in the original expansion. This multiple counting of equivalent terms is strictly forbidden since the interaction potential must be proper, so it is concluded that the interaction V in the Hartree term must not be screened, at least in the RPA sense. Said in other words, the bubble series in Fig.(5)(b) already “naturally” screens itself. We note that some studies argue that, when solving the BSE in a restricted subspace of the full Hilbert space (for example, the subspace associated with low-energy bands close to the Fermi energy), the Hartree interaction should be appropriately screened by states outside of said subspace (for example, higher energy bands and/or other physical subsystems such as substrates). This is the so called S-approximation [(28)]. As we will explain further ahead, for the case of isolated hBN, due to the simplicity of its band structure, we actually do not need to account for the screening of higher order bands. However, when dealing with the hBN-metal hetero-structure, we do screen both the Hartree and Fock term due to the proximity of the metal. For a more in dept discussion on the screening in the exchange term and the S-approximation see the article [(28), (29)]. As a compendium see also the many-body theory books [(30), (31)]. We take this discussion as an opportunity to point out that we are working within the static limit, W(ω= 0), meaning that we do not account for the dynamics (i.e the frequency dependency) of the screened potential W(ω). As a side note, in a non-equilibrium Dyson equation sense, this would be 27 equivalent to close the system of equations for the lesserand greater-GF, also know as Kadanoff-Baym equation, by setting the collision integrals to zero [(32)]. As discussed in the article [(33)], taking into account the dynamical effects is possible however, instead of obtaining a simple eigenvalue problem, one obtains a non-linear one, this is, we would need to solve the BSE self-consistently because the Fock selfenergy itself would depend on the resulting excitonic energy. Fortunately, as also pointed out in their work, taking the static limit is a reasonable approximation for most semiconductor crystals since the plasmon energies that control the dynamic of the screening are much bigger than the excitonic binding energies, effectively closing the iterative process. Figure 5: Diagrammatic approach to the iterative Dyson equation by expansion of the electron-hole propagator Gin powers of the (a) Hartree self-energy (b) Fock self-energy. (c) RPA expansion of the effective potential WRPA. On the topic of screening, we take this opportunity to lay out the explicit form of the static potentials describing the e-e interaction. We comment on the results in three dimensions but also in two dimensions as a precursor to the study of hBN. Firstly, the bare interaction of the Hartree/exchange term corresponds to the Coulomb potential, V(r) = e2 4πε0 1 |r|,(3.30) expressed in reciprocal space. Performing the 3D and 2D Fourier transforms, one obtains, respectively, V3D(k) = e2 ε0 1 |k|2and V2D(k) = e2 2ε0 1 |k|.(3.31) While for the case of a 3D dielectric the screened interaction of the Fock/direct can be obtained just by directly making the alteration ε0→εin the 3D bare potentials, the same is not true for the 2D case. As 28 derived in Appendix B.1, the suitable choice of potential to account for the repulsive screened interaction between electrons in a polarizable 2D semiconductor is the Rytova-Keldysh potential, W(k) = e2 2ε0 1 |k| 1 1 + r0|k|,(3.32) where r0=χ2D/2 is the so called effective screening radius which is material dependent. The term 1/(1 + r0k)account for the RPA dielectric screening as discussed above. 29 Chapter 4 Numerical Implementation In this section we discuss on the numerical implementation details of the generalized eigen-problem in Eq.(2.64). For this discussion, we refer to bright excitons in hBN as a concrete example and visual aid, however, we will refrain ourselves to comment beyond numerical details since this section is prior to that of excitons hBN. In due time such comments will be made and refer to this section so one gets a complete grasp of some details. We start by showing how the discretization of the eigenvalue problem is made with a k-point sampling of the electronic first Brillouin zone, tackling the complications that emerge such as the small apparent momenta divergence of the interaction and the need to implement a cutoff for large momentum. Next, in order to obtain the excitonic energies and wave-functions, we review possible numerical eigen-solvers outside and within the Tamm-Dancoff approximation. Finally, we test the convergence of the excitonic energies with the number of k-point samplings and norm of cutoff as a compromise between numerical precision and computational cost (on a 8core machine). 4.1 Discretization of the eigenvalue problem 4.1.1 k-point sampling Looking back at Eq.(3.16), although in the thermodynamic limit the sum over the electron’s momentum k′should be actually an integral, in order for us to be able to withdraw any information about the excitonic band structure and wave-function, we need to retain k′as a sum over a discretized grid of Nkequidistant points within the 1BZ separated by some interval ∆k, Z1BZ dk (2π)2→∆kX k∈1BZ .(4.1) 30 This can be done in a mid-point rule as kn1,n2(b1,b2) = 2n1−√Nk−1 2√Nk×b1+2n2−√Nk−1 2√Nk×b2,with n1, n2∈1,2, ..., pNk, (4.2) where b1and b2are reciprocal space vectors and the interval length reads ∆k=|b1|/√Nkalong both the ˆ b1and ˆ b2direction. In Fig.(6) we shown a concrete example of a discretized 1BZ grid for the case of hBN. For example, substituting n1=√Nkand n2= 1/2 √Nk+ 1into Eq.(4.2), we obtain the point near the edge of the 1BZ along the ˆ b1direction just shy of the symmetry point M(which is exactly on the edge of the 1BZ) by half a division length. Notice that, to guarantee that n2is a whole number we must set Nkas an odd perfect square. Also, with Nkbeing odd, it is guaranteed that the Γsymmetry point is exactly hit. In addition, particularly for the hBN case, to ensure that the immensely relevant valley symmetry points ±Kare exactly part of the grid, i.e ±K! =kn1,n2(bhBN 1,bhBN 2), one needs to specifically use NK∈ {3+6N}2where Nare natural numbers. Figure 6: Discrete hBN first Brillouin zone grid with Nk= 441 points (corresponding to √Nk= 21 points along each b1/b2direction) defined in a mid-point rule in the electronic k-space. After performing this k-point sampling, the linear response equation in Eq.(3.16) (shown only within the TDA for the sake of space) reads      ℏω1Nk×Nk−     R{k′ 1+Q}c3,k′ 1v4 {k1+Q}c1,k1v2R{k′ 2+Q}c3,k′ 2v4 {k1+Q}c1,k1v2 k′ → R{k′ 1+Q}c3,k′ 1v4 {k2+Q}c1,k2v2R{k′ 2+Q}c3,k′ 2v4 {k2+Q}c1,k2v2 k′ → ↓k↓k ...                ρ(1) {k′ 1+Q}c3k′ 1v4 ρ(1) {k′ 2+Q}c3k′ 2v4 ↓k′      =     J{k1+Q}c1,k1v2 J{k2+Q}c1,k2v2 ↓k      , (4.3) where kand k′run on all the points of the discretized 1BZ. From here we can construct our two-particle Hamiltonian matrix with each entry of the resonant block Rand coupling block Ccalculated using the expressions in Eqs.(3.19) through (3.24) for a fixed exciton center of mass momentum Q. The electronic 31 Figure 10: Convergence of the first eight energy levels for a bright exciton in hBN as a function of (1st column) the number of k-sampling points Nkusing |kcutoff|= 1.5Å−1and ( 2nd column ) as a function of the cutoff norm |kcutoff|using Nk≈9000 points. Each energy convergence plot is accompanied by the absolute error in a log10 scale and the relative error between energies calculated with consecutive increments of the respective parameter. There are only 5states visible because we treated degenerated states as being the same state (which we checked to behaved the same apart from some very minor variations.). The blue shaded area corresponds to the values taken as the optimal compromise between numerical precision and computational cost. The pink shaded area characterizes the low-energy regime and has radius |kcutoff|=|M−K|/2 ≈0.43Å−1. 38 Chapter 5 Excitons on hBN structures 5.1 Tight-binding model for the single-particle Bloch states Hexagonal boron nitride (hBN) is a 2D material composed of a simple layer of alternating boron and nitrogen atoms disposed in a planar honeycomb lattice, as shown in Fig.(11)(a). hBN shares a lot of similarities with graphene, also a 2D honeycomb structured material but instead composed of only carbon atoms. The most relevant distinction is that graphene behaves as a semi-metal with a zero-gap at its Dirac points while hBN, due the different electrostatic environment in the boron and in the nitrogen atom, has an opening gap of about ϵg= 5.9eV (there are actually a lot of different results for ϵgin the literature however the mentioned value is one of the more commonly reported [(37)]. Also, hBN has a slightly larger lattice constant than graphene (about 1.8%), being around a0= 2.5Å [(38)]. The planar honeycomb lattice can be described as a triangular Bravais lattice generated by the real vectors basis a1=a0 21,√3,(5.1) a2=a0 2−1,√3.(5.2) In each Wigner-Seitz cell, we have one atom of boron and one atom of nitride, which we designate as sub-lattices Aand Brespectively, having positions, sA=(0,0),(5.3) sB=a0 √3(0,1).(5.4) For each site A, the position of the nearest-neighbors (NN) in the sites Bare given by δ1=a0 √3(0,1),(5.5) δ2=a0 2√3−√3,−1,(5.6) δ3=a0 2√3√3,−1.(5.7) 39 All these vectors are shown in Fig.(11)(a) within the real space lattice. Furthermore, from the real lattice basis vectors follow the reciprocal lattice basis vectors b1=2π √3a0√3,1,(5.8) b2=2π √3a0−√3,1.(5.9) which are shown in Fig.(11)(b) together with the first zone of Brillouin, which they form. Figure 11: (a) hBN real space honeycomb lattice constructed from two superposed triangular sub-lattices of boron atoms (depicted in red), denoted as sub-lattice A, and of nitrogen atoms (depicted in blue), denoted as sub-lattice B. The vectors a1and a2are the lattice basis vectors and δ1,δ2and δ3are the nearest-neighbor vectors. (b) hBN reciprocal space lattice with b1and b2its basis vectors. The first Brillouin zone is emphasized in light gray while the remaining are only outlined. The red dots correspond to the Dirac points ±Kand the green dots correspond to the 1BZ edges Mpoints. In order to compute the excitonic energies using Eqs.(3.19) through (3.24), we first need to know the single-particle Bloch wave-functions ϕς kλ. In this work we calculate them in a nearest-neighbors tightbinding model. For this, we write the system’s single-particle NN tight-binding Hamiltonian in real space as HTB(R) = X i ϵAa† RiaRi+X i ϵBb† RibRi−tX ⟨i,j⟩a† RibRi+δj+b† RjaRi−δj,(5.10) where the operators a† Ri(aRi) create (annihilate) an electron in the sub-lattice Ain a given Bravais lattice site Riwhile the operators b† Ri(bRi) create (annihilate) an electron instead in the sub-lattice B(in a 40 given Bravais lattice site Ri). Therefore, the first two terms correspond to the isolated single-particles Hamiltonian of the site Aand B, respectively, and the last term to the hybridization between neighboring sites iand j, describing the possible hoppings from site Ato site Band vice-versa. We only assess hopping terms up to the first neighbors terms, which is denote by ⟨i, j⟩, and consider a static hopping term in either direction, i.e tRi,Rj→ −t. Notice that, contrarily to graphene, since the atoms on sites A and Bare different the single-particle energies ϵAand ϵBare inherently different. We can represent the NN TB Hamiltonian in Eq.(5.10) in reciprocal space by expressing the creation/annihilation operators as their Fourier counterparts, aRi=1 √VX k e+ik·(Ri+sA)ak,(5.11) bRi=1 √VX k e+ik·(Ri+sB)bk,(5.12) and rearranging the expression just that the identity δ(k−k′) = 1/NPie−iRi·(k−k′)is apparent. We obtain HTB(R) = X k ϵAa† kak+X k ϵBb† kbk−tX kγka† kbk+γ† kb† kak,(5.13) with the newly-defined γcomplex number, γk=X ⟨j⟩ e+ik·δj.(5.14) If we now define a row vector c† k=ha† kb† kiwe can rewrite the system’s Hamiltonian as HTB R= Pkc† kHTB kckwith the hBN NN TB Hamiltonian matrix being HTB(k) =   ϵA−tγk −tγ† kϵB .(5.15) Within this simplified tight-binding model, the expression for the electronic two-band structure can easily be obtained analytically by diagonalizing the matrix in Eq.(5.15), yielding E± TB(k) = ±v u u tϵ2+t2"3+2cos (a0kx)+4cos a0√3 2ky!cos a0 2kx#.(5.16) Here we defined the zero point energy at (ϵA+ϵB)/2 and defined ϵ≡(ϵA−ϵB)/2 at the middle of the gap such that ϵA=ϵand ϵB=−ϵ. The valence band corresponds to the E− TB(k)dispersion while the E+ TB(k)corresponds to the conduction band, as shown in Fig.(12) which is accompanied by the density of states DoS(E) = Pkδ(E−E(k)) . Notice that, if ϵA=ϵB, as is the case for 41 graphene, we obtain ϵ= 0 and the band dispersion closes in a linear fashion at the so called Dirac points, K±= (±4π/(3a0),0). In hBN, the electronic band dispersion is also at its minimum near these points but has instead a parabolic shape. In either case, this points represent a fundamental symmetry of the system, called valley parity. To see why the dispersion is parabolic at these valley points, we Taylor series expand the exponential of γkin Eq.(5.14) near k→K+pwith p→0. We obtain e+ip·δj≈1 + ip·δj. Now, since P⟨j⟩e+iK·δj= 0 we are left with γK+p≃ip·X ⟨j⟩ e+iK·δjδj=−√3a0 2(px−ipy).(5.17) Invoking the Pauli matrices definitions, from Eq.(5.15) we can write the TB Hamiltonian Hk TB in this lowenergy regime as HTB(K+p) = ϵσz+t√3a0 2(p·σ),(5.18) which clearly resembles the 2D Dirac Hamiltonian, HDirac =σzmc2+c(p·σ)with ϵtaking the role of the rest mass energy mc2and instead with a velocity vF=t√3a0/2 , termed the Fermi velocity , as a replacement to the velocity of light c. Notice that, for the case of graphene, since ϵ= 0, the electrons would behave as if they are massless. In this limit, the hBN low-energy dispersion can be written as the typical relativistic dispersion relation ETB(K+p) = ±qp2v2 F+m2 effv4 F.(5.19) where meff is the effective mass of the electron at a given point near the valleys. Although not captured in this simple tight-binding model, if one does some type of DFT to obtain a more complete electronic band structure, one could see that the hBN bands do actually cross between themselves (see, for example, Fig.(1) from [(39)]). This appears to be troublesome to our two-band timedependent Hartree-Fock mean-field theory since certain transitions could occur between bands that are not accounted for in our model. However, these other intersecting bands corresponds to electronic states that are orthogonal to the ones we use in our two-band model and thus, will not interfere (i.e, even if we accounted for this other bands in our model, the form factors in Eqs.(3.19)-(3.24) would always give zero for transitions between those bands). However, this only applies for the hBN case since, if one was dealing instead with TMDs, one would need to account for at least three bands [(26)]. Furthermore, following the context of TMDs, one should also need to account for the spin-orbit coupling (SOC) where the effective Hamiltonian for such a system could be obtained by adding to Eq.(5.18) the term HSOC =tstsz(σz−1)/2 where tsquantifies the spin-orbit coupling and sz=±1labels the spin projection of the bands [(26), (40)]. 42 Figure 12: hBN electronic band structure from a nearest-neighbor tight-binding model accompanied by the density of the states. The dispersion goes along the symmetry path k: Γ →K→M→Γ and was calculated using ϵg= 7.8eV for the energy gap, t= 3.1eV for the hopping parameter and a0= 1.42√3Å for the honeycomb lattice length. 5.2 Isolated hBN excitonic properties 5.2.1 Bright exciton: singlet state for Q= 0 Firstly, we focus on the results, in and out of the TDA, for the (optical active) bright excitons, corresponding to the singlet state for a zero center of mass momentum exciton. Bright excitons, as opposed to dark excitons, dominate the optical properties of semiconductors since they can form/recombine from a single photon absorption/emission, making them the main focus of most of the literature on excitons. The hBN NN TB electronic parameters to perform the calculation were: e2/ε0= 104/55.3eVÅ, ϵg= 7.8eV for the energy band-gap, t= 3.1eV for the hopping parameter, a0= 1.42√3Å for the honeycomb lattice length and r0= 10Å for the effective screening radius of the Rytova-Keldysh potential. We have selected these values because they correspond to those featured in the reference [(41)], with which we intend to make a comparative analysis of the results. Furthermore, following the convergence tests done in Sec.IV, we will be using a k-sampling of Nk= 8649 points with a cutoff norm of |kcutoff|= 1.5Å. We show the results for the bright exciton energies in Fig.(13) and the results for the corresponding wave-function intensities in both reciprocal and real space in Fig.(15). The real space representation was 43 obtained via ΨX cv(R+sα,R′+sβ) = 1 NkX k eik·((R+sα)−(R′−sβ)) ϕβ kv∗ϕα kcΨX cv(k),(5.20) where R+sαis the position of the electron in the real space lattice and R′+sβthe position of the hole. In particular, the results shown in Fig.(15) have the hole fixed at the center boron atom in the origin of the referential, i.e R′=0and sβ=sA. Such expression can be obtained by considering the rDM in real space, ρRα,R′β=Dc† R′,βcR,αEand inverse-Fourier transforming the creation/annihilation as in Eq.(3.2). Analyzing the energy state results in Fig.(13), we observe a 1st and a 2nd lowest energy state having both the same energy 6.31eV. This states are then separated from the next higher energy state by 0.75eV. Subsequently, we obtained a 4th state considerably above the 3rd by 0.09eV and then a succession of two pairs of states, 5th and 6th, and 7th and 8th, having also the same energy (in pairs). We see that the lowest energy state is isolated from the higher energy states by a gap of 0.75eV, compared to the 0.25eV range between the 3st and 8nd states. Furthermore, the emergence of states having the same energies is to be expected due to valley parity symmetry where, and speaking from a low-energy scheme in the vicinity of K±, each decoupled valley would contributes with its own eigenvalue. This, of course, assumes that inter-valley coupling in nearly absent which is not necessarily true as seen in the degeneracy lift between the 3rd and 4th state. Formally, we must sum the total wave-function intensities of the states that share the same energy, in order to preserve the natural symmetry of the (perfect) crystal lattice [(42)]. We do this in both the reciprocal and real intensities plots in Fig.(15) and denote the now single double-degenerated states by their pure constituents states in a square bracket notation. If we were to break the symmetry of the crystal, for example by slightly displacing one of the atoms, we would expect a splitting of these constituent states. In order to understand if our calculations are trustworthy, we compare our results with those obtained in previous studies, specifically the references [(39), (41), (43), (44)]. We note that these comparisons are not meant to be one-to-one because the electronic parameters and methodology do not exactly match between the different articles and ours. Also, we note that in [(41)] refers not to the excitonic energies EXbut instead to their binding energy Ebwith respect to the electronic energy gap, i.e Eb=EX−ϵg. Moreover, the order of appearance of the states may depend on whether the electronic calculations where based on a TB or an ab inition model. In particular, comparing with the results from [(39), (43)], we see that in our TB calculations the non-degenerated 3rd and 4th states, which are degenerate in their ab inition calculations, are switched with the degenerated 5th and 6th states, which are non-degenerate in their ab inition calculations. Beside this points, we consider the values and behavior of our results within reason 44 with the mentioned studies: a lowest energy state around ∼6eV isolated from the higher energy states by a gap of ∼1eV with three pairs of degenerated states and two non-degenerated states separated by ∼0.1eV. Figure 13: hBN bright exciton energies (corresponding to the optically active singlet state with center of mass momentum Q= 0) for the first eight excitonic states. The energy values read as: 6.31eV for the 1st and 2nd state, 7.06eV for the 3rd, 7.15eV for the 4th, 7.18eV for the 5th and 6th and 7.31eV for the 7th and 8th. The parameters utilized for these results are made explicit in the beginning of Sec.V.B.I. Also, the calculation were performed within the TDA. We now refer to the results for the wave-functions shown in Fig.(15). Since in the Wannier limit the exciton can be treated as a hydrogenoide model [(1)], we can classify and identify to some extent the excitonic states in an scheme borrowed from the 2D atomic orbitals terminology: (n, ℓ, m)with nthe principal quantum number, ℓ= 0, ..., n −1the azimuthal quantum number and m=−ℓ, ...ℓ the magnetic quantum number, denoting the states in the format 1s≡(1,0,0),2s≡(2,0,0),2p0≡ (2,1,0),2p1≡(2,1,1), etc... The probability densities of these hydrogenic atomic orbitals are shown in Fig.(14) as a visual reference guide. 45 Figure 14: 2D probability density projection onto the plane zOyof the hydrogenic atomic orbitals in a (n, ℓ, m)representation. Image altered from [(45)] Evaluating side by side the probability densities of the atomic orbitals in Fig.(14) with our results for the hBN bright excitonic in Fig.(15), one could assess through a visual comparison that the excitonic lower energy state resembles the 1satomic orbital and that the higher excited state resembles the 2satomic orbital. Although the comparison for the in-between states is not visually evident, possibly due to trigonal warping and/or hybridizing of the sand pbehavior, we can at least assess that, since they are strikingly similar, they must belong to the same orbital family. Given they are only three different states on that orbital family, #3, #4 and #[5,6], it must correspond to the 2porbital family. Thus, we assess the #[5,6] as being the 2p0state and the non-degenerated #3 and #4 states as being the 2p−1and 2p+1 states respectively. Whilst not necessarily common, the 2slevel is indeed found to be above the 2plevels. 46 Figure 15: hBN bright exciton (normalized) wave-function intensities |ψX|2for the first eight excitonic states. The representation in done in (1st row) reciprocal space and in (2nd row) real space. The single double-degenerated states are constructed by summing the intensities of the states that share the same energy, in order to preserve the natural symmetry of the crystal lattice. In the real space representation, the hole (denoted by a small black dot) is fixed at the center boron atom in the origin of the referential. While this hydrogenoide classifications of the excitonic states is useful, it its the genuine triangular point group symmetry C3vthat should formally describe the excitonic states. The C3vsymmetry group is described by three different classes: the identity E, two 3-fold rotation symmetries C3and three mirror symmetries σv, reflecting along the axis of highest rotational symmetry. It decomposes into three irreducible representations (irreps): the single-degenerated symmetric irrep A1, the single-degenerated anti - symmetric irrep A2, characterized by an odd character for the σvreflections, and the double-degenerated irrep E[(46)]. Immediately, we can identify the double degenerate 2s,2pand 2p0as being irreducible presented by E. It then only remains to make the correspondence between the 2p+1 and 2p−1states and the irreps A1or A2. For this, we do an intensity-phase representation of the excitonic wave-function where each k-point of the discretized 1BZ has an associated color value given by to the wave-function phase and an opacity proportional to the normalized intensity, as shown in Fig.(16). From here, and ignoring the numerical noise, it is clear that the 3rd state, corresponding to the 2p−1state, is symmetric with respect to a highest rotational symmetry axis, in particular the ky-axis, and therefore is irreducible represented by the A1irrep. One the other hand, the 4th state, corresponding to the 2p+1 state, is anti-symmetric with respect to the ky-axis and is instead irreducible represented by the A2irrep. However, contrary to our findings, in [(39)], which also studied the excitonic states in terms of the C3vsymmetric, found the anti-symmetric state to be instead the first of the non-degenerate states. 47 Figure 21: Metal screened potentials in a log10 scale as a function of the distance dbetween the hBN monolayer and the metal The bare potential, corresponding to d→ ∞, is colored black. The utilized effective screening radius reads r0= 10Å. In Fig.(22), we repeat the calculations for the first eight energy levels of the bright exciton, (corresponding to the singlet state with Q= 0) but instead with the metal screened potentials as a function of the distance dbetween the hBN monolayer and the bulk metal in Eqs.(5.21) and (5.22). As one would expect, the greater the metal screening, the lower the excitonic binding energies, Eb=ϵg−EX. In a crude explanation, this happens because the electric field lines that go outside the hBN and pierce the metal hamper the electron-hole interaction making them less attracted to each other. In the limiting case where the metal is so far away that its screening has no effect on the hBN monolayer, the excitonic energies tend to the isolated monolayer values. Figure 22: hBN-metal hetero-structure bright exciton energies as a function of the distance dbetween the hBN monolayer and the bulk metal, for the first eight excitonic states. The distance dincrement is 5Å. The parameters utilized for these results are made explicit in the beginning of Sec.V.B.I. Also, the calculation were performed within the TDA. 54 Chapter 6 Conclusions and future work In this work, we showed how collective excitations can emerge from a many-body system of interacting particle in a time-dependent reduced density matrix Hartree-Fock mean-field theory and then we applied this theory to the particular case of hBN structures in order to obtain its excitonic states. In Sec.II, we studied the linear response of a system of electrons from a many-body perturbation theory by inspecting the time evolution of the system’s rDM in a Hartree-Fock scheme. To this end, and as a mean to simplify the rDM’s equation of motion, two main approximations were made. Firstly, we introduced a mean-field approximation where we decoupled two-particle expectation values into a product of two one-particle expectation values. Secondly, we imposed the case of an insulator at zero temperature such that the occupation degree of freedom could only be classified as either occupied or empty. From these significant simplifications, it was shown that the first order rDM could ultimately be obtained by solving instead a generalized eigen-problem for the effective two-particle Hamiltonian of the electron-hole interaction. Furthermore, in Sec.III, we particularized this formalism to the case of a crystal system and an atomistic e-e interaction written in a Bloch basis assuming ultra-localized Wannier functions. We then derived the explicit structure of the eigen-problem in terms of the Bloch momentum and spin degree of freedom. Concerning the spin structure, we showed that (and neglecting spin-orbit effects) the eigen-problem decouples into a singlet and triplet set of solutions, which differ by the contribution of the repulsive Hartree/exchange interaction. In addition, we discussed on the role of screening on the Hartree and Fock terms. We clarified why the Hartree interaction does not need to be screened, at least in the RPA sense, and commented on why taking the screening static limit is justifiable. Subsequently, in Sec.IV we tackled some details of a possible numerical solution for the generalized eigen-problem. Finally, in Sec.V, we focused on the particular case of hBN structures and described the electronic single-particle Bloch states in a nearest-neighbor tight-binding model. We then solved the generalized eigen-problem and obtained its excitonic band structure and respective wave-functions for the first eight 55 excitonic states. Foremost, we concerned ourselves with bright excitons, corresponding to the optically active singlet state for near-zero excitonic center of mass momentum. We observed the emergence of pairs of states having the same energies, which is to be expected due to valley parity symmetry, and stated that we must consider these pair of states together in order to preserve the natural symmetry of the crystal. Thereby, for calculations within the TDA, we obtained an energy of 6.32eV for the lowest energy state #[1,2], which is isolated from the next higher energy state by a gap of 0.75eV. Subsequently, we obtained a 4th state considerably above the 3th by 0.09eV. We then classified and identified these states by their wave-functions in an scheme borrowed from the 2D atomic orbitals and through the triangular point group symmetry C3vof the real space lattice. Furthermore, concerning to the whole excitonic band dispersion, we commented on the overall structure and degeneracies of the singlet states, then on the spin-splitting between the singlet and triplet state, and finally on the validity of the Tamm-Dancoff approximation. We found that, while the spin-splitting for zero excitonic center of mass momentum was not verified in our case, for larger momentum the spin-splitting is clearly visible. Finally, concerning the reliability of the TDA, we found no difference whatsoever on the results calculated in and out of the TDA, thus corroborating its appeal. By the end of this section, we considered the case of a hBN-metal hetero-structure in order to study the effects of (bulk) metal screening on the excitonic energy levels. We found that, as expected, as the metal gets closer to the hBN monolayer and the screening gets more intense, the lower the excitonic binding energies are. As a continuation of this work, the next step would be to take the obtained excitonic energies and eigenvalues and calculate some optical properties following the discussion of Sec.II.C.b. Moreover, we could calculate the expected life-time of the excitons, for example, through decay into the electromagnetic modes of a planar laser cavity. As future lines of research, it would be interesting to apply the formalism of this work to other 2D structures, such as TMDs in a three-band model. Indeed, even within the context of hBN, there still remains numerous ideas to explore. How does the excitonic energies vary with the effective screening radius r0? Or with a dielectric environment other than the vacuum? What are the effects on the excitonic band structure due to spin-orbit coupling? What about other structures such as bilayers, or twisted bilayers, or periodic alternating metal-dieletric-hBN structures? 56 Bibliography [1] “Knox, r. s. (1983). introduction to exciton physics. in collective excitations in solids (pp. 183-245). boston, ma: Springer us.,” [2] “Dvorak, m., wei, s. h., wu, z. (2013). origin of the variation of exciton binding energy in semiconductors. physical review letters, 110(1), 016402.,” [3] “Egri, i. (1979). a simple model for the unified treatment of wannier and frenkel excitons. journal of physics c: Solid state physics, 12(10), 1843.,” [4] “N. rytova, moscow univ. phys. bull. 3, 30 (1967),” [5] “L.v. keldysh, sov. j. exp. theor. phys. lett. 29, 658 (1979),” [6] “Cudazzo, p., sponza, l., giorgetti, c., reining, l., sottile, f., gatti, m. (2016). exciton band structure in two-dimensional materials. physical review letters, 116(6), 066803.,” [7] “Mueller, t., malic, e. (2018). exciton physics and device application of two-dimensional transition metal dichalcogenide semiconductors. npj 2d materials and applications, 2(1), 29.,” [8] “Massicotte, m., vialla, f., schmidt, p., lundeberg, m. b., latini, s., haastrup, s., ... koppens, f. h. (2018). dissociation of two-dimensional excitons in monolayer wse2. nature communications, 9(1), 1633.,” [9] “Wang, g., chernikov, a., glazov, m. m., heinz, t. f., marie, x., amand, t., amp; urbaszek, b. (2018). colloquium: Excitons in atomically thin transition metal dichalcogenides. reviews of modern physics, 90(2).,” [10] “Quintela, m. f., peres, n. m. (2020). a colloquium on the variational method applied to excitons in 2d materials. the european physical journal b, 93, 1-16.,” [11] “Kalugin, n. g., rostovtsev, y. v. (2009). ” dark” and” bright” excitons in carbon nanotubes: New media for quantum optics. journal of nanoelectronics and optoelectronics, 4(3), 302-306.,” 57 [12] “Berkelbach, t. c., hybertsen, m. s., reichman, d. r. (2015). bright and dark singlet excitons via linear and two-photon spectroscopy in monolayer transition-metal dichalcogenides. physical review b, 92(8), 085413.,” [13] “Ye, z., cao, t., o’brien, k., zhu, h., yin, x., wang, y., ... zhang, x. (2014). probing excitonic dark states in single-layer tungsten disulphide. nature, 513(7517), 214-218.,” [14] “Loh, k. p. (2017). brightening the dark excitons. nature nanotechnology, 12(9), 837-838.,” [15] “Zhou, y., scuri, g., wild, d. s., high, a. a., dibos, a., jauregui, l. a., ... park, h. (2017). probing dark excitons in atomically thin semiconductors via near-field coupling to surface plasmon polaritons. nature nanotechnology, 12(9), 856-860.,” [16] “Zhang, s., li, b., chen, x., ruta, f. l., shao, y., sternbach, a. j., ... basov, d. n. (2022). nanospectroscopy of excitons in atomically thin transition metal dichalcogenides. nature communications, 13(1), 542.,” [17] “Kim, y., kim, j. (2021). near-field optical imaging and spectroscopy of 2d-tmds. nanophotonics, 10(13), 3397-3415,” [18] “Saiki, t., matsuda, k., nomura, s., mihara, m., aoyagi, y., nair, s., takagahara, t. (2004). nano-optical probing of exciton wave-functions confined in a gaas quantum dot. microscopy, 53(2), 193-201.,” [19] “Zhang, x. x., cao, t., lu, z., lin, y. c., zhang, f., wang, y., ... heinz, t. f. (2017). magnetic brightening and control of dark excitons in monolayer wse2. nature nanotechnology, 12(9), 883-888.,” [20] “Gelly, r. j., renaud, d., liao, x., pingault, b., bogdanovic, s., scuri, g., ... lončar, m. (2022). probing dark exciton navigation through a local strain landscape in a wse2 monolayer. nature communications, 13(1), 232.,” [21] “Jiang, j., pachter, r. (2022). analysis of localized excitons in strained monolayer wse 2 by first principles calculations. nanoscale, 14(31), 11378-11387.,” [22] “Darlington, t. p., carmesin, c., florian, m., yanev, e., ajayi, o., ardelean, j., ... schuck, p. j. (2020). imaging strain-localized excitons in nanoscale bubbles of monolayer wse2 at room temperature. nature nanotechnology, 15(10), 854-860.,” [23] “Blase, x., duchemin, i., jacquemin, d., loos, p. f. (2020). the bethe–salpeter equation formalism: From physics to chemistry. the journal of physical chemistry letters, 11(17), 7371-7382.,” 58 [24] “Strinati, g. (1988). application of the green’s functions method to the study of the optical properties of semiconductors. la rivista del nuovo cimento (1978-1999), 11(12), 1-86.,” [25] “Salpeter, e. e. (2008). bethe-salpeter equation–the origins,” [26] “Liu, gui-bin shan, wen-yu yao, yugui yao, wang xiao, di. (2013). three-band tight-binding model for monolayers of group-vib transition metal dichalcogenides. physical review b. 88. 10.1103/physrevb.88.085433.,” [27] “Molina-sánchez, a., sangalli, d., hummer, k., marini, a., wirtz, l. (2013). effect of spin-orbit interaction on the optical spectra of single-layer, double-layer, and bulk mos 2. physical review b, 88(4), 045412.,” [28] “Benedict, l. x. (2002). screening in the exchange term of the electron-hole interaction of the bethesalpeter equation. physical review b, 66(19).,” [29] “Qiu, d. y., da jornada, f. h., louie, s. g. (2021). solving the bethe-salpeter equation on a subspace: Approximations and consequences for low-dimensional materials. physical review b, 103(4), 045117.,” [30] “Cottam, m. g., haghshenasfard, z. (2020). many-body theory of condensed matter systems: An introductory course. cambridge university press.,” [31] “Jishi, r. a. (2013). feynman diagram techniques in condensed matter physics. cambridge university press.,” [32] “Lipavský, p., Špička, v., velický, b. (1986). generalized kadanoff-baym ansatz for deriving quantum transport equations. physical review b, 34(10), 6933.,” [33] “Rohlfing, michael louie, steven. (2000). electron-hole excitations and optical spectra from first principles. phys. rev. b. 62. 10.1103/physrevb.62.4927,” [34] “Andrilli, s., hecker, d. (2022). elementary linear algebra. academic press.,” [35] “Arnoldi, w. e. (1951). the principle of minimized iterations in the solution of the matrix eigenvalue problem. quarterly of applied mathematics, 9(1), 17-29.,” [36] “Gazzola, s., nagy, j. g. (2014). generalized arnoldi–tikhonov method for sparse reconstruction. siam journal on scientific computing, 36(2), b225-b247.,” 59 [37] “Cassabois, g., valvin, p., gil, b. (2016). hexagonal boron nitride is an indirect bandgap semiconductor. nature photonics, 10(4), 262-266.,” [38] “Ishigami, m., aloni, s., zettl, a. (2003, december). properties of boron nitride nanotubes. in aip conference proceedings (vol. 696, no. 1, pp. 94-99). american institute of physics.,” [39] “F. ferreira, a. j. chaves, n. m. r. peres, and r. m. ribeiro, ”excitons in hexagonal boron nitride singlelayer: a new platform for polaritonics in the ultraviolet,” j. opt. soc. am. b 36, 674-683 (2019),” [40] “Scharf, b., xu, g., matos-abiague, a., Žutić, i. (2017). magnetic proximity effects in transition-metal dichalcogenides: converting excitons. physical review letters, 119(12), 127403.,” [41] “Quintela, mauricio henriques, j. peres, nuno. (2022). theoretical methods for excitonic physics in two-dimensional materials: A tutorial.,” [42] “Wirtz, l., marini, a., grüning, m., attaccalite, c., kresse, g., rubio, a. (2008). comment on “huge excitonic effects in layered hexagonal boron nitride”. physical review letters, 100(18), 189701.,” [43] “Galvani, thomas paleari, fulvio miranda, henrique molina-sánchez, alejandro wirtz, ludger latil, sylvain amara, hakim ducastelle, françois. (2016). excitons in boron nitride single layer. physical review b. 94. 10.1103/physrevb.94.125303.,” [44] “Wu, f., qu, f., macdonald, a. h. (2015). exciton band structure of monolayer mos 2. physical review b, 91(7), 075310.,” [45] “Qijing zheng, matplotlib: Hydrogen wave function,” [46] “Dresselhaus, m. s., dresselhaus, g., jorio, a. (2007). group theory: application to the physics of condensed matter. springer science business media.,” [47] “Flórez, f. g., siebbeles, l. d., stoof, h. t. c. (2020). effects of material thickness and surrounding dielectric medium on coulomb interactions and two-dimensional excitons. physical review b, 102(12), 125303.,” 60 Appendices 61 Appendix A Details on the theoretical description of excitons A.1 Symmetry properties of the interaction matrix elements We show the details on how to arrive at the symmetry properties presented in Eqs.(2.10) and (2.11) of the interaction matrix elements, as described in Eq.(2.9), of the many-body system of electron Hamiltonian in Eq.(2.7). From the imposed symmetry V(r−r′) = V(r′−r), we obtain Vαβ γδ =Zd3rd3r′ϕ∗ α(r)ϕ∗ β(r′)V(r−r′)ϕγ(r′)ϕδ(r) =Zd3rd3r′ϕ∗ α(r′)ϕ∗ β(r)V(r′−r)ϕγ(r)ϕδ(r′) =Zd3rd3r′ϕ∗ α(r′)ϕ∗ β(r)V(r−r′)ϕγ(r)ϕδ(r′) =Zd3rd3r′ϕ∗ β(r)ϕ∗ α(r′)V(r−r′)ϕδ(r′)ϕγ(r) =Vβa δγ .(A.1) Similarly Vαb γδ ∗=Zd3rd3r′ϕα(r)ϕβ(r′)V(r−r′)ϕ∗ γ(r′)ϕ∗ δ(r) =Zd3rd3r′ϕ∗ δ(r)ϕ∗ γ(r′)V(r−r′)ϕδ(r′)ϕα(r) =Vδγ βα =Vγδ αβ (A.2) A.2 Commutators and Anti-commutator properties We show the details on how to arrive at the commutator presented in Eqs.(2.19) and Eqs.(2.20). We evaluate each commutator individually given the fermionic and commutator operator properties, nca, c† bo=nc† b, cao=δab (A.3) {ca, cb}=nc† b, c† ao= 0 (A.4) 62 Foremost, we derive commutator and anti-commutator properties for various numbers of operator in a general fashion and then apply then to the explicit case at hand. The commutator and anti-commutator operate, respectively as, [A, Z] = AZ −ZA (A.5) {A, Z}=AZ +ZA. (A.6) From here it follows that {A, Z}={Z, A}(A.7) [A, Z] = −[Z, A](A.8) Consequently we have that [AB, Z] = ABZ −ZAB =ABZ −ZAB +AZB −AZB =A(BZ +ZB)−(ZA +AZ)B =A{B, Z}−{Z, A}B(A.9) and [AB, Z] = ABZ −ZAB =ABZ −ZAB −AZB +AZB =A(BZ −ZB)−(ZA −AZ)B =A[B, Z]−[Z, A]B(A.10) Using both of the properties above we further obtain [ABCD, Z] = [(AB)(CD), Z] = (AB)[(CD), Z]−[E, (AB)](CD) =AB (C{D, Z}−{Z, C}D)+(A{B, Z}−{Z, A}B)CD (A.11) 63 of length η→0“around” each respective interface. Firstly, integrating from −d−ηto −d+ηyields −∂z˜ V−d+η −d−η=−χm∂z˜ V(z=−d−η)−−χd∂z˜ V(z=−d+η) ⇒−Bk||e−k||d−C(−k||)e+k||d−Ak||e−k||d=−χmAk||e−k||d+χd(Bk||e−k||d+C(−k||)e+k||d) ⇔−B+Ce2k||d+A=−χmA+χd(B−Ce2k||d) ⇔A(1 + χm)−(B−Ce2k||d)(1 + χd) = 0 (B.25) Similarly, integrating instead from −ηto +ηyields −∂z˜ V +η −η=−e ε0 +χ2D(ik||)2˜ V(z= 0) −χd∂z˜ V(z=−η)−−χd∂z˜ V(z= +η) ⇒−(D−E)k|| −(B−C)k||=−e ε0−χ2Dk||2(B+C)−χd(B−C)k|| +χd(D−E)k|| ⇔−D+E+B−C=−e ε0 1 k|| −χ2Dk|| (B+C)−χd(B−C) + χd(D−E) ⇔B(1 + χ2Dk|| +χd) + C(−1 + χ2Dk|| −χd)−(D−E)(1 + χd) = −e ε0 1 k|| (B.26) And lastly, integrating from h−ηto h+ηyields −∂z˜ V h+η h−η=−χd∂z˜ V(z=h−η) ⇒−F(−k||)e−k||h−Dk||ek||h+E(−k||)e−k||h=−χdDk||ek||h+E(−k||)e−k||h ⇔F+De2k||h−E=−χdDe2k||h+χdE ⇔F+ (De2k||h−E)(1 + χd) = 0 (B.27) Thus, to find the coefficients Athrough Fone needs to solve the system of equations A(1 + χm)−(B−Ce2k||d)(1 + χd) = 0 (B.28) B(1 + χ2Dk|| +χd) + C(−1 + χ2Dk|| −χd)−(D−E)(1 + χd) = −e ε0 1 k|| (B.29) F+ (De2k||h−E)(1 + χd) = 0,(B.30) together with the tangential continuity of the electric field at the interfaces, Ae−k||d−(Be−k||d+Cek||d) = 0,(B.31) (B+C)−(D+E) = 0,(B.32) (Dek||h+Ee−k||h)−Fe−k||h= 0.(B.33) Directly from the system of equations we have, Fe−k||h=−(1 + χd)Dek||h−Ee−k||h(B.34) Fe−k||h=Dek||h+Eek||h(B.35) 70 meaning that D=χd 2 + χd e−2k||hE≡Gd(k||)E(B.36) where we defined a new function Gd(k||). Also directly from the system of equations, we have (1 + χm)Ae−k||d= (1 + χd)Be−k||d−Cek||d(B.37) Ae−k||d=Be−k||d+Cek||d(B.38) and thus C=−χm−χd 1 + χd+χm e−2k||dB≡Gm(k||)B(B.39) where, once again, we defined another new functionGm(k||)related to Gd(k||)via 1 + Gm(k||)B=1 + Gd(k||)E(B.40) We obtain 1 + χ2Dk|| +χdB+−1 + χ2Dk|| −χ1Gm(k||)B+ (−1−χd)Gd(k||)E+ (1 + χd)E=−1 k|| e ε0 ⇒1 + χ2Dk|| +χ1+−1 + χ2Dk|| −χdGm(k||)+(1+χd)1−Gd(k||)1 + Gm(k||) 1 + Gd(k||)B=−1 k|| e ε0 (B.41) resulting in the potential ˜ V=B+C=−e k||ε0εd"χ2Dk|| εd +2−Gd(k||)Gm(k||) 1 + Gd(k||)1 + Gm(k||)#−1 (B.42) Now, in the special case of a 3D perfect metal, χm→ −∞, we have that Gm(k||) = −χm−χd 1 + χd+χm e−2k||d≈ −χm χm e−2k||d=−e−2k||d(B.43) and in the special case of the dielectric being actually a vacuum, χd= 0 we have that Gd(k||) = χd 2 + χd e−2k||h=−e−2k||h(B.44) In this cases, the effective potential, corresponding to the screened Rytova-Keldysh potential yields V(k||) = e 2ε0 1 k|| 1 r0k|| +1 2 ek||d sinh(k||d) (B.45) 71 I knowledge the support of the Center of Physics (CF-UM-UP) and the Portuguese Foundation for Science and Technology (FCT) for the research grants funded by the 43/ECUM/CFUM/2022 - QNEEX2D project.