{'fraction_non_alphanumeric': 0.028179205274053827, 'fraction_numerical': 0.01262235443088859, 'mean_word_length': 4.572804659498208, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Color centers in NaCl by hybrid functionals

We present in this work the electronic structure and transition energies (both thermodynamic and optical) of Cl vacancies in NaCl by hybrid density functionals. The underestimated transition energies by the semi-local functional inherited from the band gap problem are recovered by the PBE0 hybrid functional through the non-local exact exchange, whose amount is adjusted to reproduce the experimental band gap. The hybrid functional also gives a better account of the lattice relaxation for the defect systems arising from the reduced self-interaction. On the other hand, the quantitative agreement with experimental vertical transition energy cannot be achieved with hybrid functionals due to the inaccurate descriptions of the ionization energies of the localized defect and the positions of the band edges.


I. INTRODUCTION
The anion vacancy is a prototypical point defect in alkali halide crystals, which can be introduced by heating them in an alkali metal vapor or by x-ray irradiation and by electron bombardment. Anion vacancies in alkali halide are color centers (F centers) because they are responsible for the coloration of the otherwise transparent crystals. 1 The neutral anion monovacancy is the simplest form of the F center, with a single bound electron in the vacancy center. The electronic transition between the defect-induced levels can be trigged by the absorption of visible light. For example, the presence of the absorption peak (F band) associated with the F center at 465 nm renders the NaCl crystal yellow. [2][3][4][5] An eminent application of anion vacancies in alkali halides due to the coloration is the color center laser first demonstrated in 1965 by Fritz and Menke. 6,7 Due to their exotic features and long history, F centers in alkali halides have been studied extensively using various techniques such as optical absorption [8][9][10] , Raman spectroscopy [11][12][13][14] , the Stark effect [15][16][17] , and luminescence [18][19][20] . It is now well established that the defect states induced by anion vacancies in alkali halides are localized deep levels because of the small dielectric constant, large band gap, and high effective mass of the host crystals. These localized defect states cannot be correctly described within the effective mass approximation. 21 Early calculations employed either a linear combination of atomic orbital (LCAO) or the vacancy-centered variational wavefunction method within the Hartree-Fock (HF) theory. The F center was treated as a hydrogenic system trapped in the Coulomb potential created by the ions as point charges. [22][23][24][25] Although these calculations agree well with experimental data, it was later pointed out by Murrel and Tennyson, on basis of their LCAO self-consistent field (SCF) calculations within embedded cluster method, that the agreement could be for-tuitous, provided the results are sensitive to the selection of model Coulomb potential surrounding the anion vacancy. 26 Thus it is desirable to use a "true" potential from ab initio calculations. Numerous LCAO SCF cluster calculations have been performed at the restricted open-shell (ROHF) 26,27 and unrestricted HF (UHF) level. [28][29][30] The improvements, however, are limited and not unanimous as a result of the lack of electronic correlations and of the small cluster size limited by the computational capability.
Recent advances in density functional theory (DFT) makes it the tool of choice for studying the electronic structure of point defects in insulators. One major deficiency of DFT is that the standard local and semi-local approximations for the exchange-correlation (XC) functional suffer from the spurious electron self-interaction and the lack of the derivative discontinuities of the XC potential with respect to the particle number 31 . The calculated band gap is therefore generally smaller than experimental value, and the underestimation is more pronounced for wide band gap insulators. For NaCl, DFT calculations predict band gaps from 4.5 eV to 5.0 eV [32][33][34] , much smaller than the experimental value (8.5 eV) 35 . The band-gap problem poses serious uncertainties to the positions of the defect electronic levels and makes the credibility of the calculation questionable. One of the simplest yet the crudest way to overcome the gap problem is to apply a scissors operation to the conduction band and shift it rigidly to the experimental value relative to the valence band maximum (VBM). However, this still leaves ambiguities to the position of defect levels. In particular, for deep levels which do not follow either the characteristics of the VBM or conduction band minimum (CBM), the scissors scheme is not well-justified and not satisfactory in first-principles calculations. Other physically intuitive approaches to the problem include self-interaction corrections (SICs) 36 , and the DFT+U method where the localized states with a strong Coulomb repulsion, such as incomplete d or f shells in transitionmetal oxides, are modeled by a Hubbard-like U term 37 to open the band gap.
An alternative, and a more general approach to the gap problem is to relieve the unphysical self-interaction and derivative discontinuity problems by incorporating the non-local exact exchange (or Fock exchange) into the DFT XC functional. One of the implementation of the exact exchange is the screened exchange (sX) method, which replaces the whole local-density approximation (LDA) exchange by a Thomas-Fermi screened Coulomb exchange potential 38 . In hybrid density functionals, a fraction of exact exchange is mixed with LDA or generalized gradient approximation (GGA) exchange. The hybrid functionals generally not only improve the band gap, but also yield better results in bulk properties such as lattice parameters, bulk moduli and heats of formation for semiconductors and insulators 39,40 . A recent study comments that hybrid functionals are superior to SICs in reducing self-interaction errors since the unitary invariance is preserved. 41 Successful applications of hybrid functionals in the defect properties of various oxides and semiconductors have been recently reported. [42][43][44][45][46][47][48] In this work, we revisit the anion (Cl) vacancy in NaCl in different charge states (-1, 0 and +1) by both semilocal and hybrid functionals, aiming for a better understanding of the F center in alkali halide and the performance of the hybrid functionals for localized defects. We note that besides the band-gap problem, the finite-size effect arises within the current modeling scheme for the vacancy. The common procedure for defect energetics calculation is to embed the defect into a supercell under periodic boundary condition (PBC). The advantage of using supercells instead of cluster methods is that the band structure of the host crystal is well-defined, as the cell is bulk-like. 49 However, tractable DFT calculations are usually constrained to about 1000 atoms, and the size of system is further limited for hybrid functionals in a plane-wave basis set. A single Cl vacancy in a 1000atom NaCl supercell corresponds to a vacancy concentration of the order of 10 20 cm −3 , which is much higher than those found in experiment (10 15 to 10 19 cm −3 ). The periodic images of the point defects in a high density thus give rise to unrealistic defect-defect interactions, making the defect energetics dependent on the size of the supercells. The problem is even more serious for charged defects, as the neutralizing background slows the convergence of defect energies with respect to the supercell size. 50 Other sources of error for small supercells involve the elastic energy due to artificial relaxations of ions, and defect level dispersions introduced by defect-defect interactions. Corrections for the finite-size effect have been found indispensable in defect calculations for realistic interpretations [51][52][53][54] , and they will be discussed and applied to the present study of Cl vacancies.

II. CHOICE OF HYBRID FUNCTIONALS
We first assess the performance of the hybrid functionals on the bulk properties of NaCl, and compare it to the GGA-PBE (Perdew-Burke-Ernzerhof) functional 55 . Both unscreened and screened hybrid functionals are employed. In the unscreened PBE0 functional, the exchange part of the XC energy E xc is constructed by mixing a fraction (α) of non-local exact exchange E x with PBE exchange E PBE x , while the correlation energy is simply (1) The amount of exact exchange α is a variable from 0 to 1, although conventionally α = 0.25 is used as suggested by perturbation theory 56 . In practice α is usually varied to meet the experimental gap value. In a plane-wave basis set, the evaluation of the exact (HF) exchange is a hog to the computational resources and tends to be rather slow because of its truly non-local nature. The calculation can be accelerated by truncating the slowly decaying long-range part of the exact exchange as in the HSE (Heyd-Scuseria-Ernzerhof) hybrid functional 57 : ( The screening parameter µ in Eq. (2) determines the separation of the short-range (sr) and long-range (lr) parts: In one limit when µ = 0, the long-range term is zero and HSE reduces to the unscreened PBE0 functional. For µ → ∞, HSE is identical to GGA-PBE since the whole exact exchange is screened. Here we use the optimized µ = 0.207Å −1 , following Ref. 58 along with α = 0.25 and refer to this functional as HSE06.
In Table I selected bulk properties of NaCl calculated using the GGA-PBE and hybrid functionals are summarized together with the experimental values. The calculations are carried out in the projector augmented wave (PAW) framework with the vasp code [61][62][63] . A semicore pseudopotential of Na is used, treating the 2p3s electrons as valence electrons. The kinetic cutoff energy for the plane-wave basis set is 500 eV. A Γ-centered 8 × 8 × 8 Monkhorst-Pack k-point mesh 64 is applied to the primitive cell containing one formula unit of NaCl. For HSE06 and PBE0 calculations, a down-sampled 4 × 4 × 4 mesh is used to evaluate the non-local exact exchange. The down-sampling for the non-local exchange reduces the computing time significantly. It is generally necessary for the PBE0 to have a finer k-point mesh than for the screened HSE functional to reach convergence. 39,40 For the present case, the HF exchange using the PBE0 changes by roughly 15 meV per atom from the downsampled 4 × 4 × 4 to the full 8 × 8 × 8, while the energy is already converged within 10 −2 meV with the HSE06 functional. Nevertheless, the choice of the down-sampled k-point for the PBE0 calculations is sufficient for the bulk properties. The lattice constant a is determined when the residual force is smaller than 5 meV/Å. The high frequency macroscopic dielectric constant ǫ ∞ can be calculated within a GW scheme using the randomphase approximation (RPA) with local field effect included. Around 90 empty bands are used for calculating the dielectric constant. The dielectric constant will also be referred to later for the finite-size corrections. Finally, the formation energy ∆H f is obtained In Eq. (4), E NaCl(s) is the total energy of bulk NaCl. E Na(s) is the energy of bulk Na in a body centered cubic (bcc), which was optimized and calculated using the same k-point mesh and cutoff energy as the bulk NaCl. E Cl2(g) refers to the energy of one gas phase Cl 2 molecule in a large tetragonal cell. One immediately observes that the hybrid functionals improve not only the direct band gap (Γ 15 → Γ 1 ) but also the lattice constant and heat of formation compared to the GGA-PBE calculation in Table I, in agreement with earlier calculations. 39,40 Yet, it is found the band gaps are still underestimated for the hybrid functionals with the original fraction (0.25) of exact exchange, and the PBE0 yields a much closer value to experiment than the HSE06. This implies that for wide gap insulators, as the electronic screening is quite weak, the unscreened exact exchange in PBE0 is preferred. For defect calculations, it is customary to tune the fraction of the exact exchange so that the experimental band gap can be reproduced. 47,65 By increasing the amount of the non-local exchange from 0.25 to 0.40 within the PBE0, the band gap of NaCl recovers nearly to the experimental value, and the lattice constant and heat of formation are also reproduced best among the chosen functionals. We note that hybrid functionals tend to underestimate the dielectric constant of NaCl, a trend also found for semiconductors and other insulators. 66 An accurate description of the electronic dielectric constant with the hybrid functionals will require an explicit account of excitonic effects. 67 To this end, we face several functionals for the subsequent calculations of the Cl vacancy in NaCl. The PBE0 (α = 0.40) (we will refer it to mPBE0 hereafter) is apparently favored since it reproduces the experimental gap. However, as the choice of the fraction of the non-local exact exchange is empirical to some extent, its impact on the position of the deep defect level for wide gap insulators is still unknown. Meanwhile, the screened hybrid functional is of great interest as it shows considerable success in the prediction of defect energetics. Therefore it is plausible to also include the HSE06 functional with the original α, as well as the GGA-PBE for the defect calculations.

III. ELECTRONIC STRUCTURE OF CHLORINE VACANCIES
In this section we briefly sketch out the single-particle Kohn-Sham (KS) eigenvalues of the Cl vacancy induced electronic levels. Supercells containing 64 atoms are employed for the calculations. The Brillouin zone is sampled with a 2 × 2 × 2 k-point mesh, and a plane-wave basis set cutoff energy of 450 eV is used. Further, the k-point mesh for the non-local exact exchange is down-sampled to the Γ-point for HSE06, while a full 2 × 2 × 2 k-point mesh is necessary for well converged energies in PBE0 calculations. The convergence criterion for full relaxations is 0.01 eV/Å.
The removal of one Cl atom in a perfect NaCl crystal leaves a neutral vacancy V 0 with one electron bound to the vacancy center. The localized nature of the unpaired electron can be clearly identified in the charge density isosurface shown in Fig. 1. The 1s characteristics of the wavefunction in the vacancy is contributed equally from the six neighboring Na atoms. The negligible displacement of the neighboring atoms around V 0 (see Table II) keeps the singly occupied a 1g level unshifted after relaxation. In the +1 charge state V + , the a 1g state is unoccupied and the polaronic hole is trapped in the vacancy. The nearest neighbor Na atoms tend to relax away from the vacancy because of the positive electrostatic potential TABLE II. Displacements (in the equilibrium bond length 0.5a0) of the nearest-neighbor Na atoms around the Cl vacancy calculated using the 64-atom cell. The positive value represents an outward relaxation against the vacancy, and vice versa. The values obtained with the 216-atom cell using the GGA-PBE are shown in parentheses. inside the vacancy. The outward relaxation delocalizes the a 1g (see Fig. 1), and shifts it to higher energy towards the CBM. In the -1 charge state V − , the a 1g state becomes doubly-occupied. Upon relaxation the six nearest Na atoms show inward displacement towards the vacancy (Fig. 1) as a polaronic distortion. As a consequence, the two electrons are more localized inside the vacancy site, shifting the a 1g state to lower energy. We note that the relaxations of the neighboring atoms around the anion vacancy follow the O h symmetry for all charge states. No symmetry lowering (or Jahn-Teller distortion) is found since the defect level is either singly occupied for V 0 , or doubly occupied (unoccupied) for V − (V + ). The reduced self-interaction in hybrid functionals also results in more pronounced atomic displacements for the charged defects as seen in Table II. Table III summarizes the KS energies of the Cl vacancies in various charge states for both rigid and relaxed defect structures. The choice of k-point gives rise to a dispersion of the electronic level within the finite-size supercell scheme. The dispersion introduces a strong dependence on the supercell size of the energy level at the Γ-point ǫ Γ . Thus we average the defect level energy over the Brillouin zone ǫ since the averaged level shows a much better convergence than ǫ Γ . 51 The dispersion also slightly pushes the host CBM to higher energies. The finite-size effect will be discussed in detail in Sec. IV B. As predicted by all functionals, the a 1g states in all charged states lie within the upper half of the host band gap. We note that the absolute KS energies of the a 1g states [ǫ KS (a 1g ) + ǫ VBM ] are roughly unaffected when going from semi-local to hybrid functionals. In hybrid functionals, the VBM is lowered by 0.9 eV (HSE06) and 2.0 eV (mPBE0) with respect to the GGA-PBE as a result of the reduced self-interaction for the Cl 3p states. On the other hand, for rigid structures, the widening of the band gap in hybrid functionals tends to place the a 1g of V − further away from the CBM compared to the semilocal functional, while the unoccupied a 1g state of V + is closer to the CBM when hybrid functionals are used.

IV. THERMODYNAMIC TRANSITION ENERGIES AND FINITE-SIZE CORRECTIONS
In general the single-particle energy level of the defect as calculated from the KS equation differs from the experimentally observed transition energies. 51 A rigorous approach to the transition energies, as discussed in this section, relies on the total energy difference of the defect energetics in various charged states.

A. Formalisms of formation energies and transition energies
A central quantity for the defect energetics is the formation energy E f for a defect D in charge state q where E D and E H are the total energy of the supercell with the defect D, and the host supercell without defects, respectively. n i is the number of atoms removed from the supercell (n i = 1 for the Cl monovacancy) or the number of impurities added (n i < 0). µ i refers to the chemical potential of the associated defect particle reservoir, and is subject to equilibrium conditions. For the present study, under extreme Cl-rich (or equivalently Na-poor) conditions This places an upper limit on µ Cl . The lower bound can be deduced from the following relation: Therefore under Cl-poor (or Na-rich) conditions, which facilitate the formation of Cl vacancies, and this sets the lower limit of µ Cl . The remaining term ǫ VBM + ǫ F in the formation energy [Eq. (5)] represents the chemical potential, or Fermi energy of the electrons in charged defects. The Fermi energy ǫ F is varied within the band gap referenced to the energy of the host VBM ǫ VBM (0 ≤ ǫ F ≤ E g ). Here ǫ VBM is evaluated as the energy difference between a perfect host supercell and the same host supercell with one electron removed from the VBM: where n is the number of electron in the host supercell. In principle, one needs a sufficiently large supercell with n → ∞ corresponding to the dilute limit. In practice, a fractional charge q can be used along with a small supercell to obtain the ǫ VBM In the present case, the ǫ VBM converges well within a 64-atom cell and a fraction charge of 0.001 e.
For charged defects, it is evident from Eq. (5) that the formation energy is dependent on the chemical potential of the exchanged electron. The thermodynamic transition energy ǫ(q/q ′ ) is defined as the Fermi energy at which the charge state q and q ′ of the defect system can be transformed spontaneously from one to the other. Therefore at the transition energy ǫ(q/q ′ ) these two charge states have the same formation energy. This gives the following form of the transition energy B. Finite-size corrections Before proceeding to the results, we shall discuss the correction methods for the finite-size effect, whose causes have been already identified in Introduction. For charged defects, the simplest correction is to align the electrostatic potential in the defect supercell to that of the host supercell. This is usually done by inspecting the potential difference ∆V between the core potentials of the atoms far from the defect center and that of the bulk cell, and the energy correction term is essentially ∆E = q∆V . This correction is rationalized by the fact that in periodic supercell calculations the zero of the electrostatic potential is chosen arbitrarily for each calculation, and the charged defect gives rise to a constant shift in the potential so that the bulk VBM cannot be applied directly to the defect supercell. However, due to the small size of the supercell used even the atoms farthest from the charged defect center are not bulk-like, making such correction scheme inaccurate. Recent study reveals that the potential alignment resembles the Makov-Payne scheme, 50 whereas the latter targets the correction of the unphysical defect-defect interactions. Indeed, Komsa and Rantala found that ∆V has the form (ǫL) −1 , where L is the lattice constant of a cubic supercell. 69 This is analogous to the Makov-Payne scheme to first order.
The popular Makov-Payne scheme for a charged defect in a cubic supercell in the dilute limit (L → ∞) is expanded as (12) where α Md is the Madelung constant dependent on the lattice type, and q and Q the monopole and quadrupole moment of the defect charge, respectively 50 . The first order term, also called the Madelung energy, is thus the correction to the monopole-monopole interaction arising from the periodic image. The L −1 behavior of the artificial electrostatic interaction vanishes slowly, and this is usually the leading source of error. The higher order corrections have much smaller effects on the formation energy for ionic crystals, and it is usually accurate enough to include the corrections up to the quadrupole term. We note that in the Makov-Payne scheme the defect states are assumed to be localized, which is the case for the Cl vacancies in NaCl. For delocalized levels higher order corrections might become necessary. Although the Makov-Payne expansion is sound and accurate, it has been found that direct corrections using the Madelung energy and multipole interactions are prone to overshoot the formation energy, in particular for small supercells 70 . A more reliable approach is to employ a scaling method by performing a series of calculations using supercells of different sizes with the same symmetry. 71 The corrected formation energy E f (L → ∞) then can be extrapolated to the dilute limit by fitting the calculated formation energies within finite-size cells to where a n and E f (L → ∞) are fitting parameters. It is clear that this scaling law method requires at least 4 supercells, and is rather computationally laborious. In a recent work Freysoldt et al. proposed a general correction scheme (we will refer to it as the FNV scheme hereafter) for finite-size effect based on a single calculation of defect supercell without empirical parameters 72 : where E latt q is the macroscopically screened lattice energy of the defect charge q d with compensating background, and ∆ q/b is an alignment term referenced to the bulk supercell to account for the microscopic screening. As the long-range E latt q scales as L −1 and the short-range alignment term as L −3 , the FNV scheme can be seen as an extension to the Makov-Payne expansion. It also allows for an explicit expression for the third-order L −3 energy term. Now we apply both Makov-Payne scaling and FNV schemes to the formation energies of the charged Cl vacancies (V + and V − ) in NaCl. We refrain from including the potential alignment in these two schemes in order to avoid double-counting of the long-range L −1 term. Indeed Castleton et al. noticed that finite-size scaling with where N γ and N β are normalization constants, and x is the fraction of the relative amount of the exponential decay. In practice, the resulting corrected energy is insensitive to the choice of the specific parameters in Eq. (15). 72 In Fig. 2, we demonstrate the effects of finite-size corrections to the formation energies of Cl vacancies in +1 (V + ) and -1 (V − ) charge states. No relaxation is taken into account at this stage so as to exclude the finitesize effect of elastic energies. For V + , one first notices that the extrapolated formation energy from the Makov-Payne scaling law falls in line with that of the FNV scheme. The L −1 clearly dominates for the Makov  a rapid convergence of the V + formation energy. We see that finite-size correction is indeed mandatory for an accurate description of formation energy of charged defects. For the smallest 64-atom cell, the uncorrected formation energy is underestimated by roughly 0.6 eV. Even for the 1000-atom cell, the formation energy without correction is still 0.2 eV too low. Complexity arises when we move to the -1 charged Cl vacancy. The uncorrected formation energies in Fig. 2 for V − exhibit a zigzag evolution with respect to the increasing supercell size, making the Makov-Payne fit unreliable. Meanwhile, the FNV correction apparently yields a too high energy for small supercells, and the value does not appear to converge until we use the 512-atom cell. The source of such error is identified as the spurious dispersion of the defect levels as a result of the overlap between the wavefunctions of the defect and its periodic images. In the dilute limit, this localized defect level should be strictly a flat band. However, as shown in Fig. 3, the Cl vacancy induced a 1g level within the band gap shows a prominent dispersion for small supercells. For hybrid functionals, the dispersion is less pronounced as the exact exchange favors a more localized electronic state in the vacancy. This artificial interaction tends to push the CBM of the host crystal to higher energies, or otherwise the CBM will become populated. The defect level dispersion has a short-range characteristics, as is evidenced by the exponential fit of the a 1g bandwidth with respect to the supercell lattice constant L for the 64-, 216-and 512-atom cells (in Fig. 3). This short-ranged effect is not included in either the Makov-Payne or FNV scheme, and thus one has to take it into account explicitly. We note that the dispersion correction is not necessary for V + since its a 1g level is unoccupied. Here the correction for the dispersion is considered by calculating the energy difference between the KS energy of the a 1g level at the Γ-point ǫ Γ and the a 1g KS energy averaged over the sampled k-points in the Brillouin zone ǫ. The FNV scheme based on the dispersion corrected formation energies is again able to yield converged results for small supercells, and the results are also comparable to the Makov-Payne method including the dispersion effect in the dilute limit.
For the neutral Cl vacancy V 0 , the situation becomes trouble-free since the electron is tightly bound to the vacancy center with a strongly localized electron density distribution as seen in Fig. 1. The formation energy barely varies for the supercells considered, and therefore there is no need for finite-size corrections for V 0 .
We have shown the finite-size corrections for the charged defect supercells in rigid geometries with atoms fixed at their bulk positions. However, the introduction of a Cl vacancy inevitably changes the electrostatic potential of the local environment, resulting in atomic relaxations around the vacancy. The supercell approach, in this aspect, will lead to another error because the supercell employed in practice is usually not large enough for all local relaxations around the defect. This error can be partially alleviated by restricting the relaxations to the first two atomic shells around the defect, although it might underestimate the relaxation energy. Here we assess the finite-size effect on the elastic energy of the Cl vacancy in various charge states based on full relaxations using the GGA-PBE functional. We do not discuss the relaxations of the outer shell atoms since their displacements are much smaller than the first shell Na atoms and they contribute little to the formation energy. For the neutral vacancy V 0 , negligible inward relaxations of the six nearest neighbor Na atoms are found with supercells containing up to 216 atoms (see Table II). For V − , the net negative potential induced by the excess electron added to the vacancy gives rise to an inward displacement of the neighboring cations. As discussed in Sec. III, this results in a more localized a 1g state, which consequently suppresses the dispersion of the a 1g with respect to the rigid structure (see Fig. 3). Due to the finite-size of the supercell, Table II shows that the displacement obtained from the 64-atom cell is 0.05Å smaller than that from the 216-atom cell. The inability to fully relax in the 64atom cell consequently gives a formation energy about 0.2 eV higher than that of the larger supercells. For the positively charged vacancy V + , the six nearest neighbor Na atoms experience an outward displacement due to the positive potential in the vacancy center. In contrast to V − , the finite-size effect on the elastic energy is not significant for V + , as seen in Table II, since the atomic displacements using the 64-and 216-atom cells are of similar magnitude.
With all these comprehensive finite-size effects in mind, we now summarize the correction scheme applied in the present work. We restrict the calculations of formation energies with hybrid functionals to the use of the 64-atom supercell. Thanks to the localized nature of the defect state and negligible relaxation for V 0 , no correction is necessary. For V + we apply the FNV to the 64-atom cell and refrain from any correction for the elastic energy. For V − the dispersion correction is applied to the 64atom cell, followed by the FNV scheme. Further, the relaxation energies are aligned with those obtained from the 216-atom cell, provided the latter already yields a converged elastic energy. In practice, due to the similar amount of atomic displacement (see Table II), we use the PBE result as a reference for the hybrid functional, and subsequently lower the formation energies by 0.2 eV for the relaxed 64-atom supercells in the -1 charge state.

C. Chlorine vacancy thermodynamic transition energies
The calculated formation energies with finite-size corrections for the Cl vacancy in NaCl are shown in Table IV, with the Fermi energy ǫ F fixed at the VBM. For the neutral vacancy V 0 , the hybrid functionals yield higher formation energies than the GGA-PBE, although the energy differences are small. The functional dependence of formation energies for the charged states is much more prominent. With the hybrid functionals, we obtain higher formation energies for V − and lower formation energies for V + . In particular, the negative formation energy of V + under Na-rich conditions suggests that the F + center could be predominating when the Fermi energy is close to the VBM.
To trace the source of the functional dependence of formation energy for charged vacancies, we may first rewrite the formation energy of V + in a rigid geometry with ǫ F fixed at the VBM as where the first term on the right-hand side is the electron ionization energy of V 0 (or equivalently the affinity energy of V + ), and the second term simply the formation energy of V 0 . By decomposing the formation energy into several contributions, it is clear that the discrepancies in E f (V + ) stem mostly from the different positions of the VBM by various functionals. We note that the ionization energy of V 0 shows very small changes (within 0.05 eV) from semi-local to hybrid functionals, consistent with the similar absolute energy of the singly occupied a 1g state (see Table III). Analogously, the formation energy of V − can be rewritten as where the first term on the right-hand side is now the (negative) electron affinity energy of V 0 (or the ionization energy of V − ). In contrast to the ionization energy, it is found that hybrid functionals tend to yield a smaller affinity energy of V 0 than semi-local functionals, and that the difference reaches up to 0.9 eV. Along with the ǫ VBM , they explain the variations in the formation energy observed in Table IV. The atomic relaxation energy due to the polaronic electron or hole can be further extracted from Table IV as the difference between the rigid and relaxed structures. It is not surprising that the relaxation energies given by various functionals are consistent provided the atomic displacements are similar with these functionals (see Table II). The relaxation energy for V + is about 0.8 eV, while it ranges from 0.2 to 0.4 eV for V − .
For a charged defect, the formation energy is a function of the Fermi energy as illustrated in Fig. 4. The intersections of different charge states are the thermodynamic transition levels defined in Eq. (11). We see in Fig. 4 that the transition levels ǫ(+/0) and ǫ(0/−) are both within the band gap. Therefore, the mPBE0 functional predicts that all charge states (-1, 0 and +1) of Cl vacancy could be thermodynamically stable when the Fermi energy is varied within the band gap. We note that although this is also qualitatively predicted by the uncorrected formation energies, the Fermi energy window for the neutral V 0 vacancy is much narrower.
The density functional dependence of the thermodynamic transition levels is illustrated in Fig. 5(a). Both neutral and +1 charge states are predicted to be stable since the ǫ(+/0) levels are within the band gap for all functionals, with transition energies increasing from 3.81 eV to 5.36 eV as the band gap widens from semi-local to hybrid functionals. On the other hand, the ǫ(0/−) is placed slightly above the CBM in the PBE, while its   position falls into the band gap in HSE06 and is further shifted downwards with respect to the CBM in the mPBE0 calculations. Hence in contrast to the GGA-PBE, both the HSE06 and mPBE0 imply that the V − is stable. While the thermodynamic transition energy generally increases with respect to the VBM as the band gap enlarges, we see from Fig. 5(b) that the absolute position [ǫ(q/q ′ ) + ǫ VBM ] of the ǫ(+/0) level remains roughly unaffected from semi-local to hybrid functionals. This coincides with the findings by Alkauskas et al. that the calculated energy levels of localized defect are generally not tied to the position of the CBM. 73 The ǫ(0/−) levels are nevertheless more dispersed.

V. OPTICAL PROPERTIES OF THE COLOR CENTER
The experimentally available optical properties of the F (V 0 ) and F ′ (V − ) center in NaCl serve as a benchmark for the assessment of the performance of the functionals. The optical processes are clearly marked in the configuration coordinate diagram in Fig. 6 according to the Franck-Condon principle. 74 In the Franck-Condon approximation the electronic transition is assumed to occur very fast compared with the motion of nuclei in the lattice. Therefore, the optical excitation spectrum observed in experiment does not involve the relaxation of the defect structure, in contrast to the thermodynamic transition.
The optical absorption and emission can be described by vibronic (simultaneous vibrational and electronic) transitions, in which the lattice vibration mode is treated by a quantum harmonic oscillator. We first consider the excitation of an F center, which is a well-defined feature. [2][3][4]25 By absorption of a photon, the unpaired electron is transferred to an excited electronic state (V + state) and an excited vibrational state. The excitation of an electron into the CBM is equivalent to bringing an electron to a reservoir with a chemical potential of ǫ VBM + E g . The optical absorption energy E a , as illustrated in Fig. 6, is thus given by where the first term on the right-hand side is the formation energy of the unrelaxed +1 charge state (at coordinate R 0 ) with the Fermi energy at the CBM, and the second term the formation energy of the relaxed neutral charge state (at R 0 ). The excited state (F + center) subsequently relaxes to its zero-point vibration states. The energy gain due to the relaxation is the Stokes shift E S between the vertical absorption and the zero-phonon line (ZPL). The ZPL is the transition energy from the lowest vibrational level (zero phonon mode) of the ground state to the lowest level of the excited state (at R 1 ), without energy transfer to lattice phonons. In terms of thermodynamic transition energy, it is easy to see from Fig. 4 that the ZPL can be expressed in terms of the difference between the band gap E g and ǫ(+/0). In the present case, due to the identical formation energy between the rigid and relaxed V 0 , the Stokes shift reduces to the relaxation energy of the V + from the rigid structure. In the vertical emission (luminescence), the excited electron from the CBM recombines into the defect level and the emission energy E e is given by where the defect structure in the neutral state is kept fixed as that in the relaxed +1 charge state. Finite-size correction on this relaxation energy is taken into account according to Sec. IV B. Once the electron is in the ground state, it relaxes to the bottom of the state with the relaxation energy E AS (i.e. the anti-Stokes shift between the ZPL and the vertical emission). The pronounced Stokes and anti-Stokes shifts (see Table V) are expected due to the large polaronic distortion. The calculated vertical absorption and emission energies of the F center using the GGA-PBE and hybrid functionals are reported in Table V. The zero-point energy is not included since it is usually comparable for both the ground state and excited state. It should be borne in mind that the total energy difference scheme (∆SCF) based on the one-particle picture is not capable of describing the electron-hole coupling in the optical absorption. The presence of the exciton will usually introduce a uniform redshift of the absorption peak irrespective of the density functionals.
We find that the absorption and emission energies are underestimated by the GGA-PBE, and this is most likely related to the small band gap. At the present stage it remains unclear whether mPBE0 or HSE06 is more appropriate for the absorption energy since the excitonic binding energy of the F center is yet unknown. Meanwhile, the available experimental ZPLs for several F -aggregated centers (1.96 eV for R 2 band 78 and 1.48 eV for N band 79 ) suggest the mPBE0 might give a too large ZPL. Nevertheless, the hybrid functionals yield more realistic optical transition energies than the GGA functional. Analogously, we extend the calculation to the optical process of the F ′ center, which is formed when an electron is trapped at an F center by light absorption at low temperatures. 4 Instead of the sharp and bell-shaped curve of the F band, the F ′ center of NaCl gives rise to a broad F ′ absorption band, peaking at the longer wavelengths side of the F band. [2][3][4] It is seen in Table V that all functionals now predict smaller absorption energies with respect to the experimental F ′ band peak. This tendency is not changed even if the excitonic effect is taken into account as the electron-hole interaction will further decrease the absorption energy. In accord with the F band absorption, the F ′ E a increases from the semi-local functional to the hybrid functionals as the calculated band gap widens. The E a values given by the GGA-PBE and HSE06 are well below the experimentally observed peak, whereas the mPBE0 yields a value that is in better agreement with experiment. In addition, we find that the various functionals predict either a negative or a very small emission energy E e from the excited F ′ state to the ground state. A negative emission energy in Table V suggests that the configuration coordinate of the intersection lies between the coordinates of the minimum of the ground state and excited state. It is conceivable that in such case the excited state can return back to the ground state through a non-radiative process by vibrational relaxations, which leads to the luminescence quenching. The non-radiative path is also valid for a vibronic system with a small emission energy where the intersection is in the vicinity of the minimum of the excited state. 80 Experimentally, a radiative transition of an F ′ excited state is indeed absent. 81

VI. DISCUSSION
We have shown that, while hybrid functionals have been reported to be adequate for defects in some semiconductors, 42-48 the description of the localized anion vacancy in a wide gap insulator is less satisfactory by hybrid functionals when compared to the experimental optical absorption spectra. In this section, we aim to identify the possible origins of the failure of hybrid functionals for the description of the color centers in NaCl.
We start with the discussion of the optical absorption since it is well-defined experimentally. The absorption energy of an F center in Eq. (18) can be rewritten as where ǫ(+/0) R0 refers to the vertical transition energy occurring at the geometry for the neutral vacancy, and IP(V 0 ) is the ionization energy of V 0 . Therefore the absorption energy is solely dependent on the ionization of the neutral vacancy and the position of the band edge in a perfect supercell, and no structural relaxation is involved. As the ionization energy is not sensitive to the choice of the functional as discussed in Sec. IV C, it becomes obvious that the discrepancies in the absorption energy reported in Table V mostly stem from the variations in the CBM energy. For instance, the GGA-PBE ǫ CBM is 1.45 eV lower than the mPBE0 value as a result of the well-known band gap problem associated with the local and semi-local DFT functionals. While the energy gap can be reproduced by mixing 40% exact exchange in mPBE0, it is yet not clear whether the positions of the band edges are accurate. In principle, the band gap problem can be overcome by quasiparticle (QP) self-energy calculations based on many-body perturbation theory. Here we follow the widely adopted GW approximation for the electronic selfenergy 82 and calculate the QP corrections to the Kohn-Sham eigenvalues. The GW approximation can be understood as the Hartree-Fock theory with a dynamically screened Coulomb interaction. The QP energies are calculated in a two-atom NaCl unit cell with a Γ centered 4 × 4 × 4 k-point mesh and an energy cutoff of 200 eV for the response function, and a total of 256 bands. The dynamic dielectric matrix is constructed with a frequency grid of 200 points. 83 We note that the QP gap of NaCl is closely related to the starting wavefunction and selfconsistency. It is found that single shot G 0 W 0 correction is too small when it is applied to the GGA-PBE eigenstates, whereas G 0 W 0 on top of mPBE0 overestimates the QP gap. A fully self-consistent GW calculation also yields a too large QP gap irrespective of the initial eigenstates, as a result of the neglect of the attractive electron-hole interaction. 67 By updating the eigenvalues (four times) in the Green's function G and keeping the screened Coulomb interaction W at the RPA level within the initial PBE eigenvalues, the GW 0 @PBE scheme produces a QP gap of 8.43 eV, in agreement with experiment. The VBM is now lowered by 2.77 eV with respect to the PBE eigenvalue, and the CBM is lifted up by 0.66 eV, reaching to 4.90 eV by QP corrections. Compared to the GW 0 result, the mPBE0 CBM is placed 0.8 eV too high in energy, while the CBM energy calculated with the HSE06 functional coincides with that of the GW 0 (see Table III). A good agreement with the experiment F band absorption energy can be already obtained if the ǫ CBM in Eq. (20) is naively replaced by the GW 0 value while keeping the ionization energy untouched. Therefore, the ionization energies of the neutral Cl vacancy V 0 calculated by semi-local and hybrid functionals are well described from the total energy difference method. In contrast, the calculated ionization energy of the negative charge system V − is less satisfactory with GGA-PBE and is not much improved with the hybrid functionals based on the experimental F ′ band absorption peak and the GW 0 CBM energy. This is easily understood since the electronic correlations for the removal of a second electron from the a 1g level is beyond the scope of DFT. 91 These many-body effects are accessible from the manybody perturbation theory, e.g. in the GW approximation via the self-energy.
In the GW approximation, the QP energies of the highest occupied and lowest unoccupied level correspond to the electron removal and addition energies, respectively. It is then straightforward from Eq. (20) that the absorption energy can be calculated as the QP energy difference between the CBM and the lowest unoccupied state of the V + , provided the ionization potential of the neutral defect system is equivalent to the electron affinity of the positive charged system. For example, using a 64atom supercell with a cutoff energy of 100 eV and a Γ point for the response function and 1024 total bands, the GW 0 @PBE method yields an F band absorption energy of 2.47 eV. We note that the two-particle excitonic effect in the optical absorption is not taken into account in the GW approximation either.
For shallow defects, it has been found that the hybrid functional shows great improvement over local or semilocal functionals in the defect transition energies 48 . This is mostly likely benefited from the fact that the position of the shallow defect follows the band edge (either CBM or VBM), which can be reproduced by hybrid functionals with tunable α. For deep levels as demonstrated in this study, we find that a reproduction of a realistic band gap by an ad hoc tuning of the amount of the exact exchange in hybrid functionals does not guarantee an accurate description of the optical defect levels, and the thermodynamic charge transition levels as well. We see that the GGA-PBE is prone to an underestimation of the vertical transition energy, which is obviously impaired by a small band gap and a low conduction band edge. A significant shift for the band edges can be observed with the hybrid functionals. This leads to an increased vertical transition energy which is usually in better agreement with exper-iment, although the overestimation of the CBM energy by an increased fraction of the exact exchange in mPBE0 might give too high values (e.g. for the F band).
In addition to the electronic contributions (e.g. ionization energy and the position of the band edge), the structural relaxation also plays an important role in predicting the transition energies. The Frank-Condon shift (i.e. the difference between the absorption and emission energy) sheds light on the effect of the exact exchange on the lattice relaxations around the vacancy. By adopting the GW 0 CBM energy into the F band emission energy in Table V, we find that the lattice relaxation in presence of the electron-phonon interaction is best accounted by the mPBE0 hybrid functional as a result of the more localized electron density in the vacancy. The localization is proportional to the amount of the non-local exact exchange which reduces the self-interaction arising from the DFT XC functional. The localized nature of the trapped electron is further enhanced by a GW calculation, exhibiting an even smaller k-dispersion of the a 1g level than the mPBE0 since the GW approximation is free of self-interaction. In this context, we expect the structural relaxation will also be more realistic within many-body perturbation theory.
The discussion of the vertical transition energies finally invokes us to return to the thermodynamic transition level presented in Sec. IV C. The thermodynamic transition level can be readily decomposed into the electronic and structural contributions as where ∆E + D and ∆E 0 D are the relaxation energies from the initial geometry for the positive and neutral charge state, respectively. By incorporating the experimental vertical absorption energy and mPBE0 relaxation energy into Eq. (21) and (22), we come to the thermodynamic transition energies ǫ(+/0)=6.6 eV and ǫ(0/−)=7.5 eV. Therefore, the general picture of the Cl vacancy energetics by the adjusted hybrid functional so as to reproduce the experimental band gap remains qualitatively sound albeit not numerically accurate.

VII. SUMMARY
In this work we revisit the color centers in NaCl using hybrid density functionals. The reduced self-interaction error alleviated by the exact exchange in hybrid functionals allows for a more accurate description of the atomic relaxations in the vicinity of the vacancy compared to the DFT results. Yet, hybrid functionals are unable to achieve quantitative agreement with experimental optical absorption peaks. We show that this is closely related to the overestimations of the band edges when the fundamental gap is reproduced by an empirical amount of exact exchange. Meanwhile, when the electronic correlation comes into play during the removal (addition) of a second electron from (to) the localized defect level, the ionization (affinity) energy predicted within the framework of DFT is far from satisfactory. More elaborate methods (such as GW approximation 84-87 and two-particle Bethe-Salpeter equation [88][89][90][91] ) are thus necessary for more accurate descriptions of the optical process and defect energetics of the localized defects in wide gap insulators.