# Fileset

[ComputCondMatt41(2024)e00977.pdf](https://mdr.nims.go.jp/filesets/4cd8ae4b-88f2-4b24-9874-bdbd60fa8c75/download)

## Creator

Kaoru Ohno, [Ryoji Sahara](https://orcid.org/0000-0003-0788-2985), Takeshi Nanri, Yoshiyuki Kawazoe

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Optical properties of rutile TiO<math display="inline">  <msub>    <mrow></mrow>    <mrow>      <mn>2</mn>    </mrow>  </msub></math> with Zr, Mo, Zn, Cd impurities](https://mdr.nims.go.jp/datasets/2c3cf691-86a0-4bda-94c6-a4d293e81b22)

## Fulltext

Optical properties of rutile TiO[formula omitted] with Zr, Mo, Zn, Cd impuritiesComputational Condensed Matter 41 (2024) e00977 A2Contents lists available at ScienceDirectComputational Condensed Matterjournal homepage: www.elsevier.com/locate/cocomFull length articleOptical properties of rutile TiO2 with Zr, Mo, Zn, Cd impuritiesKaoru Ohno a,b,∗, Ryoji Sahara b, Takeshi Nanri c, Yoshiyuki Kawazoe d,e,fa Department of Physics, Graduate School of Engineering, Yokohama National University, 79-5 Tokiwadai, Hodogaya-ku, Yokohama 240-8501, Japanb Research Center for Structural Materials, National Institute for Materials Science (NIMS), 1-2-1 Sengen, Tsukuba, Ibaraki, 305-0047, Japanc Research Institute for Information Technology, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japand New Industry Creation Hatchery Center, Tohoku University, 6-6-4 Aramaki Aza Aoba, Aoba-ku, Sendai, Miyagi, 980-8579, Japane School of Physics, Institute of Science, Suranaree University of Technology, 111 University Avenue, Nakhon Ratchasima, 30000, Thailandf Physics and Nanotechnology, SRM Institute of Science and Technology, Kattankurathur, Tamil Nadu, 603203, IndiaA R T I C L E I N F OKeywords:TitaniaFirst principlesPhotoabsorption spectraQuasiparticleMany-body perturbation theoryk-point samplingA B S T R A C TTo explore quasiparticle (QP) energy gaps and photoabsorption spectra of rutile TiO2 with nonmagnetictransition metal (Zr, Mo, Zn, Cd) impurities, we conducted a 𝛤 -point only 𝐺𝑊 + Bethe–Salpeter equation(BSE) calculation on a 72 (or 71) atom supercell. Our findings reveal that Zn and Cd impurities must coexist,at least partly, with oxygen vacancies to maintain charge neutrality. Among the systems considered, Mo,Zn, or Cd doped rutile TiO2 may exhibit optical absorption and catalytic activity under visible light. Theresulting QP energy gaps (𝛥𝜀QP) and photoabsorption energies (PAEs) are fairly in good agreement with bothexperimental and theoretical data currently available. The necessary conditions for the applicability of the𝛤 -point only approach in the 𝐺𝑊 + BSE framework were found to be: (1) The 𝛤 -point only 𝐺𝑊 calculationshould reproduce a reasonable band gap. (2) The ‘‘superficial’’ exciton binding energy (the diagonal elementof 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 between 𝑣 = VBM and 𝑐 = CBM, where 𝑊 and 𝑋 are the direct and exchange terms of theBSE matrix elements, respectively) must be positive or marginally negative. (3) The ‘‘real’’ exciton bindingenergy (𝛥𝜀QP− the lowest PAE) should be positive, even if it is exceptionally small.1. IntroductionRutile and anatase are two important polymorphs of TiO2. Anataseirreversibly reverts to rutile at high temperatures, while rutile is struc-turally more stable than anatase at elevated temperatures, exhibit-ing greater strength and stability. In 1972, Fujishima and Honda [1]demonstrated that n-type rutile TiO2 can facilitate water splitting underillumination, claiming that the material does not suffer degradation.This report sparked significant interest in the semiconductor-electrolytecommunity [2]. The photocatalytic activity of TiO2 [3,4] is applicableto water [5] and environmental [6] purification, while dye-synthesizedTiO2 exhibits photovoltaic behavior relevant to solar cells [7]. More-over, TiO2 shows antibacterial activity [8], biocompatibility [9], andso on. Rutile is also used in optical materials, ultraviolet (UV)-cutcosmetics, pigments, and colorant, and rutile quartz is used in jewelry.In these applications, the optical properties of the transition metaldoped rutile TiO2 are a central focus, as pure rutile TiO2 has a wideband gap of 3.3 eV and requires UV radiation to excite electrons fromthe valence band maximum (VBM) to the conduction band minimum(CBM).∗ Corresponding author at: Department of Physics, Graduate School of Engineering, Yokohama National University, 79-5 Tokiwadai, Hodogaya-ku, Yokohama240-8501, Japan.E-mail address: ohno@ynu.ac.jp (K. Ohno).So far, numerous experimental and theoretical investigations havefocused on rutile TiO2 with transition metal impurities. Magnetic tran-sition metal impurities such as Cr, Mn, Fe, Co, Ni, and Ru exhibitinteresting behaviors, including properties of dilute magnetic semi-conductors [10–13] and enhanced conductivity [14,15], in additionto their optical characteristics [16–18]. In contrast, the optical prop-erties have been the main subject of nonmagnetic transition metalimpurities, such as Zr, Mo, Zn, and Cd. A range of studies, encom-passing experimental and theoretical (density functional theory (DFT))approaches, have been carried out for Zr [19–22] and Mo [23–26].Additionally, several experimental studies focused on Zn [27,28] andCd [28–30]. DFT calculations have been performed for Mo and Cd bySong et al. [31], for Zn by Saini et al. [32], and for Cd by Sato et al. [33]and Errico et al. [34]. However, all these calculations were based onDFT, which is not suitable for studying excited states.The Green’s function technique in many-body perturbation the-ory [35,36] provides a robust first-principles approach for determiningthe quasiparticle (QP) energies of materials, which can be preciselymeasured experimentally using photoelectron spectroscopy. In Hedin’shttps://doi.org/10.1016/j.cocom.2024.e00977Received 27 August 2024; Received in revised form 28 October 2024; Accepted 28vailable online 6 November 2024 352-2143/© 2024 The Authors. Published by Elsevier B.V. This is an open access a October 2024rticle under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). https://www.elsevier.com/locate/cocomhttps://www.elsevier.com/locate/cocommailto:ohno@ynu.ac.jphttps://doi.org/10.1016/j.cocom.2024.e00977https://doi.org/10.1016/j.cocom.2024.e00977http://creativecommons.org/licenses/by/4.0/K. Ohno et al.𝑃nhp𝐺ag(fiksacsiedor(swcwrttcpesebComputational Condensed Matter 41 (2024) e00977 set of equations [35], the exchange–correlation component of the self-energy is represented as 𝛴xc = 𝑖𝐺 𝑊 𝛤 , where 𝐺 is the one-particleGreen’s function and 𝛤 is the vertex function. The dynamically screenedCoulomb interaction is expressed as 𝑊 = 𝜀−1𝑣, where 𝑣 is the bareCoulomb interaction and 𝜀 = 1 − 𝑣𝑃 is the dielectric function, with= −𝑖𝐺 𝐺 𝛤 serving as the polarization function. Here, it is important toote that the normalization problem of the QP wave functions [37,38]as recently been solved by Nakashima et al. [39]. Furthermore, thetheory is not limited solely to the ground state; it is also applicableto excited eigenstates as the initial reference state [40,41]. The sim-lest approximation is to assume 𝛤 = 1, which is referred to as the 𝑊 approximation [35,42]. The well-known one-shot 𝐺 𝑊 (𝐺0𝑊0)pproach [43–48], which uses the Kohn–Sham (KS) wave functionsand eigenvalues [49] from DFT [50], reproduces a reasonable energyap. Accurate analysis of excitonic effects and photoabsorption spectraPAS) necessitates solving the Bethe–Salpeter equation (BSE) [36,51]ollowing the 𝐺0𝑊0 calculation [52–57]. This comprehensive procedureis recognized as the 𝐺 𝑊 + BSE approach.When addressing the BSE for periodic crystals, k-point samplingis generally employed [54–57]. However, in systems with defects ormpurities, the supercell size inevitably becomes substantial, making-point sampling cumbersome and time consuming. For such largeupercell systems, it is reasonable to justify the use of the 𝐺 𝑊 + BSEpproach without k-point sampling. Thus, a 𝛤 -point only 𝐺 𝑊 + BSEalculation [52,53,57–60] would be applicable to extensive supercellystems.In this paper, we apply the 𝛤 -point only 𝐺 𝑊 + BSE approach tonvestigate QP energy gaps and PAS of rutile TiO2 with nonmagnetictransition metal impurities. The main objective is to assess the feasi-bility of a 𝛤 -point only 𝐺 𝑊 + BSE calculation for crystalline rutileTiO2, both with and without impurity atoms, using a 72 (or 71) atomsupercell. We also identify a typo in the key expression of the electron–hole interaction 𝑊𝑣𝑐;𝑣′𝑐′ in the BSE in the pioneering papers [51,57]. So,we present a detailed derivation of the correct equation in Appendix.The impurities under consideration include Zr, Mo, Zn, and Cd, withach impurity atom substituting one of the Ti atoms in the supercell.The impurity concentration is 𝑥 = 1∕24 = 4.2 at% in the chemical formTi1−𝑥M𝑥O2 (where M = Zr, Mo, Zn, Cd). If the impurity concentrationincreases further, the PAS may change drastically, although we havenot performed such calculations. Throughout the course of our calcula-tions, we observed that an oxygen vacancy (OV) must accompany theZn and Cd impurity atoms. This is due to the instability of Zn2+ andCd2+ when placed in the Ti4+ substitutional site of TiO2 without thepresence of an OV.2. MethodAfter a careful convergence check (see Section 4) for pure rutileTiO2, we use a 2 × 2 × 3 tetragonal supercell (derived from the six-atom tetragonal primitive cell) in the present study. This supercell hasimensions 𝑎 = 𝑏 = 9.1453 Å and 𝑐 = 8.7767 Å, including a totalf 72 atoms. For the impurity systems, one of the titanium atom iseplaced with a Zr, Mo, Zn, or Cd atom. For simplicity, spin polarizationi.e., local magnetic moment) is not considered. Additionally, for theystems with a Zn or Cd impurity atom, we also examine configurationsith an OV near the impurity atom.We use Vienna ab initio simulation package (VASP) [61] to optimizethe entire supercell and the atomic positions within it, employing thelocal density approximation (LDA) [49,62] with a plane wave (PW)utoff energy of 500 eV. For the calculation of optical properties,e adopt the all-electron mixed basis approach [58–60,63–68] usingTOMBO program [69], where the one-particle wave functions areexpressed as a linear combination of both PWs and atomic orbitals(AOs). All core and (truncated) valence numerical AOs (confined insidethe nonoverlapping atomic spheres) are utilized alongside the PWs to 𝑊2 describe the electronic states. The eigenvalue problem is solved self-consistently within the LDA [49,62] using approximately 11,000 PWscorresponding to a cutoff energy of 25.9 Ry.In the one-shot 𝐺 𝑊 approach [43–48], the QP wave functionsare approximated by the Kohn–Sham orbitals, denoted as 𝜙𝑐 (𝒓) forconduction (𝑐) and 𝜙𝑣(𝒓) for core–valence (𝑣) levels. We can separatethe 𝐺 𝑊 self-energy 𝛴𝑥𝑐 (𝜔) (excluding the Hartree term) into the Fockexchange term 𝛴𝑥 and the correlation term 𝛴𝑐 (𝜔), such that 𝛴𝑥𝑐 (𝜔) =𝛴𝑥+𝛴𝑐 (𝜔). Their diagonal matrix elements are expressed, respectively,as⟨𝜙𝑛 |𝛴𝑥 |𝜙𝑛 ⟩ = −occ∑𝑣 ∫ 𝜙∗𝑛(𝒓)𝜙𝑣(𝒓)1|𝒓 − 𝒓′|𝜙∗𝑣(𝒓′)𝜙𝑛(𝒓′) 𝑑𝒓𝑑𝒓′ (1)and⟨𝜙𝑛 |𝛴𝑐 (𝜔) |𝜙𝑛 ⟩ = ∫ 𝑑𝒓𝑑𝒓′𝜙∗𝑛(𝒓)𝜙𝑛(𝒓′) 𝑖2𝜋 ∫∞−∞𝑑 𝜔′𝑒−𝑖𝜔′0+×[emp∑𝑐𝜙𝑐 (𝒓)𝜙∗𝑐 (𝒓′)𝜔 − 𝜔′ − 𝜀LDA𝑐 + 𝑖0++occ∑𝑣𝜙𝑣(𝒓)𝜙∗𝑣(𝒓′)𝜔 − 𝜔′ − 𝜀LDA𝑣 − 𝑖0+]×[𝑊 (𝒓, 𝒓′;𝜔′) − 1|𝒓 − 𝒓′|], (2)where 𝜀LDA𝑐 and 𝜀LDA𝑣 are the LDA energy eigenvalues, and 𝑊 (𝒓, 𝒓′;𝜔) isthe dynamically screened Coulomb interaction. In calculating the corre-lation term and the polarization function 𝑃 = −𝑖𝐺 𝐺, the core contribu-tion is ignored in the occupied (occ) summation, while a large numberof empty states are included in the empty (emp) summation: 531 forCd@rutile, 534 for Mo@rutile, 535 for Zr@rutile and Cd+OV@rutile,540 for Zn@rutile, 544 for pure rutile and Zn+OV@rutile, totaling1000 levels including core and valence levels. Hereafter, we refer toutile TiO2 doped with impurity A as ‘‘A@rutile.’’ The full 𝜔′ integra-ion [67] is performed at the 𝐺 𝑊 level, as it has been demonstratedthat this integration is particularly necessary for TiO2 [70]. Then,he QP energies are evaluated using (𝜇LDA𝑥𝑐 being the LDA exchange–orrelation potential) [44]𝜀QP𝑛 = 𝜀LDA𝑛 +𝑍𝑛 ⟨𝜙𝑛 |𝛴𝑥𝑐 (𝜀LDA𝑛 ) − 𝜇LDA𝑥𝑐 |𝜙𝑛 ⟩ (3)with a renormalization factor𝑍𝑛 =[1 − 𝜕⟨𝜙𝑛 |𝛴𝑥𝑐 (𝜔) |𝜙𝑛 ⟩𝜕 𝜔|||||𝜔=𝜀LDA𝑛]−1. (4)In the BSE [51,57,58] within the Tamm–Dancoff approximation [51,57], the singlet excitation energies 𝛺 (s)𝑟 are obtained by solving thematrix eigenvalue equation∑𝑣′ ,𝑐′(𝐷𝑣𝑐;𝑣′𝑐′ + 2𝑋𝑣𝑐;𝑣′𝑐′ )𝐴(s)𝑟 (𝑣′, 𝑐′) = 𝛺 (s)𝑟 𝐴(s)𝑟 (𝑣, 𝑐), (5)while the triplet excitation energies 𝛺 (t )𝑟 are obtained from∑𝑣′ ,𝑐′𝐷𝑣𝑐;𝑣′𝑐′ 𝐴(t )𝑟 (𝑣′, 𝑐′) = 𝛺 (t )𝑟 𝐴(t )𝑟 (𝑣, 𝑐). (6)Here, 𝐴(s)𝑟 (𝑣, 𝑐) and 𝐴(t )𝑟 (𝑣, 𝑐) denote the singlet and triplet BSE am-litudes, respectively, representing the expansion coefficients of thelectron–hole (exciton) wave functions (see Appendix):𝜒𝑟(𝒓, 𝒓′) = −∑𝑣,𝑐𝐴𝑟(𝑣, 𝑐)𝜙𝑐 (𝒓)𝜙∗𝑣(𝒓′). (7)The matrix 𝐷 represents a diagonal term minus a direct term, given by𝐷𝑣𝑐;𝑣′𝑐′ = (𝜀QP𝑐 − 𝜀QP𝑣 )𝛿𝑣𝑣′𝛿𝑐 𝑐′ −𝑊𝑣𝑐;𝑣′𝑐′ , (8)where 𝑊𝑣𝑐;𝑣′𝑐′ denote the matrix elements of the direct term, repre-enting the screened electron–hole interaction, while 𝑋𝑣𝑐;𝑣′𝑐′ in Eq. (5)denote the matrix elements of the exchange term. Therefore, in singletexcitation, the elements 𝑊𝑣𝑐;𝑣′𝑐′− 2𝑋𝑣𝑐;𝑣′𝑐′ act as the exciton binding en-rgy. Since the most important matrix element is the diagonal elementetween 𝑣 = VBM and 𝑐 = CBM, we will refer to the corresponding− 2𝑋 as the ‘‘superficial’’ exciton binding energy. The matrix𝑣𝑐;𝑣𝑐 𝑣𝑐;𝑣𝑐K. Ohno et al.CtftitiqamTsu𝛥3smopsrcbOComputational Condensed Matter 41 (2024) e00977 elements 𝑊𝑣𝑐;𝑣′𝑐′ are the frequency integral of the dynamically screenedoulomb interaction 𝑊 (𝒓, 𝒓′;𝜔) multiplied by two energy denominatorerms (see Appendix for the detailed derivation and a typo in thepioneering papers [71]),𝑊𝑣𝑐;𝑣′𝑐′ =∫ 𝑑𝒓𝑑𝒓′𝜙∗𝑐 (𝒓)𝜙𝑐′ (𝒓)𝜙𝑣(𝒓′)𝜙∗𝑣′ (𝒓′) 𝑖2𝜋 ∫∞−∞𝑑 𝜔𝑊 (𝒓, 𝒓′;𝜔)𝑒−𝑖𝜔0+×⎡⎢⎢⎣1𝛺𝑟 − 𝜔 − (𝜀QP𝑐′ − 𝜀QP𝑣 ) + 𝑖0++ 1𝛺𝑟 + 𝜔 − (𝜀QP𝑐 − 𝜀QP𝑣′ ) + 𝑖0+⎤⎥⎥⎦,(9)where 𝛺𝑟 denotes the photoabsorption energy (PAE) obtained as ainal result. By decomposing 𝑊 (𝒓, 𝒓′;𝜔) as 1∕|𝒓 − 𝒓′| + [𝑊 (𝒓, 𝒓′;𝜔) −1∕|𝒓 − 𝒓′| ], Eq. (9) can be separated into two terms as𝑊𝑣𝑐;𝑣′𝑐′ = 𝑈𝑣𝑐;𝑣′𝑐′ + 𝑆𝑣𝑐;𝑣′𝑐′ , (10)where the matrix elements 𝑈𝑣𝑐;𝑣′𝑐′ represent the bare Coulomb term,𝑈𝑣𝑐;𝑣′𝑐′ = ∫ 𝜙∗𝑐 (𝒓)𝜙𝑐′ (𝒓)1|𝒓 − 𝒓′|𝜙𝑣(𝒓′)𝜙∗𝑣′ (𝒓′) 𝑑𝒓𝑑𝒓′, (11)and the matrix elements 𝑆𝑣𝑐;𝑣′𝑐′ represent the remaining correlationerm. In evaluating 𝑆𝑣𝑐;𝑣′𝑐′ , the generalize plasmon pole (GPP) model[44] is used to avoid the 𝜔 integration; see Appendix. On the otherhand, the matrix elements 𝑋𝑣𝑐;𝑣′𝑐′ of the exchange term are given by𝑋𝑣𝑐;𝑣′𝑐′ = ∫ 𝜙∗𝑐 (𝒓)𝜙𝑣(𝒓)1|𝒓 − 𝒓′|𝜙𝑐′ (𝒓′)𝜙∗𝑣′ (𝒓′) 𝑑𝒓𝑑𝒓′. (12)We evaluate 𝑈𝑣𝑐;𝑣′𝑐′ , 𝑆𝑣𝑐;𝑣′𝑐′ , and 𝑋𝑣𝑐;𝑣′𝑐′ as well as ⟨𝜙𝑛 |𝛴𝑥 |𝜙𝑛 ⟩ and⟨𝜙𝑛 |𝛴𝑐 (𝜔) |𝜙𝑛 ⟩ in Eqs. (1) and (2) in Fourier space as (A.32)–(A.34)n Appendix; for ⟨𝜙𝑛 |𝛴𝑥 |𝜙𝑛 ⟩ and ⟨𝜙𝑛 |𝛴𝑐 (𝜔) |𝜙𝑛 ⟩, see Ref. [44]. How-ever, for the evaluation of 𝑈𝑣𝑐;𝑣′𝑐′ and 𝑋𝑣𝑐;𝑣′𝑐′ as well as ⟨𝜙𝑛 |𝛴𝑥 |𝜙𝑛 ⟩,we exclude the fully on-site AO-only contributions, which are sep-arately calculated in a single integral along the radial direction, asaccurately as possible, using analytic integral forms of the Poissonequation for one-center problem. Finally, to obtain PAS, contributionsfrom three mutually orthogonal polarization directions — (100), (010)and (001) — are averaged using the formula provided in Ref. [54]. Wehave utilized all these formulas in our previous studies [58–60].The exchange terms ⟨𝜙𝑛 |𝛴𝑥 |𝜙𝑛 ⟩ and 𝑋𝑣𝑐;𝑣′𝑐′ and the bare Coulomberm 𝑈𝑣𝑐;𝑣′𝑐′ are evaluated using approximately 88,000 𝑮 vectors, cor-responding to a cutoff energy of 103.6 Ry. In contrast, the polarizationfunction 𝑃 = −𝑖𝐺 𝐺 and the correlation terms ⟨𝜙𝑛 |𝛴𝑐 (𝜔) |𝜙𝑛 ⟩ and𝑆𝑣𝑐;𝑣′𝑐′ are evaluated with approximately 4000 𝑮,𝑮′ vectors, corre-sponding to a cutoff energy of 13.2 Ry. In solving the BSE, we include50 empty levels along with all occupied valence levels. We confirmedthat these parameters are sufficient to achieve convergence in theresulting QP and PAEs within 0.1 eV.Lastly, we briefly comment on the tuning of the code within thehybrid MPI-OpenMP parallel architecture. The 𝐺 𝑊 + BSE portion ofTOMBO requires extensive computations, making performance opti-mization crucial. To reduce computational time, several techniqueswere implemented in this study. First, the conditional branches at thennermost loop nesting structure was removed and the computation se-uence was adjusted to ensure continuous memory access. This changellowed for more efficient use of the SIMD arithmetic unit and cachememory in the CPU. Additionally, the data structure was modified sothat intermediate data generated during computation is retained inemory rather than being written to a file, which reduced file I/O time.hese efficiency improvements resulted in approximately a twofoldpeed-up. Furthermore, by minimizing the number of computing nodessed in each case, it became possible to run a large number of casessimultaneously.3. ResultsFirst, the optimized parameters of the tetragonal supercell, alongwith the resulting LDA energy gap 𝛥𝜀LDA and QP (𝐺 𝑊 ) energy gap𝜀QP, are listed in Table 1. The QP energy gap of pure rutile is 3.343 Table 1Optimized parameters of the tetragonal supercell and resulting LDA gap 𝛥𝜀LDA and QP(𝐺 𝑊 ) energy gap 𝛥𝜀QP.System 𝑎 = 𝑏 𝑐 𝛥𝜀LDA 𝛥𝜀QPpure rutile 9.1453 Å 8.7767 Å 1.78 eV 3.34 eVZr@rutile 9.1623 Å 8.8197 Å 1.80 eV 3.39 eVMo@rutile 9.1654 Å 8.7718 Å 0.30 eV 0.97 eVZn@rutile 9.1499 Å 8.7958 Å 0.40 eV 1.08 eVCd@rutile 9.1801 Å 8.8269 Å 0.24 eV 0.810 eVZn+OV@rutile 9.1363 Å 8.8223 Å 1.78 eV 3.37 eVCd+OV@rutile 9.1449 Å 8.8639 Å 1.76 eV 3.39 eVeV, which is consistent with previous 𝐺 𝑊 results of 3.3 eV [66,67],3.34 eV [70], and 3.59 eV [72], as well as the experimental value of.3 ± 0.5 eV [73]. For Zr doped rutile, both the LDA and QP energy gapsare nearly the same as those of pristine rutile. This finding aligns witha previous study [21] that indicated minimal changes in the LDA bandstructure due to the substitution of Zr, which is expected given that Zrand Ti are both group 4 elements. In contrast, all other systems dopedwith Mo, Zn, or Cd exhibit lower LDA and QP energy gaps. Our resultfor Mo doped rutile is consistent with the DFT-level spectra calculatedby Soussi et al. [26].Next, the resulting singlet PAS are shown in Fig. 1. The low-energyregions of the PAS for (a) pure rutile, (c) Zr@rutile, and (e) Mo@rutileare magnified in the right panels (b), (d) and (f), respectively. For purerutile (Fig. 1(a), (b)), the photoabsorption begins around 3 eV, whichaligns well with the experimental evidence [74]. Although there is aside peak at approximately 3.75 eV and no peak around 8.5 eV inFig. 1(b) in contrast to the previous 𝐺 𝑊 + BSE calculations [70,75,76], the overall behavior is in good agreement with the experimentalPAS [77]. Similarly, for Zr doped rutile (Fig. 1(c), (d)), the pho-toabsorption starting around 3 eV is consistent with the experimentalabsorption edge [19]. In the case of Mo (Fig. 1(e), (f)), the low-energypikes at 1.5 eV and around 2 and 3 eV in the PAS of Fig. 1(f)ay contribute to the long absorption tail and the observed red shiftf the absorption edge by approximately 0.5 eV in the experimentalhotoabsorption spectra [23].In contrast, for Zn and Cd, large peaks appear around 1 eV inFig. 1(g) and (h). However, it is noteworthy that there is no such evi-dence in the experimental PAS, although small red shifts of 0.05–0.1 eVhave been reported for the absorption edge in the experiments [27–30]. We found that this significant discrepancy arises because Zn andCd impurity atoms must coexist with OV to compensate for the chargedifference between Zn2+/Cd2+ and Ti4+. To address this, we introducedan OV in the vicinity of the Zn/Cd atom in the supercell and performeda full optimization of the system. The relaxed supercell parameters,along with the resulting LDA and 𝐺 𝑊 energy gaps, are listed in the lasttwo rows of Table 1. The LDA and 𝐺 𝑊 energy gaps for these systemsare clearly larger than those for systems without OV, suggesting thetability of the Zn/Cd-OV pair in rutile TiO2. The singlet PAS for thesesystems are shown in Fig. 2. The low energy portions of the PAS for (a)Zn+OV@rutile and (c) Cd+OV@rutile are magnified in the right panels(b) and (d). As expected, the peaks around 1 eV have disappeared,and the overall behavior has become similar to that of pristine rutile.Now, the absorption intensity below 4 eV decreases smoothly towardzero, aligning with the experimental PAS [27–30]. The small red shifteported for the absorption edge may be attributed to an imperfectoexistence of Zn/Cd and OV. Specifically, this small red shift coulde observed in the absorption edge if some Zn/Cd impurities withoutV were also present in the sample.4. DiscussionFor pure rutile, we examined the convergence of results with respectto the supercell size. We performed 𝐺 𝑊 + BSE calculations usingseveral different supercells, ranging from 2 × 1 × 1 to 2 × 2 × 4,K. Ohno et al. Computational Condensed Matter 41 (2024) e00977 Fig. 1. Resulting PAS of (a), (b) pure rutile, (c), (d) Zr@rutile, (e), (f) Mo@rutile, (g) Zn@rutile, and (h) Cd@rutile. Panels (b), (d) and (f) are magnified plots of (a), (c) and (e).while maintaining consistent cutoff energies and the number of levelsas used in the 2 × 2 × 3 supercell. (Note that the 2 × 2 × 4 supercellutilized a slightly lower cutoff energy of 69.9 Ry for ⟨𝜙𝑛 |𝛴𝑥 |𝜙𝑛 ⟩ and𝑈𝑣𝑐;𝑣′𝑐′ due to computing resource limitations.) The resulting LDA andQP energy gaps, along with the first three lower PAEs (i.e., the energyeigenvalues from the singlet Eq. (5) of the BSE), are shown in Table 2.As indicated in this table, both the energy gaps and PAEs approachcertain values as the supercell size increases, although the changes arenot so significant. This suggests that the 2 × 2 × 3 supercell is sufficientwithin the accuracy of 0.1 eV, yielding a QP (𝐺 𝑊 ) energy gap 𝛥𝜀QP ofapproximately 3.3 eV.In addition to this check, we also evaluated the 𝐺 𝑊 calculationof the 6-atom primitive cell using k-point sampling with various k-point meshes, ranging from 2 × 1 × 1 to 3 × 3 × 6 for pure rutile.The resulting LDA and QP energy gaps are presented in Table 3.4 The dependence of the energy gaps on the k-point grid is notablysignificant and scattered compared to their dependence on supercellsize in Table 2. The convergence with respect to the k-point mesh isslow, and the 3 × 3 × 6 mesh does not appear to be fully converged.Indeed, previous studies [67,70,72,75,76] have utilized 4 × 4 × 6 k-points or more for better convergence. In our prior work [67] usingthe same TOMBO code with a 5 × 5 × 8 k-point grid, we obtained aQP energy gap 𝛥𝜀QP of 3.30 eV.Below, we will discuss the validity of using only the 𝛤 -point inthe 𝐺 𝑊 + BSE approach in more detail. The primary condition forasserting the validity of this approach is that the 𝛤 -point only 𝐺 𝑊calculation in the large supercell should reasonably reproduce theavailable experimental and theoretical band gaps. This criterion isessential for judging whether the 𝛤 -point only 𝐺 𝑊 + BSE approachis applicable in a supercell system. Next, we consider the specifics ofK. Ohno et al. Computational Condensed Matter 41 (2024) e00977 Fig. 2. Resulting PAS of (a), (b) Zn+OV@rutile and (c), (d) Cd+OV@rutile. Panels (b) and (d) are magnified plots of (a) and (c).Table 2Change in the LDA gap 𝛥𝜀LDA, QP (𝐺 𝑊 ) energy gap 𝛥𝜀QP, and first three lower PAEs(i.e., the energy eigenvalues of the singlet Eq. (5) of the BSE) with respect to thesupercell size.Supercell 𝛥𝜀LDA (eV) 𝛥𝜀QP (eV) PAEs (eV)2 × 1 × 1 1.46 2.92 2.85 2.92 3.002 × 1 × 2 1.57 3.31 3.28 3.37 3.372 × 1 × 3 1.74 3.49 3.47 3.61 3.802 × 1 × 4 1.72 3.24 3.22 3.22 3.352 × 2 × 1 1.61 3.42 3.23 3.25 3.252 × 2 × 2 1.61 3.21 3.17 3.20 3.212 × 2 × 3 1.78 3.34 3.32 3.37 3.372 × 2 × 4 1.77 3.29 3.18 3.18 3.18Table 3Change in the LDA gap 𝛥𝜀LDA and QP (𝐺 𝑊 ) energy gap 𝛥𝜀QP (in units of eV) withrespect to the k-points.k-points 𝛥𝜀LDA 𝛥𝜀QP k-points 𝛥𝜀LDA 𝛥𝜀QP k-points 𝛥𝜀LDA 𝛥𝜀QP2 × 1 × 1 1.61 3.40 2 × 2 × 1 1.67 3.70 3 × 3 × 1 1.45 3.772 × 1 × 2 1.83 3.87 2 × 2 × 2 1.87 4.09 3 × 3 × 2 1.82 4.062 × 1 × 3 1.52 3.10 2 × 2 × 3 1.52 2.91 3 × 3 × 3 1.45 3.142 × 1 × 4 1.73 3.48 2 × 2 × 4 1.72 3.45 3 × 3 × 4 1.63 3.382 × 1 × 5 1.68 3.26 2 × 2 × 5 1.72 3.36 3 × 3 × 5 1.69 3.342 × 1 × 6 1.73 3.52 2 × 2 × 6 1.74 3.45 3 × 3 × 6 1.70 3.38Table 4Diagonal elements between 𝑣 = VBM and 𝑐 = CBM of direct Coulomb term 𝑈𝑣𝑐;𝑣𝑐 ,direct correlation term 𝑆𝑣𝑐;𝑣𝑐 , direct (screened Coulomb) term 𝑊𝑣𝑐;𝑣𝑐 = 𝑈𝑣𝑐;𝑣𝑐 + 𝑆𝑣𝑐;𝑣𝑐 ,exchange term 𝑋𝑣𝑐;𝑣𝑐 , and ‘‘superficial’’ exciton binding energy 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 =𝑈𝑣𝑐;𝑣𝑐 + 𝑆𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 in units of eV.System 𝑈𝑣𝑐;𝑣𝑐 𝑆𝑣𝑐;𝑣𝑐 𝑊𝑣𝑐;𝑣𝑐 𝑋𝑣𝑐;𝑣𝑐 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐pure rutile −0.0377 0.0212 −0.0165 0.0018 −0.0201Zr@rutile −0.0316 0.0190 −0.0126 0.0023 −0.0172Mo@rutile 1.8773 −1.4117 0.4656 0.1395 0.1866Zn@rutile 0.1880 −0.0687 0.1193 0.1253 −0.1313Cd@rutile 0.2349 −0.0970 0.1379 0.0957 −0.0535Zn+OV@rutile −0.0629 0.0441 −0.0188 0.0036 −0.0260Cd+OV@rutile −0.0624 0.0442 −0.0182 0.0044 −0.0270the BSE calculation using only the 𝛤 -point. In the singlet BSE (5), themost important matrix elements are the diagonal elements between5 Table 5First eight lower PAEs, i.e., the energy eigenvalues of the singlet Eq. (5) of theBSE.System PAEs (eV)pure rutile 3.32 3.37 3.37 3.57 3.64 3.64 3.65 3.65Zr@rutile 3.36 3.37 3.41 3.59 3.60 3.62 3.63 3.65Mo@rutile 0.71 0.75 1.32 1.42 1.44 1.46 1.50 1.54Zn@rutile 0.91 1.08 1.10 1.21 1.53 1.54 1.54 1.58Cd@rutile 0.808 0.89 0.98 1.01 1.28 1.42 1.50 1.53Zn+OV@rutile 3.35 3.42 3.45 3.50 3.56 3.56 3.64 3.66Cd+OV@rutile 3.36 3.41 3.45 3.48 3.55 3.58 3.60 3.64𝑣 = VBM and 𝑐 = CBM for the direct Coulomb term 𝑈𝑣𝑐;𝑣𝑐 (Eq. (11)),the direct correlation term 𝑆𝑣𝑐;𝑣𝑐 (Eq. (10)), and the exchange term𝑋𝑣𝑐;𝑣𝑐 (Eq. (12)). These values are summarized in Table 4. Additionally,the direct (screened Coulomb) term 𝑊𝑣𝑐;𝑣𝑐 = 𝑈𝑣𝑐;𝑣𝑐 + 𝑆𝑣𝑐;𝑣𝑐 , and the‘‘superficial’’ exciton binding energy 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 = 𝑈𝑣𝑐;𝑣𝑐 + 𝑆𝑣𝑐;𝑣𝑐 −2𝑋𝑣𝑐;𝑣𝑐 are presented in the fourth and sixth (last) columns of this table.For the Mo@rutile system, 𝑈𝑣𝑐;𝑣𝑐 is a large positive value, while𝑆𝑣𝑐;𝑣𝑐 is a moderate negative value. Consequently, the screenedelectron–hole interaction 𝑊𝑣𝑐;𝑣𝑐 = 𝑈𝑣𝑐;𝑣𝑐 + 𝑆𝑣𝑐;𝑣𝑐 remains positive. Onthe other hand, 𝑋𝑣𝑐;𝑣𝑐 is small and positive. As a result, the ‘‘superficial’’exciton binding energy 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 is also positive, making theVBM–CBM diagonal element of 𝐷𝑣𝑐;𝑣𝑐 + 2𝑋𝑣𝑐;𝑣𝑐 in the singlet Eq. (5)(i.e., ‘‘superficial’’ PAE) smaller than the QP energy gap 𝛥𝜀QP = 𝜀QP𝑐 −𝜀QP𝑣 . Indeed, as shown in Table 5, the lower PAEs, i.e., lower energyeigenvalues of the singlet Eq. (5), for Mo@rutile are smaller than theQP energy gap 𝛥𝜀QP presented in Table 1. The ‘‘real’’ exciton bindingenergy, which can be defined as 𝛥𝜀QP− (the lowest PAE), is listed inTable 6. This value is 0.26 eV in the case of the Mo@rutile system.This explains the appearance of peaks at small PAEs in Fig. 1(f) for theMo impurity case.Conversely, for pure rutile, Zr@rutile, Zn+OV@rutile, andCd+OV@rutile, 𝑈𝑣𝑐;𝑣𝑐 is marginally negative, while 𝑆𝑣𝑐;𝑣𝑐 is smallbut positive, and 𝑋𝑣𝑐;𝑣𝑐 is very small and positive. Consequently, thescreened electron–hole interaction 𝑊𝑣𝑐;𝑣𝑐 = 𝑈𝑣𝑐;𝑣𝑐 + 𝑆𝑣𝑐;𝑣𝑐 and the‘‘superficial’’ exciton binding energy 𝑊𝑣𝑐;𝑣𝑐− 2𝑋𝑣𝑐;𝑣𝑐 are both marginallynegative. However, the lowest PAE listed in Table 5 is at least slightlysmaller than the QP energy gap 𝛥𝜀QP in Table 1. Therefore, as shownK. Ohno et al.sh2(Titk𝑊tiiucpppPaspatepceebTbtcComputational Condensed Matter 41 (2024) e00977 Table 6‘‘real’’ exciton binding energy, i.e., 𝛥𝜀QP− (the lowest PAE).System 𝛥𝜀QP− (the lowest PAE) (eV)pure rutile 0.02Zr@rutile 0.03Mo@rutile 0.26Zn@rutile 0.17Cd@rutile 0.002Zn+OV@rutile 0.02Cd+OV@rutile 0.03in Table 6, the ‘‘real’’ exciton binding energy, 𝛥𝜀QP− (the lowest PAE),is positive, even if it is very small.The Zn@rutile and Cd@rutile systems behave somewhere betweenthese two extremes. As shown in Table 4, 𝑈𝑣𝑐;𝑣𝑐 is a positive butmall, while 𝑆𝑣𝑐;𝑣𝑐 is marginally negative. Thus, the screened electron–ole interaction 𝑊𝑣𝑐;𝑣𝑐 remains positive, although small. If we subtract𝑋𝑣𝑐;𝑣𝑐 from this 𝑊𝑣𝑐;𝑣𝑐 , the resulting ‘‘superficial’’ exciton bindingenergy becomes marginally negative similar to the cases for pure rutile,Zr@rutile, Zn+OV@rutile, and Cd+OV@rutile. Similarly, as shown inTable 6, the ‘‘real’’ exciton binding energy becomes positive, althoughthe value is quite small for Cd@rutile. However, the origin of the lowerenergy spikes in Fig. 1(g) and (h) is different from the VBE (𝑣) and CBE𝑐) levels and related to the levels far away from the VBE and CBE.his is, of course, an exceptional case, due to the inherent chemicalnstability of these systems.The marginally negative values of the screened electron–hole in-eraction 𝑊𝑣𝑐;𝑣𝑐 may also arise when 𝑣 and 𝑐 correspond to the same-point in k-point samplings. The vanishing 𝑮 = 0 component of𝑣𝑐;𝑣𝑐 is the reason for its potential negativity. Typically, the Fourierransform of the bare Coulomb interaction 𝑣(𝑮) = 4𝜋∕𝛺|𝑮|2 (where 𝛺s the volume of the supercell) is set to a nonzero value at 𝑮 = 0 [78–80]. Then, in a plasmon pole model [44,81,82], the 𝑮 = 0 componentof the screened electron–hole interaction 𝑊𝑣𝑐;𝑣𝑐 automatically vanishes.When we express 𝑊𝑣𝑐;𝑣𝑐 as 𝑊𝑣𝑐;𝑣𝑐 = 𝑈𝑣𝑐;𝑣𝑐 +𝑆𝑣𝑐;𝑣𝑐 (see Eq. (10)), the𝑮 = 0 component of 𝑆𝑣𝑐;𝑣𝑐 in a plasmon pole model also vanishes for thesame reason; see Appendix. The 𝑮 = 0 component of the bare Coulombnteraction 𝑈𝑣𝑐;𝑣𝑐 must similarly vanish in the following reason. Insidethe unit cell, the 𝑣𝑣 and 𝑐 𝑐 charges in Eq. (11) are exactly equal tonity:∫ 𝜙∗𝑐 (𝒓)𝜙𝑐 (𝒓) 𝑑𝒓 = ∫ 𝜙𝑣(𝒓)𝜙∗𝑣(𝒓) 𝑑𝒓 = 1. (13)Considering their signs, the 𝑣𝑣 and 𝑐 𝑐 charges are plus and minus one,respectively. Their images repeat across all unit cells, ensuring overallharge neutrality for the whole system. Therefore, the 𝑮 = 0 componentof their Coulomb interaction 𝑈𝑣𝑐;𝑣𝑐 must be set to zero. (In the 𝛤 -pointonly calculation, if we sum up all 𝑐 𝑐 charges independently of all 𝑣𝑣charges over all unit cells, the Coulomb integral will obviously diverge.Instead, to take into account the excitonic effect correctly, we have tocount not only the interaction between the 𝑐 𝑐 and 𝑣𝑣 charges but alsothe interaction between 𝑐 𝑐 and 𝑐 𝑐 (or 𝑣𝑣 and 𝑣𝑣) charges in differentunit cells, at least in the 𝐺 = 0 limit. Therefore, we can assume chargeneutrality in this calculation. Note that this issue does not arise in k-oint sampling, because 𝑐 𝑐 and 𝑣𝑣 charges are only assumed in thearticular unit cell under consideration.) The 𝑮 = 0 component of theexchange term 𝑋𝑣𝑐;𝑣𝑐 also vanishes, because∫ 𝜙∗𝑐 (𝒓)𝜙𝑣(𝒓) 𝑑𝒓 = ∫ 𝜙𝑐 (𝒓′)𝜙∗𝑣(𝒓′) 𝑑𝒓′ = 0 (14)clearly holds in the 𝑮 = 0 limit (see Appendix).In the context of k-point samplings, marginally negative values ofthe ‘‘superficial’’ exciton binding energy (𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐) can lead to aositive value for the ‘‘real’’ exciton binding energy (𝛥𝜀QP− the lowestAE) as follows. Consider, for example, the inclusion of the 𝑋-pointlongside the 𝛤 -point, assuming no mixing between different levels forimplicity. In this case, the singlet Eq. (5) becomes o6 [𝐷𝑣𝛤 𝑐𝛤 ;𝑣𝛤 𝑐𝛤 + 2𝑋𝑣𝛤 𝑐𝛤 ;𝑣𝛤 𝑐𝛤 −𝑊𝑣𝛤 𝑐𝛤 ;𝑣𝑋 𝑐𝑋 + 2𝑋𝑣𝛤 𝑐𝛤 ;𝑣𝑋 𝑐𝑋−𝑊𝑣𝑋 𝑐𝑋 ;𝑣𝛤 𝑐𝛤 + 2𝑋𝑣𝑋 𝑐𝑋 ;𝑣𝛤 𝑐𝛤 𝐷𝑣𝑋 𝑐𝑋 ;𝑣𝑋 𝑐𝑋 + 2𝑋𝑣𝑋 𝑐𝑋 ;𝑣𝑋 𝑐𝑋][𝐴 (s)(𝑣𝛤 , 𝑐𝛤 )𝐴 (s)(𝑣𝑋 , 𝑐𝑋 )]= 𝛺 (s)±[𝐴 (s)(𝑣𝛤 , 𝑐𝛤 )𝐴 (s)(𝑣𝑋 , 𝑐𝑋 )]. (15)If 𝐷𝑣𝑐;𝑣𝑐 = 𝛥𝜀QP −𝑊𝑣𝑐;𝑣𝑐 and 𝑋𝑣𝑐;𝑣𝑐 are independent of the 𝛤 - and 𝑋-oints (i.e., behave as constants), the secular equation can be expresseds(𝛥𝜀QP −𝑊𝑣𝑐;𝑣𝑐 + 2𝑋𝑣𝑐;𝑣𝑐 −𝛺(s)±)2− |||−𝑊𝑣𝛤 𝑐𝛤 ;𝑣𝑋 𝑐𝑋 + 2𝑋𝑣𝛤 𝑐𝛤 ;𝑣𝑋 𝑐𝑋|||2= 0,(16)which yields the following ± solutions:𝛺 (s)± = 𝛥𝜀QP −𝑊𝑣𝑐;𝑣𝑐 + 2𝑋𝑣𝑐;𝑣𝑐 ±|||−𝑊𝑣𝛤 𝑐𝛤 ;𝑣𝑋 𝑐𝑋 + 2𝑋𝑣𝛤 𝑐𝛤 ;𝑣𝑋 𝑐𝑋|||(17)Therefore, when 𝑊𝑣𝑐;𝑣𝑐 is marginally negative, the − solution can offsethe negative ‘‘superficial’’ exciton binding energy, shifting the ‘‘real’’xciton binding energy to a positive value. This effect becomes moreronounced with the inclusion of additional k-points, of course.Now, we demonstrate that a similar explanation applies to supercellalculations without k-point samplings. Assume that a unit cell (eitherprimitive or supercell) is doubled in one direction, and the number oflevels at the 𝛤 -point is also doubled. If we ignore all the 𝑣-𝑣′ and 𝑐-𝑐′interlevel interactions except for the interaction between CBM (0) andCBM+1 (1), the singlet Eq. (5) becomes[𝐷𝑣𝑐0;𝑣𝑐0 + 2𝑋𝑣𝑐0;𝑣𝑐0 −𝑊𝑣𝑐0;𝑣𝑐1 + 2𝑋𝑣𝑐0;𝑣𝑐1−𝑊𝑣𝑐1;𝑣𝑐0 + 2𝑋𝑣𝑐1;𝑣𝑐0 𝐷𝑣𝑐1;𝑣𝑐1 + 2𝑋𝑣𝑐1;𝑣𝑐1][𝐴 (s)(𝑣, 𝑐0)𝐴 (s)(𝑣, 𝑐1)]= 𝛺 (s)±[𝐴 (s)(𝑣, 𝑐0)𝐴 (s)(𝑣, 𝑐1)]. (18)Again, if 𝐷𝑣𝑐;𝑣𝑐 = 𝛥𝜀QP −𝑊𝑣𝑐;𝑣𝑐 and 𝑋𝑣𝑐;𝑣𝑐 are independent of 𝑐0 and 𝑐1,and behave as constants, the secular equation is obtained as(𝛥𝜀QP −𝑊𝑣𝑐;𝑣𝑐 + 2𝑋𝑣𝑐;𝑣𝑐 −𝛺(s)±)2− |||−𝑊𝑣𝑐0;𝑣𝑐1 + 2𝑋𝑣𝑐0;𝑣𝑐1|||2= 0, (19)which results in the following ± solutions:𝛺 (s)± = 𝛥𝜀QP −𝑊𝑣𝑐;𝑣𝑐 + 2𝑋𝑣𝑐;𝑣𝑐 ±|||−𝑊𝑣𝑐0;𝑣𝑐1 + 2𝑋𝑣𝑐0;𝑣𝑐1|||. (20)Thus, even if the screened electron–hole interaction 𝑊𝑣𝑐;𝑣𝑐 is marginallynegative, the − solution can compensate for the negative ‘‘superficial’’xciton binding energy, potentially shifting the ‘‘real’’ exciton bindingnergy (𝛥𝜀QP − 𝛺 (s)− ) to a positive value. This effect becomes moreprominent with larger supercell. Therefore, the marginally negativevalues of 𝑊𝑣𝑐;𝑣𝑐 and the ‘‘superficial’’ exciton binding energy need note a major concern.We present the resulting exciton wave functions (EWFs) forMo@rutile at PAEs of 0.71 eV and 1.54 eV in Fig. 3. The position ofthe hole is fixed near the impurity atom. As illustrated in the figures,the EWFs are clearly localized around the impurity. This supercellcalculation can be likened to using some k-point mesh in the calculationof the six-atom primitive unit cell. If this analogy holds, similar EWFswould also emerge in such k-point samplings. In general, to fullyconfine EWFs within an 𝑛1 × 𝑛2 × 𝑛3 supercell, it would be necessary toeither increase the supercell size or the number of the k-points in theprimitive unit cell, both of which are computationally demanding tasks.The spatial extension of EWFs is related to the ‘‘real’’ exciton bindingenergy. For Mo@rutile, the ‘‘real’’ exciton binding energy 𝛥𝜀QP− (thelowest PAE) = 0.26 eV (see Table 6), which is reasonably large.his value is further greater than the positive ‘‘superficial’’ excitoninding energy 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 = 0.19 eV in Table 4, indicating thatit corresponds to a Frenkel exciton.On the other hand, for pure@rutile, Zr@rutile, Zn+OV@rutile, andCd+OV@rutile, the ‘‘real’’ exciton binding energies (𝛥𝜀QP− the lowestPAE) listed in Table 6 are very small. These values fall well belowhe accuracy threshold of the 𝐺 𝑊 + BSE calculation. Therefore, wean reasonably expect that there is no issue with depicting the PASf these systems using only the 𝛤 -point calculation, even though theK. Ohno et al. Computational Condensed Matter 41 (2024) e00977 Fig. 3. Exciton wave functions (EWFs) of (a) Mo@rutile at PAE = 0.71 eV and (b) Mo@rutile at PAE = 1.54 eV. The position of hole is fixed at the vicinity of the impurity atom.Small sky blue, purple and red balls represent, respectively, Ti, Mo, and O atoms. Yellow and blue regions correspond, respectively, to plus and minus sign of the EWF. (Forinterpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)EWF may correspond to a widely spread Wannier exciton that cannotbe fully confined within the supercell. A similar scenario would arisein primitive cell calculations using k-point sampling. When the ‘‘real’’exciton binding energy is smaller than the error bar of the 𝐺 𝑊 +BSE calculation using 𝑛1 × 𝑛2 × 𝑛3 k-points, one can still accuratelydepict the PAS in k-point sampling, even if the EWF is not confinedwithin the unfolded 𝑛1 × 𝑛2 × 𝑛3 supercell. According to Sham andRice [83], this situation relates to the effective-mass equation for theWannier exciton and corresponds to a weak electron–hole interaction.It is crucial to confirm that the ‘‘real’’ exciton binding energy (𝛥𝜀QP−the lowest PAE) is positive in both the 𝛤 -point only BSE calculationand the BSE calculation utilizing k-point sampling. Even if this value isexceptionally small, we can still compute the PAS reliably.Finally, we would like to comment on the calculations involvingmagnetic impurities. We attempted to calculate the PAS of Cr@rutileand Ru@rutile but encountered difficulties due to the emergence ofnegative PAEs in the solution of Eq. (5) of the BSE. This failure is likelyattributed to our neglect of spin polarization effects for Cr and Ru.Previous DFT studies [31,84,85] have highlighted the significance ofspin polarization effects for these impurities. Therefore, the 𝐺 𝑊 + BSEcalculation for magnetic impurities will be addressed in future work.5. ConclusionWe have employed the 𝛤 -point only 𝐺 𝑊 + BSE approach on a 72(or 71) atom supercell to obtain the quasiparticle (QP) energy gapsand photoabsorption spectra (PAS) of pure rutile TiO2 and those dopedwith nonmagnetic transition metal (Zr, Mo, Zn, Cd) impurities. The full𝜔 integration was performed at the 𝐺 𝑊 level, while the generalizedplasmon pole (GPP) model was utilized at the BSE level. In Appendix,we identified a typo in the key expression of the screened electron–holeinteraction 𝑊𝑣𝑐;𝑣′𝑐′ in the BSE in the pioneering papers [71]. Duringour calculations, we found that Zn and Cd impurities should coexistwith oxygen vacancies (OVs), as Zn2+ and Cd2+ are unstable at the Ti4+substitutional site in TiO2 without OVs. Among the systems studied,Mo, Zn, or Cd doped rutile TiO2 demonstrate visible light absorbanceand may exhibit catalytic activity under visible light. The resulting QPenergies and PAS align roughly well with available experimental andtheoretical data. We thoroughly discussed the validity of our approachand summarized the necessary conditions for the applicability of usingonly the 𝛤 -point in the 𝐺 𝑊 + BSE framework as follows: First, the𝛤 -point only 𝐺 𝑊 calculation in the assumed supercell should reason-ably reproduce the available experimental and theoretical band gaps.7 Second, the ‘‘superficial’’ exciton binding energy — i.e., the diagonalelement of 𝑊𝑣𝑐;𝑣𝑐 − 2𝑋𝑣𝑐;𝑣𝑐 between 𝑣 = VBM and 𝑐 = CBM — must bepositive or marginally negative. Third, physically, the ‘‘real’’ excitonbinding energy (𝛥𝜀QP− the lowest PAE) should be positive, even ifit is exceptionally small. The 𝛤 -point only 𝐺 𝑊 + BSE approach caneffectively handle a wide range of periodic systems, including thosewith impurities and so on, while saving computational time.CRediT authorship contribution statementKaoru Ohno:Writing – original draft, Validation, Software, Method-ology, Investigation, Funding acquisition, Formal analysis. Ryoji Sa-hara: Writing – review & editing, Validation, Resources, Project ad-ministration, Investigation, Formal analysis. Takeshi Nanri: Writing– original draft, Software, Resources, Methodology. Yoshiyuki Kawa-zoe: Writing – review & editing, Validation, Project administration,Funding acquisition, Conceptualization.Declaration of competing interestThe authors declare the following financial interests/personal rela-tionships which may be considered as potential competing interests:Kaoru Ohno reports financial support was provided by Japan Societyfor the Promotion of Science. Yoshiyuki Kawazoe reports financialsupport was provided by Air Force Asian Office of Aerospace Researchand Development. If there are other authors, they declare that theyhave no known competing financial interests or personal relationshipsthat could have appeared to influence the work reported in this paper.AcknowledgmentsIn this study, we utilized supercomputers at the Research Institutefor Information Technology in Kyushu University (supported by the JH-PCN project No.: jh240026), the Institute for Materials Research (IMR)in Tohoku University (Subject Nos.: 2211SC0501 and 2308SC0216),and the Institute for Solid State Physics in the University of Tokyo.One of the authors (K.O.) received support from the Japan Societyfor the Promotion of Science (JSPS) through KAKENHI Grants-in-Aidfor Scientific Research B (Grant Nos. 21H01877 and 23K21094). Theauthors would like to thank Dr. Ruth Pachter from the Air ForceResearch Laboratory for suggesting this study and providing valuableinsights. This work was also partially supported by the Asian Office ofAerospace Research and Development (AOARD) under project numberK. Ohno et al.(ott𝑊ttCbFEh𝜏t𝜒Computational Condensed Matter 41 (2024) e00977 FA2386-22-1-4011. Y.K. acknowledges support from (i) the NationalScience, Research and Innovation Fund (NSRF) (NRIIS Project Number90465), (ii) Thailand Science Research and Innovation (TSRI), and (iii)Suranaree University of Technology (SUT).Appendix. Derivation of 𝑾𝒗𝒄;𝒗′𝒄′One-particle and two-particle Green’s functions are defined, respec-tively, as [86]𝐺(𝑥1, 𝑥2) = −𝑖⟨𝛷𝑁𝐺 | 𝑇 [𝜓(𝑥1)𝜓†(𝑥2)] |𝛹𝑁𝐺 ⟩, (A.1)𝐺(𝑥1, 𝑥2; 𝑥′1, 𝑥′2) = (−𝑖)2⟨𝛷𝑁𝐺 | 𝑇 [𝜓(𝑥1)𝜓(𝑥2)𝜓†(𝑥′2)𝜓†(𝑥′1)] |𝛹𝑁𝐺 ⟩, (A.2)where 𝑇 means the time ordered product; 𝑥𝑖 and 𝑥′𝑖 are the 4-vectors𝒓𝑖, 𝑡𝑖) and (𝒓′𝑖 , 𝑡′𝑖); 𝜓(𝑥) and 𝜓†(𝑥) represent the annihilation and creationperators in the Heisenberg representation. The two-particle correla-ion function is defined as𝐿(𝑥1, 𝑥2; 𝑥′1, 𝑥′2) = − 𝐺(𝑥1, 𝑥2; 𝑥′1, 𝑥′2) + 𝐺(𝑥1, 𝑥′1)𝐺(𝑥2, 𝑥′2), (A.3)which excludes completely disconnected diagrams from the two-particle Green’s function. It satisfies the Bethe–Salpeter equation (BSE)𝐿(𝑥1, 𝑥2; 𝑥′1, 𝑥′2) = 𝐺(𝑥1, 𝑥′2)𝐺(𝑥2, 𝑥′1)+ ∫ 𝐺(𝑥1, 𝑥3)𝛯(𝑥3, 𝑥4; 𝑥′3, 𝑥′4)𝐺(𝑥′3, 𝑥′1)× 𝐿(𝑥′4, 𝑥2; 𝑥4, 𝑥′2)𝑑 𝑥3𝑑 𝑥′3𝑑 𝑥4𝑑 𝑥′4. (A.4)From the 𝐺 𝑊 self-energy 𝛴, the kernel 𝛯 = 𝛿 𝛴∕𝛿 𝐺 is𝛯(𝑥3, 𝑥4; 𝑥′3, 𝑥′4) = − 𝑖𝛿4(𝑥3 − 𝑥′3)𝛿4(𝑥4 − 𝑥′+4 )𝑣(𝑥3, 𝑥4)+ 𝑖𝛿4(𝑥3 − 𝑥′4)𝛿4(𝑥′+3 − 𝑥4)𝑊 ′(𝑥3, 𝑥′3), (A.5)where 𝑣(𝑥3, 𝑥4) = 𝛿(𝑡3 − 𝑡4)∕|𝒓3 − 𝒓4| is the exchange term coming fromhe functional derivative of the Hartree term of the self-energy and′(𝑥3, 𝑥′3) is the direct term coming from the functional derivative ofhe exchange–correlation component of the self-energy. 𝑊 ′ includeshe second exchange term 𝐺 𝛿 𝑊 ∕𝛿 𝐺 ∝ 𝐺 𝑊 𝑊 𝐺 [58,60], but thisterm is usually ignored. Then, 𝑊 ′ reduces to the dynamically screenedoulomb interaction 𝑊 .In cases where 𝑡1, 𝑡′1 ≶ 𝑡2, 𝑡′2 hold, the two-particle Green’s function(A.2) becomes𝐺(𝑥1, 𝑥2; 𝑥′1, 𝑥′2)= − ⟨𝛷𝑁𝐺 | 𝑇 [𝜓(𝑥1)𝜓†(𝑥′1)]𝑇 [𝜓(𝑥2)𝜓†(𝑥′2)] |𝛹𝑁𝐺 ⟩𝜃(min(𝑡1, 𝑡′1) − max(𝑡2, 𝑡′2))− ⟨𝛷𝑁𝐺 | 𝑇 [𝜓(𝑥2)𝜓†(𝑥′2)]𝑇 [𝜓(𝑥1)𝜓†(𝑥′1)] |𝛹𝑁𝐺 ⟩𝜃(min(𝑡2, 𝑡′2) − max(𝑡1, 𝑡′1)).(A.6)By inserting the complete set condition ∑𝑆 |𝛹𝑁𝑆 ⟩⟨𝛹𝑁𝑆 | = 1 for the𝑁-electron intermediate states |𝛹𝑁𝑆 ⟩, which represent configurationswith one particle and one hole (i.e., an exciton) excited from the groundstate, in between two 𝑇 products, the ⟨𝛷𝑁𝐺 | ... |𝛹𝑁𝐺 ⟩ terms in Eq. (A.6)ecome∑𝑆⟨𝛹𝑁𝐺 | 𝑇 [𝜓(𝑥1)𝜓†(𝑥′1)] |𝛹𝑁𝑆 ⟩⟨𝛹𝑁𝑆 | 𝑇 [𝜓(𝑥2)𝜓†(𝑥′2)] |𝛹𝑁𝐺 ⟩, (A.7a)∑𝑆⟨𝛹𝑁𝐺 | 𝑇 [𝜓(𝑥2)𝜓†(𝑥′2)] |𝛹𝑁𝑆 ⟩⟨𝛹𝑁𝑆 | 𝑇 [𝜓(𝑥1)𝜓†(𝑥′1)] |𝛹𝑁𝐺 ⟩. (A.7b)Following Strinati [51], we introduce the averaged and differencetimes as 𝑡1 = (𝑡1 + 𝑡′1)∕2, 𝑡2 = (𝑡2 + 𝑡′2)∕2, 𝜏1 = 𝑡1 − 𝑡′1, and 𝜏2 = 𝑡2 − 𝑡′2.rom⟨𝛹𝑁𝐺 | 𝑇 [𝜓(𝑥1)𝜓†(𝑥′1)] |𝛹𝑁𝑆 ⟩ =[⟨𝛹𝑁𝐺 |𝜓(𝒓1)𝑒−𝑖𝜏1𝜓†(𝒓′1) |𝛹𝑁𝑆 ⟩𝜃(𝜏1)− ⟨𝛹𝑁𝐺 |𝜓†(𝒓′1)𝑒𝑖𝜏1𝜓(𝒓1) |𝛹𝑁𝑆 ⟩𝜃(−𝜏1)]𝑒𝑖[(𝐸𝐺−𝐸𝑆 )𝑡1+(𝐸𝐺+𝐸𝑆 )|𝜏1|∕2]= 𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1)𝑒𝑖(𝐸𝐺−𝐸𝑆 )𝑡1 , (A.8a)⟨𝛹𝑁𝑆 | 𝑇 [𝜓(𝑥2)𝜓†(𝑥′2)] |𝛹𝑁𝐺 ⟩ =[⟨𝛹𝑁𝑆 |𝜓(𝒓2)𝑒−𝑖𝜏2𝜓†(𝒓′2) |𝛹𝑁𝐺 ⟩𝜃(𝜏2)− ⟨𝛹𝑁 |𝜓†(𝒓′ )𝑒𝑖𝜏2𝜓(𝒓 ) |𝛹𝑁 ⟩𝜃(−𝜏 )]𝜀𝑖[(𝐸𝑆−𝐸𝐺 )𝑡2+(𝐸𝐺+𝐸𝑆 )|𝜏2|∕2]𝑆 2 2 𝐺 28 = 𝜒𝑆 (𝒓2, 𝒓′2; 𝜏2)𝑒𝑖(𝐸𝑆−𝐸𝐺 )𝑡2 , (A.8b)qs. (A.7a) and (A.7b) become ∑𝑆 𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1)𝜒𝑆 (𝒓2, 𝒓′2; 𝜏2)𝑒−𝑖(𝐸𝑆−𝐸𝐺 )(𝑡1−𝑡2) and ∑𝑆 𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1)𝜒𝑆 (𝒓2, 𝒓′2; 𝜏2)𝑒𝑖(𝐸𝑆−𝐸𝐺 )(𝑡1−𝑡2), respec-tively.Noting that min(𝑡1, 𝑡′1) = 𝑡1 − |𝜏1|∕2 and max(𝑡2, 𝑡′2) = 𝑡2 + |𝜏2|∕2, wefind 𝜃(min(𝑡1, 𝑡′1) − max(𝑡2, 𝑡′2)) = 𝜃(𝑡1 − 𝑡2 − |𝜏1|∕2 − |𝜏2|∕2). Similarly, weave 𝜃(min(𝑡2, 𝑡′2) − max(𝑡1, 𝑡′1)) = 𝜃(𝑡2 − 𝑡1 − |𝜏1|∕2 − |𝜏2|∕2). In the limit2 → −0+, both 𝑡2 and 𝑡′2 become 𝑡2, 𝑡′2 → 𝑡2. Then, the time Fourierransformation of the two-particle Green’s function (A.6) from 𝑡2 to 𝜔can be calculated as∫∞−∞𝑒𝑖(𝐸𝑆−𝐸𝐺 )(𝑡2−𝑡1) 𝜃(𝑡1 − 𝑡2 − |𝜏1|∕2 − |𝜏2|∕2) 𝑒−𝑖𝜔𝑡2𝑑 𝑡2= ∫∞−∞𝑒𝑖(𝐸𝑆−𝐸𝐺 )(𝑡2−𝑡1)−𝑖𝜔𝑡2 𝜃(𝑡1 − 𝑡2 − |𝜏1|∕2) 𝑑 𝑡2= ∫𝑡1−|𝜏1|∕2−∞𝑒𝑖{(𝐸𝑆−𝐸𝐺 )−𝜔−𝑖0+}𝑡2𝑑 𝑡2 𝑒−𝑖(𝐸𝑆−𝐸𝐺 )𝑡1= −𝑖 𝑒𝑖{(𝐸𝑆−𝐸𝐺 )−𝜔}(𝑡1−|𝜏1|∕2)(𝐸𝑆 − 𝐸𝐺) − 𝜔 − 𝑖0+𝑒−𝑖(𝐸𝑆−𝐸𝐺 )𝑡1= 𝑖 𝑒−𝑖𝜔𝑡1𝑒−𝑖(𝐸𝑆−𝐸𝐺 )|𝜏1|∕2𝑒𝑖𝜔|𝜏1|∕2𝜔 − (𝐸𝑆 − 𝐸𝐺) + 𝑖0+, (A.9a)∫∞−∞𝑒𝑖(𝐸𝑆−𝐸𝐺 )(𝑡1−𝑡2) 𝜃(𝑡2 − 𝑡1 − |𝜏1|∕2 − |𝜏2|∕2) 𝑒−𝑖𝜔𝑡2𝑑 𝑡2= ∫∞−∞𝑒𝑖(𝐸𝑆−𝐸𝐺 )(𝑡1−𝑡2)−𝑖𝜔𝑡2 𝜃(𝑡2 − 𝑡1 − |𝜏1|∕2) 𝑑 𝑡2= ∫∞𝑡1+|𝜏1|∕2𝑒𝑖{−(𝐸𝑆−𝐸𝐺 )−𝜔+𝑖0+}𝑡2𝑑 𝑡2 𝑒𝑖(𝐸𝑆−𝐸𝐺 )𝑡1= −𝑖 𝑒𝑖{−(𝐸𝑆−𝐸𝐺 )−𝜔}(𝑡1+|𝜏1|∕2)(𝐸𝑆 − 𝐸𝐺) + 𝜔 − 𝑖0+𝑒𝑖(𝐸𝑆−𝐸𝐺 )𝑡1= −𝑖 𝑒−𝑖𝜔𝑡1𝑒−𝑖(𝐸𝑆−𝐸𝐺 )|𝜏1|∕2𝑒−𝑖𝜔|𝜏1|∕2𝜔 + (𝐸𝑆 − 𝐸𝐺) − 𝑖0+. (A.9b)Thus, we have [51]𝐿(𝒓1, 𝑡1, 𝒓2; 𝒓′1, 𝑡1 − 𝜏1, 𝒓′2, 𝜔)≡ −∫∞−∞𝐿(𝒓1, 𝑡1, 𝒓2, 𝑡2; 𝒓′1, 𝑡1 − 𝜏1, 𝑡2, 𝒓′2, 𝑡2 + 0+)𝑒−𝑖𝜔𝑡2𝑑 𝑡2= − 𝑖𝑒−𝑖𝜔(𝑡1−|𝜏1|∕2)∑𝑆𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1)𝜒𝑆 (𝒓2, 𝒓′2; −0+)𝜔 − (𝐸𝑆 − 𝐸𝐺) + 𝑖0+𝑒−𝑖(𝐸𝑆−𝐸𝐺 )|𝜏1|∕2+ 𝑖𝑒−𝑖𝜔(𝑡1+|𝜏1|∕2)∑𝑆𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1)𝜒𝑆 (𝒓2, 𝒓′2; −0+)𝜔 + (𝐸𝑆 − 𝐸𝐺) − 𝑖0+𝑒−𝑖(𝐸𝑆−𝐸𝐺 )|𝜏1|∕2. (A.10)Eq. (A.10) is still an exact expression without introducing any approx-imation. The first term on the right-hand side comprises the positiveenergy terms associated with exciton creations and the second termencompasses the negative energy terms linked to exciton recombina-tions. In this context, 𝜔 represents the photon energy, which is positiveand can be approximated as being close to some resonant energy 𝛺𝑟 =𝐸𝑟 −𝐸𝐺. As the photon energy approaches this resonant condition, oneof the positive energy terms becomes dominant due to its denominatornearing zero. This allows us to neglect the negative energy terms,i.e., the second term of Eq. (A.10), and rewrite Eq. (A.10) as𝐿(𝒓1, 𝑡1, 𝒓2; 𝒓′1, 𝑡1 − 𝜏1, 𝒓′2, 𝜔) ≃ − 𝑖𝑒−𝑖𝜔𝑡1 1𝜔 −𝛺𝑟 + 𝑖0+×∑𝑆𝜒̄𝑆 (𝒓1, 𝒓′1; 𝜏1)̄̃𝜒𝑆 (𝒓2, 𝒓′2; −0+). (A.11)Here, we put 𝛺𝑟 = 𝐸𝑟 − 𝐸𝐺, and the functions with bar 𝜒̄𝑆 (𝒓1, 𝒓′1; 𝜏1),̄̃𝑆 (𝒓2, 𝒓′2; 𝜏2) are 𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1), 𝜒𝑆 (𝒓2, 𝒓′2; 𝜏2) with the excitation energy𝐸𝑆 − 𝐸𝐺 replaced by 𝜔 in Eqs. (A.8) and (A.12) (below).According to Strinati [51], the particle–hole amplitudes 𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1)and 𝜒𝑆 (𝒓2, 𝒓′2; 𝜏2) are derived from their definition (A.8) by expand-ing the field operators in terms of the single-particle wave functionsK. Ohno et al.cbEi(tiiBbta𝐿abTComputational Condensed Matter 41 (2024) e00977 𝜙𝑐 (𝒓) and 𝜙𝑣(𝒓) (by noting that 𝑒−𝑖𝜏𝑗 and 𝑒𝑖𝜏𝑗 in Eq. (A.8) become𝑒−𝑖(𝐸𝐺+𝜀QP𝑐 )𝜏𝑗 and 𝑒𝑖(𝐸𝐺−𝜀QP𝑣 )𝜏𝑗 , respectively, for both 𝑗 = 1, 2) as𝜒𝑆 (𝒓1, 𝒓′1; 𝜏1) = −𝑒𝑖(𝐸𝑆−𝐸𝐺 )|𝜏1|∕2∑𝑣,𝑐𝐴𝑆 (𝑣, 𝑐)𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)×[𝜃(𝜏1)𝑒−𝑖𝜀QP𝑐 𝜏1 + 𝜃(−𝜏1)𝑒−𝑖𝜀QP𝑣 𝜏1], (A.12a)𝜒𝑆 (𝒓2, 𝒓′2; 𝜏2) = −𝑒𝑖(𝐸𝑆−𝐸𝐺 )|𝜏2|∕2∑𝑣,𝑐𝐴∗𝑆 (𝑣, 𝑐)𝜙𝑣(𝒓2)𝜙∗𝑐 (𝒓′2)×[𝜃(𝜏2)𝑒−𝑖𝜀QP𝑐 𝜏2 + 𝜃(−𝜏2)𝑒−𝑖𝜀QP𝑣 𝜏2], (A.12b)where 𝑣 and 𝑐 stand for valence and conduction QP levels, respec-tively, and 𝑒QP𝑣 and 𝜀QP𝑐 denote the corresponding QP energies. Thisorresponds to the Tamm–Dancoff approximation, neglecting the terms,where 𝐴𝑆 (𝑣, 𝑐) is replaced by 𝐵𝑆 (𝑣, 𝑐), and 𝑐 and 𝑣 are interchanged inEq. (A.12) [57]. Multiplying Eq. (A.12b) with 𝜏2 = −0+ by𝜒𝑟(𝒓′2, 𝒓2; −0+) = −∑𝑣′ ,𝑐′𝐴𝑟(𝑣′, 𝑐′)𝜙𝑐′ (𝒓′2)𝜙∗𝑣′ (𝒓2), (A.13)integrating it with respect to 𝒓2 and 𝒓′2, and using ∫ 𝜙∗𝑐 (𝒓2)𝜙𝑐′ (𝒓2)𝑑𝒓2 =𝛿𝑐 𝑐′ and ∫ 𝜙𝑣(𝒓2)𝜙∗𝑣′ (𝒓2)𝑑𝒓2 = 𝛿𝑣𝑣′ , we have the orthonormality relation∫ 𝜒𝑟(𝒓2, 𝒓′2; −0+)𝜒𝑆 (𝒓2, 𝒓′2; −0+)𝑑𝒓2𝑑𝒓′2 =∑𝑣,𝑐𝐴𝑟(𝑣, 𝑐)𝐴 ∗𝑆 (𝑣, 𝑐) = 𝛿𝑟𝑆 .(A.14)To obtain the matrix form of the BSE, we will finally multiplythe BSE (A.4) by 𝜙∗𝑐 (𝒓1)𝜙𝑣(𝒓′1) and integrate with respect to 𝒓1 and𝒓′1. Keeping this in mind, we start by considering the first term onthe right-hand side of the BSE (A.4), 𝐺(𝑥1, 𝑥′2)𝐺(𝑥2, 𝑥′1), in the limit𝜏1, 𝜏2 → −0+ (𝑡𝑖 = 𝑡′𝑖 − 0+ for 𝑖 = 1, 2). Keeping only the contributionswhere 𝒓1 corresponds to the conduction band (𝑐) and 𝒓′1 corresponds tothe valence band (𝑣), we have𝐺(𝑡1 − 𝑡2)𝐺(𝑡2 − 𝑡1)=∑𝑐 ∫𝑑 𝜔12𝜋𝑒−𝑖𝜔1(𝑡1−𝑡2)𝜔1 − 𝜀QP𝑐 + 𝑖0+∑𝑣 ∫𝑑 𝜔22𝜋𝑒−𝑖𝜔2(𝑡2−𝑡1)𝜔2 − 𝜀QP𝑣 − 𝑖0+=∑𝑐 ,𝑣 ∫𝑑 𝜔1𝑑 𝜔2(2𝜋)2𝑒−𝑖(𝜔1−𝜔2)𝑡1𝑒𝑖(𝜔1−𝜔2)𝑡2(𝜔1 − 𝜀QP𝑐 + 𝑖0+)(𝜔2 − 𝜀QP𝑣 − 𝑖0+). (A.15)Fourier transforming this from 𝑡2 to 𝜔 (note that the 𝑡2 integrationecomes 2𝜋 𝛿(𝜔1 − 𝜔2 − 𝜔), yielding 𝜔2 = 𝜔1 − 𝜔), multiplying it byq. (A.13), and integrating it with respect to 𝒓2 and 𝒓′2, we obtain∫∞−∞𝑑 𝑡2𝑒−𝑖𝜔𝑡2 ∫ 𝐺(𝑥1, 𝑥′2)𝐺(𝑥2, 𝑥′1)𝜒𝑟(𝒓′2, 𝒓2; −0+)𝑑𝒓2𝑑𝒓′2= −𝑒−𝑖𝜔𝑡1∑𝑐 ,𝑣 ∫𝑑 𝜔12𝜋1(𝜔1 − 𝜀QP𝑐 + 𝑖0+)(𝜔1 − 𝜔 − 𝜀QP𝑣 − 𝑖0+)× 𝐴𝑟(𝑣, 𝑐)𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)=∑𝑑 ,𝑐−𝑖𝑒−𝑖𝜔𝑡1𝜔 −(𝜀QP𝑐 − 𝜀QP𝑣)𝐴𝑟(𝑣, 𝑐)𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1). (A.16)We have to further multiply this by −1, because 𝐿(𝜔) in Eq. (A.10)s defined with a minus sign. However, we multiply the entire BSEA.4) by an additional factor −1. Then, this −1 factor is absorbed inthe definition of 𝐿(𝜔) on the left-hand side and in the second term ofhe right-hand side of the BSE. Thus, the explicit multiplication by −1s only necessary for the first term on the right-hand side.Also, for the left-hand side of the BSE, we multiply 𝐿(𝜔) (i.e., Eq.(A.11) with Eq. (A.12a), replacing 𝐸𝑆 − 𝐸𝐺 with 𝜔) by Eq. (A.13) andntegrate over 𝒓2 and 𝒓′2. This gives us, from Eq. (A.14),−(−𝑖𝑒−𝑖𝜔𝑡1𝑒𝑖𝜔|𝜏1|∕2𝜔 −𝛺𝑟 + 𝑖0+)∑𝑣,𝑐𝐴𝑟(𝑣, 𝑐)𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)×[𝜃(𝜏 )𝑒−𝑖𝜀QP𝑐 𝜏1 + 𝜃(−𝜏 )𝑒−𝑖𝜀QP𝑣 𝜏1]. (A.17)1 19 Here, the first negative sign comes from Eq. (A.12a), and the secondnegative sign inside the parenthesis comes from Eq. (A.11). Eq. (A.17)in the limit 𝜏1 = −0+ forms the contribution from the left-hand side ofthe BSE (A.4)𝑖𝑒−𝑖𝜔𝑡1𝜔 −𝛺𝑟 + 𝑖0+∑𝑣,𝑐𝐴𝑟(𝑣, 𝑐)𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1). (A.18)Lastly, we focus on the second term of the right-hand side of theSE (A.4). We multiply it by Eq. (A.13) and integrate with respect to𝒓2 and 𝒓′2. At 𝑡1 = 𝑡′1, the time-dependent part of 𝐺(𝑥1, 𝑥3)𝐺(𝑥′3, 𝑥′1) cane expressed (by keeping only the contributions where 𝒓1 correspondso the conduction band (𝑐) and 𝒓′1 corresponds to the valence band (𝑣))s follows:𝐺(𝑡1 − 𝑡3)𝐺(𝑡′3 − 𝑡1)=∑𝑐 ∫𝑑 𝜔12𝜋𝑒−𝑖𝜔1(𝑡1−𝑡3)𝜔1 − 𝜀QP𝑐 + 𝑖0+∑𝑣 ∫𝑑 𝜔22𝜋𝑒−𝑖𝜔2(𝑡′3−𝑡1)𝜔2 − 𝜀QP𝑣 − 𝑖0+=∑𝑐 ,𝑣 ∫𝑑 𝜔1𝑑 𝜔2(2𝜋)2𝑒−𝑖(𝜔1−𝜔2)𝑡1𝑒𝑖(𝜔1−𝜔2)𝑡3𝑒𝑖(𝜔1+𝜔2)𝜏3∕2(𝜔1 − 𝜀QP𝑐 + 𝑖0+)(𝜔2 − 𝜀QP𝑣 − 𝑖0+). (A.19)The kernel 𝛯(𝑥3, 𝑥4; 𝑥′3, 𝑥′4) depend on 𝑡3 and 𝑡′3 only through 𝜏3 = 𝑡3−𝑡′3,and, due to the 𝛿-functions of this kernel, the correlation function(𝑥′4, 𝑥2; 𝑥4, 𝑥′2) of the form (A.10) depends on 𝑡3 and 𝑡′3 only through𝜏3 = 𝑡3−𝑡′3 except for the phase factor 𝑒−𝑖𝜔𝑡3 . Therefore, after multiplyingll these factors together, we can perform the integral with respectto the intermediate time 𝑡3 = (𝑡3 + 𝑡′3)∕2 (by noting that this integralecomes 2𝜋 𝛿(𝜔1 − 𝜔2 − 𝜔), which leaves only the 𝜔2 = 𝜔1 − 𝜔 term inEq. (A.19), and by performing the 𝜔1 contour integration separately for𝜏3 > 0 and 𝜏3 < 0) as∫∞−∞𝑒−𝑖𝜔𝑡3𝐺(𝑥1, 𝑥3)𝐺(𝑥′3, 𝑥′1)𝑑 𝑡3= 𝑖𝑒−𝑖𝜔𝑡1𝑒𝑖𝜔|𝜏3|∕2∑𝑐 ,𝑣𝜃(𝜏3)𝑒𝑖𝜀QP𝑣 𝜏3 + 𝜃(−𝜏3)𝑒𝑖𝜀QP𝑐 𝜏3𝜔 − (𝜀QP𝑐 − 𝜀QP𝑣 )𝜙𝑐 (𝒓1)𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓′3)𝜙∗𝑣(𝒓′1).(A.20)For the exchange term, the kernel −𝑖𝑣(𝑥3, 𝑥4) acts at the same time𝑡3 = 𝑡4 only, and the left arms of the correlation function are closed.o proceed, we multiply 𝐿(𝑥′4, 𝑥2; 𝑥4, 𝑥′2) by Eq. (A.13). After setting𝑥4 = 𝑥′4, we can rewrite it as Eq. (A.17), replacing both 𝒓1 and 𝒓′1with 𝒓4, while also substituting 𝜏1 with 𝜏3 = −0+. Next, we multiplythis modified expression by Eq. (A.20) with 𝒓3 = 𝒓′3. By integrating theresulting expression with respect to 𝒓3 and 𝒓4, we obtain the exchangeterm𝑖𝑒−𝑖𝜔𝑡1∑𝑐 ,𝑣𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)𝜔 − (𝜀QP𝑐 − 𝜀QP𝑣 )(−𝑖)∫ 𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓3)1|𝒓3 − 𝒓4|× (−1)(−𝑖𝜔 −𝛺𝑟 + 𝑖0+)∑𝑣′ ,𝑐′𝐴𝑟(𝑣′, 𝑐′)𝜙𝑐′ (𝒓4)𝜙∗𝑣′ (𝒓4)𝑑𝒓3𝑑𝒓4. (A.21)For the direct term, we start with Eq. (A.17), in which 𝑥1 and 𝑥′1are replaced by 𝑥3 and 𝑥′3. We multiply this expression by Eq. (A.20)and by the dynamically screened Coulomb interaction (𝑖𝑊 (𝒓3, 𝒓′3; 𝜏3) =𝑖 ∫ 𝑒−𝑖𝜔′𝜏3𝑊 (𝒓3, 𝒓′3;𝜔′) 𝑑 𝜔′∕2𝜋). Next, integrating it with respect to 𝒓3,𝒓′3, and 𝜏3, we obtain the direct term𝑖𝑒−𝑖𝜔𝑡1∑𝑐 ,𝑣𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)𝜔 − (𝜀QP𝑐 − 𝜀QP𝑣 ) ∫𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓′3)∫∞−∞𝑑 𝜏3 𝑒𝑖𝜔|𝜏3|∕2× 𝑖2𝜋 ∫ 𝑑 𝜔′𝑒−𝑖𝜔′𝜏3𝑊 (𝒓3, 𝒓′3;𝜔′) (−1)(−𝑖𝑒𝑖𝜔|𝜏3|∕2𝜔 −𝛺𝑟 + 𝑖0+)×∑𝑣′ ,𝑐′𝐴𝑟(𝑣′, 𝑐′)𝜙𝑐′ (𝒓3)𝜙∗𝑣′ (𝒓′3)×[𝜃(𝜏3)𝑒𝑖(𝜀QP𝑣 −𝜀QP𝑐′)𝜏3 + 𝜃(−𝜏3)𝑒𝑖(𝜀QP𝑐 −𝜀QP𝑣′)𝜏3]𝑑𝒓3𝑑𝒓′3= 𝑖𝑒−𝑖𝜔𝑡1∑ 𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)QP QP𝑖∫𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓′3)+𝑐 ,𝑣 𝜔 − (𝜀𝑐 − 𝜀𝑣 ) 𝜔 −𝛺𝑟 + 𝑖0K. Ohno et al.ameTbvectfa4gisEComputational Condensed Matter 41 (2024) e00977 × 𝑖2𝜋 ∫ 𝑑 𝜔′𝑊 (𝒓3, 𝒓′3;𝜔′)∑𝑣′ ,𝑐′𝐴𝑟(𝑣′, 𝑐′)𝜙𝑐′ (𝒓3)𝜙∗𝑣′ (𝒓′3)𝑑𝒓3𝑑𝒓′3×[∫∞0𝑑 𝜏3𝑒𝑖(𝜔−𝜔′+𝜀QP𝑣 −𝜀QP𝑐′+𝑖0+)𝜏3 + ∫0−∞𝑑 𝜏3𝑒𝑖(−𝜔−𝜔′+𝜀QP𝑐 −𝜀QP𝑣′−𝑖0+)𝜏3]= 𝑖𝑒−𝑖𝜔𝑡1∑𝑐 ,𝑣𝜙𝑐 (𝒓1)𝜙∗𝑣(𝒓′1)𝜔 − (𝜀QP𝑐 − 𝜀QP𝑣 )𝑖∫𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓′3)𝜔 −𝛺𝑟 + 𝑖0+× 𝑖2𝜋 ∫ 𝑑 𝜔′𝑊 (𝒓3, 𝒓′3;𝜔′)∑𝑣′ ,𝑐′𝐴𝑟(𝑣′, 𝑐′)𝜙𝑐′ (𝒓3)𝜙∗𝑣′ (𝒓′3)𝑑𝒓3𝑑𝒓′3.×⎡⎢⎢⎣𝑖𝜔 − 𝜔′ − (𝜀QP𝑐′ − 𝜀QP𝑣 ) + 𝑖0++ 𝑖𝜔 + 𝜔′ − (𝜀QP𝑐 − 𝜀QP𝑣′ ) + 𝑖0+⎤⎥⎥⎦(A.22)In the last two denominators, 𝜔 can be replaced by the resulting PAE,𝛺𝑟.To derive the matrix eigenvalue equation for the BSE, we combinell these terms on both sides of the BSE (A.4). Multiplying both sides by𝜙∗𝑐 (𝒓1)𝜙𝑣(𝒓′1), and integrating them with respect to 𝒓1 and 𝒓′1, we obtainthe final equation(𝜀QP𝑐 − 𝜀QP𝑣 )𝐴𝑟(𝑣, 𝑐) = 𝛺𝑟𝐴𝑟(𝑣, 𝑐) −∑𝑣′ ,𝑐′𝑋𝑣𝑐;𝑣′𝑐′𝐴𝑟(𝑣′, 𝑐′)+∑𝑣′ ,𝑐′𝑊𝑐 𝑣;𝑐′𝑑′𝐴𝑟(𝑣′, 𝑐′), (A.23)where the exchange term is given by Eq. (12), i.e.,𝑋𝑣𝑐;𝑣′𝑐′ = ∫ 𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓3)1|𝒓3 − 𝒓4|𝜙𝑐′ (𝒓4)𝜙∗𝑣′ (𝒓4)𝑑𝒓3𝑑𝒓4, (A.24)and the direct term is given by Eq. (9) (Note that there is a typo in thepioneering papers [71]), i.e.,𝑊𝑣𝑐;𝑣′𝑐′ = 𝑖∫𝑑 𝜔′2𝜋𝜙∗𝑐 (𝒓3)𝜙𝑣(𝒓′3)𝑊 (𝒓3, 𝒓′3;𝜔′)[1𝛺𝑟 − 𝜔′ − (𝜀QP𝑐′ − 𝜀QP𝑣 ) + 𝑖0++ 1𝛺𝑟 + 𝜔′ − (𝜀QP𝑐 − 𝜀QP𝑣′ ) + 𝑖0+]𝜙𝑐′ (𝒓3)𝜙∗𝑣′ (𝒓′3)𝑑𝒓3𝑑𝒓′3. (A.25)To avoid the 𝜔′ integration in Eq. (A.25), we can utilize the GPPodel [71] or other plasmon pole models [81,82], as well as staticapproximations [54]. However, prior to applying these approximations,it is better to separate the bare Coulomb term from Eq. (A.25) asxpressed in Eq. (10). This is because the bare Coulomb term must betreated with high precision, especially in all-electron calculations. Thedynamically screened Coulomb interaction 𝑊 is Fourier transformedinto the reciprocal space as𝑊 (𝒓, 𝒓′;𝜔′) =∑𝑮 𝑮′𝑒𝑖𝑮⋅𝒓𝑊𝑮 𝑮′ (𝜔′)𝑒−𝑖𝑮⋅𝒓′ . (A.26)Hereafter, we always assume that the momentum transfer is 𝒒 = 0.hen, the symbolic relation 𝑊 = 𝜀−1𝑣 can be written explicitly as𝑊𝑮 𝑮′ (𝜔′) = 𝜀−1𝑮 𝑮′ (𝜔′)𝑣(𝑮′), (A.27)where 𝑣(𝑮) = 4𝜋∕𝛺|𝑮|2 is the Fourier transform of the bare Coulombinteraction with 𝛺 being the volume of the supercell. Under the GPPmodel at 𝒒 = 0, we separate 𝜀−1𝑮 𝑮′ (𝜔′) as [44]𝜀−1𝑮 𝑮′ (𝜔′) = 𝛿𝑮 𝑮′ +𝛺2𝑮 𝑮′2𝜔̃𝑮 𝑮′[1𝜔′ − 𝜔̃𝑮 𝑮′ + 𝑖0+− 1𝜔′ + 𝜔̃𝑮 𝑮′ − 𝑖0+]. (A.28)Here, the first term 𝛿𝑮 𝑮′ is the Fourier transformation of 1, whichrepresents the bare Coulomb term 𝑣 in 𝑊 = 𝜀−1𝑣 = 𝑣 + (𝜀−1 − 1)𝑣,while the second term corresponds to (𝜀−1− 1)𝑣. The symbols 𝛺2𝑮 𝑮′ and𝜔̃𝑮 𝑮′ are defined, respectively, as [44]𝛺2𝑮 𝑮′ = 4𝜋𝑮 ⋅𝑮′|𝑮|2𝜌(𝑮 −𝑮′), (A.29)𝜔̃2𝑮 𝑮′ =𝛺2𝑮 𝑮′−1. (A.30)𝛿𝑮 𝑮′ − 𝜀𝑮 𝑮′ (0)10 Then, in the 𝜔′ integration of Eq. (A.25), the first term becomes theare Coulomb term 𝑈𝑣𝑐;𝑣′𝑐′ given by Eq. (11),𝑈𝑣𝑐;𝑣′𝑐′ = ∫ 𝜙∗𝑐 (𝒓)𝜙𝑐′ (𝒓)1|𝒓 − 𝒓′|𝜙𝑣(𝒓′)𝜙∗𝑣′ (𝒓′) 𝑑𝒓𝑑𝒓′, (A.31)which is expressed in the reciprocal space as𝑈𝑣𝑐;𝑣′𝑐′ =′∑𝑮⟨ 𝑐 | 𝑒𝑖𝑮⋅𝒓| 𝑐′ ⟩⟨ 𝑣′ | 𝑒−𝑖𝑮⋅𝒓′| 𝑣 ⟩𝑣(𝑮), (A.32)where ∑′𝑮 means that 𝑮 = 0 is excluded from the summation in thediagonal terms 𝑐 = 𝑐′ and 𝑣 = 𝑣′ for the reason discussed aroundEq. (13) in Section 4. In the off-diagonal terms 𝑐 ≠ 𝑐′ or 𝑣 ≠ 𝑣′, the𝑮 = 0 term automatically vanishes because ⟨ 𝑐 | 𝑐′ ⟩ = 0 or ⟨ 𝑣 | 𝑣′ ⟩ = 0.The second term becomes the remaining correlation term 𝑆𝑣𝑐;𝑣′𝑐′ , whichis given within the GPP model by [71]𝑆𝑣𝑐;𝑣′𝑐′ =∑𝑮∑𝑮′⟨ 𝑐 | 𝑒𝑖𝑮⋅𝒓| 𝑐′ ⟩⟨ 𝑣′ | 𝑒−𝑖𝑮′⋅𝒓′| 𝑣 ⟩×𝛺2𝑮 𝑮′2𝜔̃𝑮 𝑮′[1𝛺𝑟 − (𝜀QP𝑐′ − 𝜀QP𝑣 ) − 𝜔̃𝑮 𝑮′+ 1𝛺𝑟 − (𝜀QP𝑐 − 𝜀QP𝑣′ ) − 𝜔̃𝑮 𝑮′]𝑣(𝑮′). (A.33)If 1∕|𝑮|2 and 1∕|𝑮′|2 in 𝛺2𝑮 𝑮′ and 𝑣(𝑮′) are replaced with some nonzeroalues in the limit 𝑮,𝑮′ → 0 [78], the 𝑮 = 0 and/or 𝑮′ = 0contributions to 𝑆𝑣𝑐;𝑣′𝑐′ automatically vanish due to the form of 𝛺2𝑮 𝑮′in Eq. (A.29). The same situation happens also when the GPP modelis directly used in 𝑊𝑣𝑐;𝑣′𝑐′ . In the reciprocal space, the exchange term(A.24) is expressed as𝑋𝑣𝑐;𝑣′𝑐′ =∑𝑮⟨ 𝑐 | 𝑒𝑖𝑮⋅𝒓| 𝑣 ⟩⟨ 𝑣′ | 𝑒−𝑖𝑮⋅𝒓′| 𝑐′ ⟩𝑣(𝑮). (A.34)Note that, here again, the 𝑮 = 0 term in the summation vanishesbecause ⟨ 𝑐 | 𝑣 ⟩ = ⟨ 𝑣′ | 𝑐′ ⟩ = 0.For the evaluation of 𝑈𝑣𝑐;𝑣′𝑐′ and 𝑋𝑣𝑐;𝑣′𝑐′ terms as well as the Fockxchange term in the 𝐺 𝑊 level, we excluded the fully on-site AO-onlyontributions, which are separately calculated in a single integral alonghe radial direction, as accurately as possible, using analytic integralorms of the Poisson equation for one-center problem.Up to now, we have ignored electron spins. However, it is requiredto consider the spin state of each level, in order to distinguish thesinglet and triplet excited (exciton) states. Indeed, the direct 𝑊𝑣𝑐;𝑣′𝑐′term exists only when the spin states of 𝑣 and 𝑣′ are the same (↑↑ or↓↓) and the spin states of 𝑐 and 𝑐′ are the same, while the exchange𝑋𝑣𝑐;𝑣′𝑐′ term exists only when the spin states of 𝑣 and 𝑐 are the samend the spin states of 𝑣′ and 𝑐′ are the same. Therefore, considering× 4 matrix consisting of ↑↑, ↑↓, ↓↑, and ↓↓ rows and columns, thelobal matrix form of the BSE becomes (note that each matrix elements also the matrix)⎡⎢⎢⎢⎢⎣𝐷 +𝑋 0 0 𝑋0 𝐷 0 00 0 𝐷 0𝑋 0 0 𝐷 +𝑋⎤⎥⎥⎥⎥⎦⎡⎢⎢⎢⎢⎣𝐴↑↑𝑟𝐴↑↓𝑟𝐴↓↑𝑟𝐴↓↓𝑟⎤⎥⎥⎥⎥⎦= 𝛺𝑟⎡⎢⎢⎢⎢⎣𝐴↑↑𝑟𝐴↑↓𝑟𝐴↓↑𝑟𝐴↓↓𝑟⎤⎥⎥⎥⎥⎦, (A.35)where the matrix elements of 𝐷 (𝐷𝑣𝑐;𝑣′𝑐′ ) are defined in Eq. (8). Thisglobal matrix eigenvalue equation can readily be solved to yield onesinglet solution (5) and three degenerate triplet solutions (6). Eacholution is, of course, the matrix eigenvalue equation, either the singletq. (5) or the (triply degenerate) triplet Eq. (6). The full 8 × 8 matrixform including not only 𝐴𝜎 𝜎′𝑟 but also 𝐵𝜎 𝜎′𝑟 beyond the Tamm–Dancoffapproximation can be found in Ref. [58].Data availabilityData will be made available on request.K. Ohno et al. Computational Condensed Matter 41 (2024) e00977 References[1] A. Fujishima, K. Honda, Electrochemical photolysis of water at a semiconductorelectrode, Nature 238 (1972) 37–38, http://dx.doi.org/10.1038/238037a0.[2] H.P. Maruska, A.K. Ghosh, Photocatalytic decomposition of water at semiconduc-tor electrodes, Sol. Energy 20 (1978) 443–458, http://dx.doi.org/10.1016/0038-092X(78)90061-0.[3] Hisashi Harada, T. Ueda, Photocatalytic activity of ultra-fine rutile in methanol-water solution and dependence of activity on particle size, Chem. Phys. Lett. 106(1984) 229–231, http://dx.doi.org/10.1016/0009-2614(84)80231-6.[4] A. Fujishima, T.N. Rao, D.A. Tryk, Titanium dioxide photocatalysis, J. Pho-tochem. Photobiol. C Photochem. Rev. 1 (2000) 1–21, http://dx.doi.org/10.1016/S1389-5567(00)00002-2.[5] R.W. Matthews, Solar-electric water purification using photocatalytic oxidationwith TiO2 as a stationary phase, Sol. Energy 38 (1987) 405–413, http://dx.doi.org/10.1016/0038-092X(87)90021-1.[6] R. Wang, K. Hashimoto, A. Fujishima, M. Chikuni, E. Kojima, A. Kitamura,M. Shimohigoshi, T. Watanabe, Photogeneration of highly amphiphilic TiO2surfaces, Adv. Mater. 10 (1998) 135–138, http://dx.doi.org/10.1002/(SICI)1521-4095(199801)10:2<135::AID-ADMA135>3.0.CO;2-M.[7] C. Longo, A.F. Nogueira, M.-A. De Paoli, H. Cachet, Solid-state and flexible dye-sensitized TiO2 solar cells: a study by electrochemical impedance spectroscopy,J. Phys. Chem. B 106 (2002) 5925–5930, http://dx.doi.org/10.1021/jp014456u.[8] K. Sunada, T. Watanabe, K. Hashimoto, Studies on photokilling of bacteria onTiO2 thin film, J. Photochem. Photobiol. A 156 (2003) 227–233, http://dx.doi.org/10.1016/S1010-6030(02)00434-3.[9] R. Carbone, I. Marangi, A. Zanardi, L. Giorgetti, E. Chierici, G. Berlanda, A.Podest/‘a, F. Fiorentini, G. Bongiorno, P. Piseri, P.G. Pelicci, P. Milani, Bio-compatibility of cluster-assembled nanostructured TiO2 with primary and cancercells, Biomater. 27 (2006) 3221–3229, http://dx.doi.org/10.1016/j.biomaterials.2006.01.056.[10] J.M.D. Coey, S.A. Chambers, Oxide dilute magnetic semiconductors—fact orfiction? MRS Bull. 33 (2008) 1053–1058, http://dx.doi.org/10.1557/mrs2008.225.[11] H. Saadaoui, X. Luo, Z. Salman, X.Y. Cui, N.N. Bao, P. Bao, R.K. Zheng, L.T.Tseng, Y.H. Du, T. Prokscha, A. Suter, T. Liu, Y.R. Wang, S. Li, J. Ding, S.P.Ringer, E. Morenzoni, J.B. Yi, Intrinsic ferromagnetism in the diluted magneticsemiconductor Co:TiO2, Phys. Rev. Lett. 117 (2016) 227202, http://dx.doi.org/10.1103/PhysRevLett.117.227202.[12] B. Prajapati, S. Kumar, M. Kumar, S. Chatterjeeb, A.K. Ghosh, Investigation of thephysical properties of Fe:TiO2-diluted magnetic semiconductor nanoparticles, J.Mater. Chem. C 5 (2017) 4257–4267, http://dx.doi.org/10.1039/C7TC00233E.[13] N. Fajariah, W.A.E. Prabowo, F. Fathurrahman, A. Melati, H.K. Dipojono, Theinvestigation of electronic structure of transition metal doped TiO2 for dilutedmagnetic semiconductor applications: A first principle study, Procedia Eng. 170(2017) 141–147, http://dx.doi.org/10.1016/j.proeng.2017.03.032.[14] R.C. da Silva, E. Alves, M.M. Cruz, Conductivity behaviour of Cr implanted TiO2,Nucl. Instrum. Methods Phys. Res. B 191 (2002) 158–162, http://dx.doi.org/10.1016/S0168-583X(02)00541-4.[15] E. Dy, R. Hui, J. Zhang, Z.-S. Liu, Z. Shi, Electronic conductivity and stabilityof doped titania (Ti1−𝑥M𝑥O2, M = Nb, Ru, and Ta) – A density functionaltheory-based comparison, J. Phys. Chem. C 114 (2010) 13162–13167, http://dx.doi.org/10.1021/jp100826g.[16] P. Triggs, F. Lévy, Optical and electronic properties of ruthenium-doped TiO2,Phys. Status Solidi (b) 129 (1985) 363–374, http://dx.doi.org/10.1002/pssb.2221290136.[17] L. Kernazhitsky, V. Shymanovska, T. Gavrilko, V. Naumov, V. Kshnyakin,T. Khalyavka, A comparative study of optical absorption and photocatalyticproperties of nanocrystalline single-phase anatase and rutile TiO2 doped withtransition metal cations, J. Solid State Chem. 198 (2013) 511–519, http://dx.doi.org/10.1016/j.jssc.2012.11.015.[18] K.A. Rahman, N. Sharma, A.J. Atanacio, T. Bak, E.D. Wachsman, M. Moffitt, J.Nowotny, Chromium segregation in Cr-doped TiO2 (rutile): impact of oxygenactivity, Ionics 25 (2019) 3363–3372, http://dx.doi.org/10.1007/s11581-018-2828-4.[19] C. Gionco, A. Battiato, E. Vittone, M.C. Paganini, E. Giamello, Structural andspectroscopic properties of high temperature prepared ZrO2-TiO2 mixed oxides,J. Solid State Chem. 201 (2013) 222–228, http://dx.doi.org/10.1016/j.jssc.2013.02.040.[20] M.A.R. Sarker, Optical properties of Al- and Zr-doped rutile single crystals grownby tilting-mirror-type floating zone method and study of structure–propertyrelationships by first principle calculations, J. Inorg. Chem. 2014 (2014) 274165,http://dx.doi.org/10.1155/2014/274165.[21] F. Ayedun, P.O. Adebambo, B.I. Adetunji, V.C. Ozebo, J.A. Oguntuase, G.A.Adebayo, Increased malleability in tetragonal Zr𝑥Ti1−𝑥O2 ternary alloys: First-principles approach, Z. Nat.forsch. A 72 (2017) 567–572, http://dx.doi.org/10.1515/zna-2017-0036.[22] A. Iwaszuk, M. Nolan, Electronic structure and reactivity of Ce- and Zr-dopedTiO2: Assessing the reliability of density functional theory approaches, J. Phys.Chem. C 115 (2011) 12995–13007, http://dx.doi.org/10.1021/jp203112p.11 [23] H. Maleki-Ghaleh, M.S. Shakeri, Z. Dargahi, M. Kavanlouei, H.K. Garabagh,E. Moradpur-Tari, A. Yourdkhani, A. Fallah, A. Zarrabi, B. Koc, M.H. Sia-dati, Characterization and optical properties of mechanochemically synthesizedmolybdenum-doped rutile nanoparticles and their electronic structure studiesby density functional theory, Mater. Today Chem. 24 (2022) 100820, http://dx.doi.org/10.1016/j.mtchem.2022.100820.[24] X. Yu, C. Li, Y. Ling, T.-A. Tang, Q. Wu, J. Kong, First principles calculations ofelectronic and optical properties of Mo-doped rutile TiO2, J. Alloys Compd. 507(2010) 33–37, http://dx.doi.org/10.1016/j.jallcom.2010.07.195.[25] X. Lu, T. Zhao, X. Gao, J. Ren, X. Yan, P. La, Investigation of Mo-, Pt-, andRh-doped rutile TiO2 based on first-principles calculations, AIP Adv. 8 (2018)075014, http://dx.doi.org/10.1063/1.5038776.[26] A. Soussi, A. Ait hssi, L. Boulkaddat, M. Boujnah, K. Abouabassi, R. Haounati,A. Asbayou, A. Elfanaoui, R. Markazi, A. Ihlal, K. Bouabid, N. El Biaze, Firstprinciple study of electronic, optical and electrical properties of Mo doped TiO2,Comput. Condens. Matter 29 (2021) e00606, http://dx.doi.org/10.1016/j.cocom.2021.e00606.[27] L.-j. Bai, G. Kou, Z.-y. Gong, Z.-m. Zhao, Effect of Zn and Ti mole ratioon microstructure and photocatalytic properties of magnetron sputtered TiO2-ZnO heterogeneous composite film, Trans. Nonferr. Met. Soc. China 23 (2013)3643–3649, http://dx.doi.org/10.1016/S1003-6326(13)62912-X.[28] N.A. Erfan, M.S. Mahmoud, H.Y. Kim, N.A.M. Barakat, Synergistic dopingwith Ag, CdO, and ZnO to overcome electron–hole recombination in TiO2photocatalysis for effective water photo splitting reaction, Front. Chem. 11(2023) 1301172, http://dx.doi.org/10.3389/fchem.2023.1301172.[29] G.H. Mohammed, A.M. Savore, M.H. Shinen, K.A. Adem, Structural and opticalproperties of CdO doped TiO2 thin films prepared by pulsed laser deposition,Eng. Tech. J. 33 (B) (2015) 919–931, http://dx.doi.org/10.30684/etj.33.5B.16.[30] K. Sahbeni, M. Jlassi, S. Khamlich, M. Kandyla, M. Kompitsas, W. Dimassi,Effect of CdO ratios on the structural and optical properties of CdO-TiO2nanocomposite thin films, J. Mater. Sci.: Mater. Electron. 31 (2020) 3387–3396,http://dx.doi.org/10.1007/s10854-020-02887-w.[31] K. Song, X. Han, G. Shao, Electronic properties of rutile TiO2 doped with 4dtransition metals: First-principles study, J. Alloys Compd. 551 (2013) 118–124,http://dx.doi.org/10.1016/j.jallcom.2012.09.077.[32] M. Saini, M. Kumar, T. Som, Ab initio study of 3d transition metal-doping effectsin rutile-TiO2: Role of bandgap tunability in conductivity behaviour, Appl. Surf.Sci. 418 (2017) 302–307, http://dx.doi.org/10.1016/j.apsusc.2017.01.262.[33] K. Sato, H. Akai, T. Minamisono, Ab initio calculations of electric field gradientsfor transition metal impurities in TiO2, Z. Nat.forsch. 53a (1998) 396–403,http://dx.doi.org/10.1515/zna-1998-6-720.[34] L.A. Errico, G. Fabricius, Mario Rentería, Metal impurities in an oxide: Ab initiostudy of electronic and structural properties of Cd in rutile TiO2, Phys. Rev. B67 (2003) 144104, http://dx.doi.org/10.1103/PhysRevB.67.144104.[35] L. Hedin, New method for calculating the one-particle Green’s function withapplication to the electron-gas problem, Phys. Rev. 139 (1965) A796–A823,http://dx.doi.org/10.1103/PhysRev.139.A796.[36] G. Strinati, Application of the Green’s functions method to the study of theoptical properties of semiconductors, Nuovo Cimento 11 (1988) 1–86, http://dx.doi.org/10.1007/BF02725962.[37] A.J. Layzer, Properties of the one-particle Green s function for nonuniform many-Fermion systems, Phys. Rev. 129 (1963) 897–907, http://dx.doi.org/10.1103/PhysRev.129.897.[38] F. Hüser, T. Olsen, K.S. Thygesen, Quasiparticle GW calculations for solids,molecules, and two-dimensional materials, Phys. Rev. B 87 (2013) 235132,http://dx.doi.org/10.1103/PhysRevB.87.235132.[39] T. Nakashima, H. Raebiger, K. Ohno, Normalization of exact quasiparticle wavefunctions in the Green’s function method guaranteed by the Ward identity, Phys.Rev. B 104 (2021) L201116, http://dx.doi.org/10.1103/PhysRevB.104.L201116.[40] K. Ohno, S. Ono, T. Isobe, A simple derivation of the exact quasiparticle theoryand its extension to arbitrary initial excited eigenstates, J. Chem. Phys. 146(2017) 084108, http://dx.doi.org/10.1063/1.4976553.[41] K. Ohno, K. Esfarjani, Y. Kawazoe, Computational Materials Science: FromAb Initio to Monte Carlo Methods, second ed., Springer-Verlag, Berlin, ISBN:9783662565407, 2018.[42] G. Strinati, H.J. Mattausch, W. Hanke, Dynamical aspects of correlation correc-tions in a covalent crystal, Phys. Rev. B 25 (1982) 2867–2888, http://dx.doi.org/10.1103/PhysRevB.25.2867.[43] M.S. Hybertsen, S.G. Louie, First-principles theory of quasiparticles: Calculationof band gaps in semiconductors and insulators, Phys. Rev. Lett. 55 (1985)1418–1421, http://dx.doi.org/10.1103/PhysRevLett.55.1418.[44] M.S. Hybertsen, S.G. Louie, Electron correlation in semiconductors and insulators:Band gaps and quasiparticle energies, Phys. Rev. B 34 (1986) 5390–5413,http://dx.doi.org/10.1103/PhysRevB.34.5390.[45] R.W. Godby, M. Schlüter, L.J. Sham, Accurate exchange–correlation potentialfor silicon and its discontinuity on addition of an electron, Phys. Rev. Lett. 56(1986) 2415–2418, http://dx.doi.org/10.1103/PhysRevLett.56.2415.[46] R.W. Godby, M. Schlüter, L.J. Sham, Self-energy operatorsand exchange–correlation potentials in semiconductors, Phys. Rev. B 37 (1988) 10159–10179,http://dx.doi.org/10.1103/PhysRevB.37.10159.http://dx.doi.org/10.1038/238037a0http://dx.doi.org/10.1016/0038-092X(78)90061-0http://dx.doi.org/10.1016/0038-092X(78)90061-0http://dx.doi.org/10.1016/0038-092X(78)90061-0http://dx.doi.org/10.1016/0009-2614(84)80231-6http://dx.doi.org/10.1016/S1389-5567(00)00002-2http://dx.doi.org/10.1016/S1389-5567(00)00002-2http://dx.doi.org/10.1016/S1389-5567(00)00002-2http://dx.doi.org/10.1016/0038-092X(87)90021-1http://dx.doi.org/10.1016/0038-092X(87)90021-1http://dx.doi.org/10.1016/0038-092X(87)90021-1http://dx.doi.org/10.1002/(SICI)1521-4095(199801)10:2<135::AID-ADMA135>3.0.CO;2-Mhttp://dx.doi.org/10.1002/(SICI)1521-4095(199801)10:2<135::AID-ADMA135>3.0.CO;2-Mhttp://dx.doi.org/10.1002/(SICI)1521-4095(199801)10:2<135::AID-ADMA135>3.0.CO;2-Mhttp://dx.doi.org/10.1021/jp014456uhttp://dx.doi.org/10.1016/S1010-6030(02)00434-3http://dx.doi.org/10.1016/S1010-6030(02)00434-3http://dx.doi.org/10.1016/S1010-6030(02)00434-3http://dx.doi.org/10.1016/j.biomaterials.2006.01.056http://dx.doi.org/10.1016/j.biomaterials.2006.01.056http://dx.doi.org/10.1016/j.biomaterials.2006.01.056http://dx.doi.org/10.1557/mrs2008.225http://dx.doi.org/10.1557/mrs2008.225http://dx.doi.org/10.1557/mrs2008.225http://dx.doi.org/10.1103/PhysRevLett.117.227202http://dx.doi.org/10.1103/PhysRevLett.117.227202http://dx.doi.org/10.1103/PhysRevLett.117.227202http://dx.doi.org/10.1039/C7TC00233Ehttp://dx.doi.org/10.1016/j.proeng.2017.03.032http://dx.doi.org/10.1016/S0168-583X(02)00541-4http://dx.doi.org/10.1016/S0168-583X(02)00541-4http://dx.doi.org/10.1016/S0168-583X(02)00541-4http://dx.doi.org/10.1021/jp100826ghttp://dx.doi.org/10.1021/jp100826ghttp://dx.doi.org/10.1021/jp100826ghttp://dx.doi.org/10.1002/pssb.2221290136http://dx.doi.org/10.1002/pssb.2221290136http://dx.doi.org/10.1002/pssb.2221290136http://dx.doi.org/10.1016/j.jssc.2012.11.015http://dx.doi.org/10.1016/j.jssc.2012.11.015http://dx.doi.org/10.1016/j.jssc.2012.11.015http://dx.doi.org/10.1007/s11581-018-2828-4http://dx.doi.org/10.1007/s11581-018-2828-4http://dx.doi.org/10.1007/s11581-018-2828-4http://dx.doi.org/10.1016/j.jssc.2013.02.040http://dx.doi.org/10.1016/j.jssc.2013.02.040http://dx.doi.org/10.1016/j.jssc.2013.02.040http://dx.doi.org/10.1155/2014/274165http://dx.doi.org/10.1515/zna-2017-0036http://dx.doi.org/10.1515/zna-2017-0036http://dx.doi.org/10.1515/zna-2017-0036http://dx.doi.org/10.1021/jp203112phttp://dx.doi.org/10.1016/j.mtchem.2022.100820http://dx.doi.org/10.1016/j.mtchem.2022.100820http://dx.doi.org/10.1016/j.mtchem.2022.100820http://dx.doi.org/10.1016/j.jallcom.2010.07.195http://dx.doi.org/10.1063/1.5038776http://dx.doi.org/10.1016/j.cocom.2021.e00606http://dx.doi.org/10.1016/j.cocom.2021.e00606http://dx.doi.org/10.1016/j.cocom.2021.e00606http://dx.doi.org/10.1016/S1003-6326(13)62912-Xhttp://dx.doi.org/10.3389/fchem.2023.1301172http://dx.doi.org/10.30684/etj.33.5B.16http://dx.doi.org/10.1007/s10854-020-02887-whttp://dx.doi.org/10.1016/j.jallcom.2012.09.077http://dx.doi.org/10.1016/j.apsusc.2017.01.262http://dx.doi.org/10.1515/zna-1998-6-720http://dx.doi.org/10.1103/PhysRevB.67.144104http://dx.doi.org/10.1103/PhysRev.139.A796http://dx.doi.org/10.1007/BF02725962http://dx.doi.org/10.1007/BF02725962http://dx.doi.org/10.1007/BF02725962http://dx.doi.org/10.1103/PhysRev.129.897http://dx.doi.org/10.1103/PhysRev.129.897http://dx.doi.org/10.1103/PhysRev.129.897http://dx.doi.org/10.1103/PhysRevB.87.235132http://dx.doi.org/10.1103/PhysRevB.104.L201116http://dx.doi.org/10.1063/1.4976553http://refhub.elsevier.com/S2352-2143(24)00099-6/sb41http://refhub.elsevier.com/S2352-2143(24)00099-6/sb41http://refhub.elsevier.com/S2352-2143(24)00099-6/sb41http://refhub.elsevier.com/S2352-2143(24)00099-6/sb41http://refhub.elsevier.com/S2352-2143(24)00099-6/sb41http://dx.doi.org/10.1103/PhysRevB.25.2867http://dx.doi.org/10.1103/PhysRevB.25.2867http://dx.doi.org/10.1103/PhysRevB.25.2867http://dx.doi.org/10.1103/PhysRevLett.55.1418http://dx.doi.org/10.1103/PhysRevB.34.5390http://dx.doi.org/10.1103/PhysRevLett.56.2415http://dx.doi.org/10.1103/PhysRevB.37.10159K. Ohno et al. Computational Condensed Matter 41 (2024) e00977 [47] M. Rohlfing, P. Krüger, J. Pollmann, Quasiparticle band-structure calculationsfor C, Si, Ge, GaAs, and SiC using Gaussian-orbital basis sets, Phys. Rev. B 48(1993) 17791–17803, http://dx.doi.org/10.1103/PhysRevB.48.17791.[48] O. Zakharov, A. Rubio, X. Blase, M.L. Cohen, S.G. Louie, Quasiparticle bandstructures of six II-VI compounds: ZnS, ZnSe, ZnTe, CdS, CdSe, and CdTe, Phys.Rev. B 50 (1994) 10780–10787, http://dx.doi.org/10.1103/PhysRevB.50.10780.[49] W. Kohn, L.J. Sham, Self-consistent equations including exchange and correla-tion effects, Phys. Rev. 140 (1965) A1133–A1138, http://dx.doi.org/10.1103/PhysRev.140.A1133.[50] P. Hohenberg, W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136 (1964)B864–B871, http://dx.doi.org/10.1103/PhysRev.136.B864.[51] G. Strinati, Effects of dynamical screening on resonances at inner-shell thresholdsin semiconductors, Phys. Rev. B 29 (1984) 5718–5726, http://dx.doi.org/10.1103/PhysRevB.29.5718.[52] G. Onida, L. Reining, R.W. Godby, R. Del Sole, W. Andreoni, Ab Initio calcula-tions of the quasiparticle and absorption spectra of clusters: The sodium tetramer,Phys. Rev. Lett. 75 (1995) 818–821, http://dx.doi.org/10.1103/PhysRevLett.75.818.[53] M. Rohlfing, S.G. Louie, Excitonic effects and the optical absorption spectrum ofhydrogenated Si clusters, Phys. Rev. Lett. 80 (1998) 3320–3323, http://dx.doi.org/10.1103/PhysRevLett.80.3320.[54] S. Albrecht, L. Reining, R. Del Sole, G. Onida, Ab initio calculation of excitoniceffects in the optical spectra of semiconductors, Phys. Rev. Lett. 80 (1998)4510–4513, http://dx.doi.org/10.1103/PhysRevLett.80.4510.[55] M. Rohlfing, S.G. Louie, Electron-hole excitations in semiconductors and in-sulators, Phys. Rev. Lett. 81 (1998) 2312–2315, http://dx.doi.org/10.1103/PhysRevLett.81.2312.[56] L.X. Benedict, E.L. Shirley, Ab initio calculation of 𝜀2(𝜔) including the electron–hole interaction: Application to GaN and CaF2, Phys. Rev. B 59 (1999)5441–5451, http://dx.doi.org/10.1103/PhysRevB.59.5441.[57] M. Rohlfing, S.G. Louie, Electron–hole excitations and optical spectra fromfirst principles, Phys. Rev. B 62 (2000) 4927–4944, http://dx.doi.org/10.1103/PhysRevB.62.4927.[58] K. Ohno, Optical properties of alkali-earth atoms and Na2 calculated by GWand Bethe–Salpeter equations, Sci. Technol. Adv. Mater. 5 (2004) 603–607,http://dx.doi.org/10.1016/j.stam.2004.02.018.[59] Y. Noguchi, K. Ohno, All-electron first-principles 𝐺 𝑊 + Bethe–Salpeter calcula-tion for optical absorption spectra of sodium clusters, Phys. Rev. A 81 (2010)045201, http://dx.doi.org/10.1103/PhysRevA.81.045201.[60] S. Yamada, Y. Noguchi, K. Ishii, D. Hirose, O. Sugino, K. Ohno, Developmentof the Bethe–Salpeter method considering second-order corrections for a 𝐺 𝑊electron–hole interaction kernel, Phys. Rev. B 106 (2022) 045113, http://dx.doi.org/10.1103/PhysRevB.106.045113.[61] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev.B 47 (1993) 558–561, http://dx.doi.org/10.1103/PhysRevB.47.558.[62] J.P. Perdew, A. Zunger, Self-interaction correction to density-functional ap-proximations for many-electron systems, Phys. Rev. B 23 (1981) 5048–5079,http://dx.doi.org/10.1103/PhysRevB.23.5048.[63] S. Ishii, K. Ohno, Y. Kawazoe, S.G. Louie, Ab initio 𝐺 𝑊 quasiparticle energiesof small sodium clusters by an all-electron mixed-basis approach, Phys. Rev. B63 (2001) 155104, http://dx.doi.org/10.1103/PhysRevB.63.155104.[64] S. Ishii, K. Ohno, V. Kumar, Y. Kawazoe, Breakdown of time-reversal symmetryof photoemission and its inverse in small silicon clusters, Phys. Rev. B 68 (2003)195412, http://dx.doi.org/10.1103/PhysRevB.68.195412.[65] S. Ishii, S. Iwata, K. Ohno, All-electron 𝐺 𝑊 calculations of silicon, diamond, andsilicon carbide, Mater. Trans. 51 (2010) 2150–2156, http://dx.doi.org/10.2320/matertrans.M2010303.[66] M. Zhang, S. Ono, K. Ohno, All-electron 𝐺 𝑊 calculation of rutile TiO2 with andwithout Nb impurities, Phys. Rev. B 92 (2015) 035205, http://dx.doi.org/10.1103/PhysRevB.92.035205.[67] M. Zhang, S. Ono, N. Nagatsuka, K. Ohno, All-electron mixed basis 𝐺 𝑊calculations of TiO2 and ZnO crystals, Phys. Rev. B 93 (2016) 155116, http://dx.doi.org/10.1103/PhysRevB.93.155116.12 [68] T. Ishikawa, R. Sahara, K. Ohno, K. Ueda, T. Narushima, Electronic structureanalysis of light-element-doped anatase TiO2 using all-electron 𝐺 𝑊 approach,Comput. Mater. Sci. 220 (2023) 112059, http://dx.doi.org/10.1016/j.commatsci.2023.112059.[69] S. Ono, Y. Noguchi, R. Sahara, Y. Kawazoe, K. Ohno, TOMBO: All-electron mixed-basis approach to condensed matter physics, Comput. Phys. Comm. 189 (2015)20–30, http://dx.doi.org/10.1016/j.cpc.2014.11.012.[70] W. Kang, M.S. Hybertsen, Quasiparticle and optical properties of rutile andanatase TiO2, Phys. Rev. B 82 (2010) 085203, http://dx.doi.org/10.1103/PhysRevB.82.085203.[71] There is a typo in the last denominator term both in Eq. (2.18) of Ref. [51] andin Eq. (20) of Ref. [57]: 𝜀𝑑 (𝐸QP𝑣 ) in their equation should read 𝜀𝑑′ (𝐸QP𝑣′ ). In thecorresponding expression using the GPP model, there is also a typo in the lastdenominator term in Eq. (23) of Ref. [57]: 𝐸QP𝑣 in their equation should read𝐸QP𝑣′ .[72] L. Chiodo, J.M. Garcia-Lastra, A. Iacomino, S. Ossicini, J. Zhao, H. Petek, A.Rubio, Self-energy and excitonic effects in the electronic and optical propertiesof TiO2 crystalline phases, Phys. Rev. B 82 (2010) 045207, http://dx.doi.org/10.1103/PhysRevB.82.045207.[73] Y. Tezuka, S. Shin, T. Ishii, T. Ejima, S. Suzuki, S. Sato, Photoemission andbremsstrahlung isochromat spectroscopy studies of TiO2 (Rutile) and SrTiO3, J.Phys. Soc. Japan 63 (1994) 347–357, http://dx.doi.org/10.1143/JPSJ.63.347.[74] K. Vos, H.J. Krusemeyer, Low temperature electroreflectance of TiO2, Solid StateCommun. 15 (1974) 949–952, http://dx.doi.org/10.1016/0038-1098(74)90701-7.[75] M. Landmann, E. Rauls, W.G. Schmidt, The electronic structure and opticalresponse of rutile, anatase and brookite TiO2, J. Phys.: Condens. Matter. 24(2012) 195503, http://dx.doi.org/10.1088/0953-8984/24/19/195503.[76] A. Thatribud, Electronic and optical properties of TiO2 by first-principle calcu-lation (DFT-GW and BSE), Mater. Res. Express 6 (2019) 095021, http://dx.doi.org/10.1088/2053-1591/ab2cad.[77] M. Cardona, G. Harbeke, Optical properties and band structure of wurtzite-typecrystals and rutile, Phys. Rev. 137 (1965) A1467–A1476, http://dx.doi.org/10.1103/PhysRev.137.A1467.[78] See, for example, Appendix B of Ref. [44]. Also, there are discussions after Eq.(2) of Ref. [79] and after Eq. (9) of [80].[79] S. Albrecht, G. Onida, L. Reining, Ab initio calculation of the quasiparticlespectrum and excitonic effects in Li2O, Phys. Rev. B 55 (1997) 10278–10281,http://dx.doi.org/10.1103/PhysRevB.55.10278.[80] T. Rangel, M.D. Ben, D. Varsano, G. Antonius, F. Bruneval, F.H. da Jornada,M.J. van Setten, O.K. Orhan, D. O’Regan, A. Canning, A. Ferretti, A. Marini,G.-M. Rignanese, J. Deslippe, S.G. Louie, J.B. Neaton, Reproducibility in 𝐺0𝑊0calculations for solids, Comput. Phys. Comm. 255 (2020) 107242, http://dx.doi.org/10.1016/j.cpc.2020.107242.[81] G.E. Engel, B. Farid, Generalized plasmon-pole model and plasmon band struc-tures of crystals, Phys. Rev. B 47 (1993) 15931–15934, http://dx.doi.org/10.1103/PhysRevB.47.15931.[82] W. von der Linden, P. Horsch, Precise quasiparticle energies and Hartree–Fockbands of semiconductors and insulators, Phys. Rev. B 37 (1988) 8351–8362,http://dx.doi.org/10.1103/PhysRevB.37.8351.[83] L.J. Sham, T.M. Rice, Many-particle derivation of the effective-mass equationfor the Wannier exciton, Phys. Rev. 144 (1966) 708–714, http://dx.doi.org/10.1103/PhysRev.144.708.[84] Z. Zarhri, A. Benyoussef, A. El Kenz, Theoretical study of TiO2 doped withsingle and double impurities, J. Supercond. Nov. Magn. 27 (2014) 1323–1328,http://dx.doi.org/10.1007/s10948-013-2439-2.[85] G.C. Vásquez, D. Maestre, A. Cremades, J. Ramírez-Castellanos, E. Magnano,S. Nappini, S.Zh. Karazhanov, Understanding the effects of Cr doping in rutileTiO2 by DFT calculations and X-ray spectroscopy, Sci. Rep. 8 (2018) 8740,http://dx.doi.org/10.1038/s41598-018-26728-3.[86] G.Y. Csanak, H.S. Taylor, R. Yaris, Green’s function technique in atomic andmolecular physics, Adv. At. Mol. Phys. 7 (1971) 287–361, http://dx.doi.org/10.1016/S0065-2199(08)60363-2.http://dx.doi.org/10.1103/PhysRevB.48.17791http://dx.doi.org/10.1103/PhysRevB.50.10780http://dx.doi.org/10.1103/PhysRev.140.A1133http://dx.doi.org/10.1103/PhysRev.140.A1133http://dx.doi.org/10.1103/PhysRev.140.A1133http://dx.doi.org/10.1103/PhysRev.136.B864http://dx.doi.org/10.1103/PhysRevB.29.5718http://dx.doi.org/10.1103/PhysRevB.29.5718http://dx.doi.org/10.1103/PhysRevB.29.5718http://dx.doi.org/10.1103/PhysRevLett.75.818http://dx.doi.org/10.1103/PhysRevLett.75.818http://dx.doi.org/10.1103/PhysRevLett.75.818http://dx.doi.org/10.1103/PhysRevLett.80.3320http://dx.doi.org/10.1103/PhysRevLett.80.3320http://dx.doi.org/10.1103/PhysRevLett.80.3320http://dx.doi.org/10.1103/PhysRevLett.80.4510http://dx.doi.org/10.1103/PhysRevLett.81.2312http://dx.doi.org/10.1103/PhysRevLett.81.2312http://dx.doi.org/10.1103/PhysRevLett.81.2312http://dx.doi.org/10.1103/PhysRevB.59.5441http://dx.doi.org/10.1103/PhysRevB.62.4927http://dx.doi.org/10.1103/PhysRevB.62.4927http://dx.doi.org/10.1103/PhysRevB.62.4927http://dx.doi.org/10.1016/j.stam.2004.02.018http://dx.doi.org/10.1103/PhysRevA.81.045201http://dx.doi.org/10.1103/PhysRevB.106.045113http://dx.doi.org/10.1103/PhysRevB.106.045113http://dx.doi.org/10.1103/PhysRevB.106.045113http://dx.doi.org/10.1103/PhysRevB.47.558http://dx.doi.org/10.1103/PhysRevB.23.5048http://dx.doi.org/10.1103/PhysRevB.63.155104http://dx.doi.org/10.1103/PhysRevB.68.195412http://dx.doi.org/10.2320/matertrans.M2010303http://dx.doi.org/10.2320/matertrans.M2010303http://dx.doi.org/10.2320/matertrans.M2010303http://dx.doi.org/10.1103/PhysRevB.92.035205http://dx.doi.org/10.1103/PhysRevB.92.035205http://dx.doi.org/10.1103/PhysRevB.92.035205http://dx.doi.org/10.1103/PhysRevB.93.155116http://dx.doi.org/10.1103/PhysRevB.93.155116http://dx.doi.org/10.1103/PhysRevB.93.155116http://dx.doi.org/10.1016/j.commatsci.2023.112059http://dx.doi.org/10.1016/j.commatsci.2023.112059http://dx.doi.org/10.1016/j.commatsci.2023.112059http://dx.doi.org/10.1016/j.cpc.2014.11.012http://dx.doi.org/10.1103/PhysRevB.82.085203http://dx.doi.org/10.1103/PhysRevB.82.085203http://dx.doi.org/10.1103/PhysRevB.82.085203http://dx.doi.org/10.1103/PhysRevB.82.045207http://dx.doi.org/10.1103/PhysRevB.82.045207http://dx.doi.org/10.1103/PhysRevB.82.045207http://dx.doi.org/10.1143/JPSJ.63.347http://dx.doi.org/10.1016/0038-1098(74)90701-7http://dx.doi.org/10.1016/0038-1098(74)90701-7http://dx.doi.org/10.1016/0038-1098(74)90701-7http://dx.doi.org/10.1088/0953-8984/24/19/195503http://dx.doi.org/10.1088/2053-1591/ab2cadhttp://dx.doi.org/10.1088/2053-1591/ab2cadhttp://dx.doi.org/10.1088/2053-1591/ab2cadhttp://dx.doi.org/10.1103/PhysRev.137.A1467http://dx.doi.org/10.1103/PhysRev.137.A1467http://dx.doi.org/10.1103/PhysRev.137.A1467http://dx.doi.org/10.1103/PhysRevB.55.10278http://dx.doi.org/10.1016/j.cpc.2020.107242http://dx.doi.org/10.1016/j.cpc.2020.107242http://dx.doi.org/10.1016/j.cpc.2020.107242http://dx.doi.org/10.1103/PhysRevB.47.15931http://dx.doi.org/10.1103/PhysRevB.47.15931http://dx.doi.org/10.1103/PhysRevB.47.15931http://dx.doi.org/10.1103/PhysRevB.37.8351http://dx.doi.org/10.1103/PhysRev.144.708http://dx.doi.org/10.1103/PhysRev.144.708http://dx.doi.org/10.1103/PhysRev.144.708http://dx.doi.org/10.1007/s10948-013-2439-2http://dx.doi.org/10.1038/s41598-018-26728-3http://dx.doi.org/10.1016/S0065-2199(08)60363-2http://dx.doi.org/10.1016/S0065-2199(08)60363-2http://dx.doi.org/10.1016/S0065-2199(08)60363-2 Optical properties of rutile TiO2 with Zr, Mo, Zn, Cd impurities Introduction Method Results Discussion Conclusion CRediT authorship contribution statement Declaration of competing interest Acknowledgments Derivation of Wvc;v'c' Appendix. Derivation of Wvc;v'c' Data availability Appendix . Data availability References