Theory#
This chapter presents the theoretical formalism behind the computation of the luminescence properties of point defects in solids. It starts with a section dedicated to the theory of luminescence, showing how Fermi’s golden rule is approximated to obtain a tractable way of computing the luminescence lineshape, following either an effective phonon mode model or a multi-phonon mode model. Then, we present the computational methodology used to obtain the parameters entering the lineshape equation, and discuss how the use of forces instead of displacements allows one to increase the phonon supercell size. Finally, we present the IFCs embedding approach, which allows one to compute the phonon modes of a defect system in large supercells.
From Fermi’s golden rule to the luminescence lineshape#
Following Fermi’s golden rule, the absolute luminescence intensity \(I(\hbar\omega)\) (number of photons per unit time per unit energy) associated with one emitting center with two states, e and g, is expressed as a function of the photon energy \(\hbar\omega\) as
where \(n_D\) is the refractive index of the material and \(\boldsymbol{\mu_{eg}}=\langle\Psi_{e}|\boldsymbol{\mu}|\Psi_{g}\rangle\) is the total dipole matrix element. It is then assumed that the electronic part of the dipole matrix element between the excited and ground states depends weakly on the nuclear coordinates (Franck-Condon approximation). We also suppose that the nuclear motions can be expressed as a superposition of 3N harmonic normal modes of vibration \(\nu\), with normal coordinates \(Q_{\nu}\) and frequencies \(\omega_{\nu}\), N being the number of atoms in the cell. The total vibrational state \(\chi_{\boldsymbol{n}}\) is the product of 3N independent harmonic oscillator eigenfunctions \(\chi_{\boldsymbol{n_{\nu}}}\) with \(n_{\nu}\) the vibrational state of the \(\nu\)-th harmonic oscillator. We write the normalized luminescence intensity as
with C a normalization constant. In the case of the absolute intensity \(I(\hbar\omega)=C\omega^3A(\hbar\omega)\), \(C=\frac{n_D}{3 \pi \hbar \epsilon_0 c^3}|\boldsymbol{\mu}_{eg}^{el}|^2\), with \(\boldsymbol{\mu}_{eg}^{el}\) the purely electronic dipole moment (inversely proportional to the radiative lifetime of the emitting center). The emission lineshape function \(A(\hbar \omega)\), which dictates the shape of the spectrum, reads
In equation (26), \(\boldsymbol{n}\) denotes the set of 3N vibrational states {\(n_1\),\(n_2\),…,\(n_{3N}\)}, \(E_{ZPL}\) is the so-called zero-phonon line energy (purely electronic energy difference), \(E_{g/e,\boldsymbol{n/m}}=\sum{_\nu}n_{\nu}\hbar\omega_{g/e,\nu}\) is the vibrational energy of state \(\chi_{g/e,\boldsymbol{n/m}}\), and \(p_{m_{\nu}}(T)\) is a Bose-Einstein thermal occupation factor.
The vibrational modes of the excited and ground states can, in principle, be different, making equation (26) computationally very expensive because of the highly multidimensional integrals to be evaluated (Franck-Condon overlaps). Some approximations are therefore needed. We first assume that the vibrational modes in the excited electronic state are identical to those in the ground state. We also assume that we work at T=0 K, to simplify the equations (\(p_{m_{\nu}}(T=0)\) = 0 if \(m_{\nu}\), a vibrational state of the electronic excited state, is different from 0). Temperature-dependent luminescence spectra will be discussed later.
Effective phonon mode model#
The simplest model assumes that the 3N phonon modes of the system can be reduced to a single effective phonon mode. If a single effective mode couples to the electronic transition, the Franck-Condon overlaps reduce to \(|\langle\chi_{e,0}|\chi_{g,n}\rangle|^2\). In the harmonic approximation, the vibrational eigenfunctions \(\chi\) are expressed in terms of Hermite polynomials, and their overlaps can be computed analytically:
with \(S\) the Huang-Rhys factor, a dimensionless constant giving the mean number of effective phonons of frequency \(\Omega\) involved in the electronic transition,
and \(\Delta Q\) the offset between the harmonic potential energy surfaces of the ground and excited states. In practice, this offset is computed as the total mass-weighted atomic displacement induced by the electronic transition:
where \(i\) labels the Cartesian axes, \(\kappa\) the atoms, \(m_{\kappa}\) the atomic masses, and \(R_{e;\kappa i}\) and \(R_{g;\kappa i}\) are the atomic positions in the excited and ground states, respectively.
This allows us to simplify equation (26) to
A schematic representation of this model for a single mode \(\nu\) is presented in Fig. 1.
Fig. 1 Schematic representation of the origin of photoluminescence (PL) spectra. On the left, ground state (GS) and excited state (ES) energy curves are projected along phonon normal coordinate \(Q_{\nu}\) (here only for one mode) and are approximated by harmonic functions with the same frequency \(\omega_{\nu}\). Vibrational energy levels and corresponding eigenfunctions are shown as horizontal lines and colored areas. GS and ES are displaced by \(\Delta Q_{\nu}\), the mass-weighted displacement between the minimum of the GS and the ES curves. On the right, the PL spectrum is formed. The zero-phonon line (ZPL) comes from the transition between the first vibrational level of the ES to the first vibrational level of the GS. Other transitions give the phonon sideband (PSB). The intensity of each peak is computed with the overlap between corresponding eigenfunctions.#
Multi-phonon mode model and generating function#
We now consider all the 3N vibrational modes. The quantities of the previous section are now computed “per phonon mode”. The \(\Delta Q_{\nu}\) (mass-weighted atomic displacements projected along phonon mode \(\nu\))
are used to compute the partial Huang-Rhys factors \(S_\nu=\frac{\omega_\nu\Delta Q_\nu^2}{2\hbar}\), where \(e_{\nu,\kappa\alpha}\) are the phonon eigenvectors. The partial Huang-Rhys factors define the Huang-Rhys spectral function:
When dealing with 3N phonons, a direct evaluation of equation (26) is impractical, and the so-called generating function approach is used instead. The lineshape function \(A(\hbar\omega)\) is evaluated as the Fourier transform of the generating function \(G(t)\):
with
where \(S(t)=\sum_\nu S_\nu e^{i\omega_\nu t}\) is the Fourier transform of the Huang-Rhys spectral function, and \(\gamma\) is a constant homogeneous Lorentzian broadening associated with each vibronic transition.
When working in the time domain, the effect of temperature (transitions involving initial vibrational states \(n_{\nu} \ne 0\), weighted by the Bose-Einstein occupation) can be included straightforwardly by rewriting the generating function as
with \(C(t,T)=\sum_\nu\overline{n}_\nu(T)S_\nu e^{i\omega_\nu t}\) the Fourier transform of the temperature-weighted Huang-Rhys spectral function, and \(\overline{n}_\nu(T)\) the average occupation number of the \(\nu\)-th phonon mode:
One can connect the multi-phonon mode methodology to the simpler effective phonon mode model presented in the previous section. Thanks to the orthonormality of the phonon eigenvectors, the total normal coordinate change \(\Delta Q\) is linked to the partial \(\Delta Q_{\nu}\) through:
which allows one to define the weight with which each mode contributes to the total atomic relaxation:
It is then possible to define an effective frequency as
and the total Huang-Rhys factor as:
Semi-classical approach#
Within a semi-classical formulation [Henderson and Imbusch, 2006], one can derive a formula for the full width at half maximum (FWHM) of the emission lineshape:
where \(S_{\mathrm{em}}\) and \(S_{\mathrm{abs}}\) are the Huang-Rhys factors associated with the emission and absorption processes, respectively, and \(\Omega_{\mathrm{g}}\) is the effective frequency of the ground state.
Averaging over the Bose-Einstein occupations of the initial vibrational states allows one to compute the temperature-dependent FWHM:
Computational methodology#
From the above approximations, we see that the photoluminescence lineshape of a point defect in a solid is uniquely defined by the Huang-Rhys spectral function \(S(\hbar\omega)\) and the zero-phonon line. To compute it from first principles, we need to obtain:
The zero-phonon line energy \(E_{ZPL}\), which is the total energy difference between the electronic excited and ground states of the system.
The atomic relaxation associated with the change in electronic state, obtained as the difference between the relaxed atomic positions of the excited state and those of the ground state.
The vibrational modes of the system. Note that, in the case of the effective phonon model, these are not required, since the effective mode is obtained from the atomic relaxation.
One way to obtain the first two ingredients is to use DFT with the \(\Delta\)SCF constrained occupation method. The vibrational modes can be obtained with finite differences or DFPT, or with the embedding methodology presented later in this chapter.
\(\Delta\)SCF constrained occupation method#
This method uses DFT with non-Aufbau electronic occupations to mimic the electron-hole interaction. Transition energies are computed as differences between total DFT energies obtained with different occupations. Note that the excited-state occupations are specific to the system under study. In the case of Eu\(^{2+}\), one of the seven 4f electrons of the spin-up channel is promoted to the next spin-up 5d state, as shown in Fig. 2. This figure also illustrates the connection with the configuration coordinate model.
Fig. 2 Schematic representation of the constrained occupation \(\Delta\)SCF method, as used in this work for simulating the luminescent properties of Eu\(^{2+}\) phosphor.#
The procedure is as follows:
Start with a relaxed ground state (labeled \(A_g\)) with energy \(E_{g}\).
Promote an electron to obtain the excited state configuration (\(A_g^*\), energy \(E_{g}^{*}\)). The difference in energy with the \(A_g\) state provides an estimate of the semi-classical absorption energy:
(43)#\[ E_{\mathrm{abs}} = E_{g}^{*} - E_{g} \]After structural relaxation in the excited state, obtain the relaxed excited state (\(A_e^*\), energy \(E_e^*\)), which allows computation of the zero-phonon line (ZPL) energy:
(44)#\[ E_{\mathrm{ZPL}} = E_e^* - E_g \]The Franck-Condon relaxation energy of the excited state is:
(45)#\[ E_{\mathrm{FC,e}} = E_{g}^{*} - E_{e}^{*} \]De-promoting the electron with the relaxed atomic positions of the excited state yields the ground state with excited geometry (\(A_e\), energy \(E_e\)). The emission energy is:
(46)#\[ E_{\mathrm{em}} = E_e^* - E_e \]The Franck-Condon relaxation energy of the ground state is:
(47)#\[ E_{\mathrm{FC,g}} = E_{e} - E_{g} \]
Assuming harmonic behaviour, the effective frequencies of the ground and excited states can be extracted as
where \(\Delta Q\) is the mass-weighted atomic displacement between the ground and excited states, as defined in equation (29). The corresponding Huang-Rhys factors can be obtained from (28).
Forces vs displacements: increasing the supercell size#
For a given phonon mode \(\nu\), the corresponding partial Huang-Rhys factor \(S_{\nu}=\frac{\omega_{\nu}\Delta Q_{\nu}^2}{2\hbar}\) is computed from the mass-weighted displacement between the ground and excited states projected onto this phonon mode:
where \(\kappa\) labels the atoms (the Cartesian index \(\alpha\) is implicit). Within the harmonic approximation,
with the interatomic force constants (IFCs) \(C_{\kappa,\kappa'}=\frac{\Delta F_{\kappa'}}{\Delta R_{\kappa}}\), Eq. (49) can be rewritten as:
The use of Eq. (51) allows one to increase the supercell size and hence minimize finite-size effects.
Indeed, the forces decay to zero faster than the displacements as a function of the distance from the defect, as illustrated in Fig. 3, panels a-b.
Because of this rapid decay, we assume that the forces computed within the small red supercell (here 288 atoms, see Fig. 3, panel c) are already converged and free of finite-size effects, and that the forces outside this red supercell are essentially zero.
This means that we can use the forces computed in the red supercell and deduce the displacements at much larger distances in the blue supercell via the IFCs computed within this large blue supercell.
This subtle point is automatically taken into account when using Eq. (51) in the blue supercell. It also means that the IFCs/phonon modes should be computed in the same blue supercell, as discussed in the next section.
Fig. 3 (a) Norm of the displacement induced by an electronic transition as a function of the distance from the defect atom (Eu). (b) Norm of the forces in the ground state at the excited state equilibrium atomic positions. The decay of forces (b) is much faster than the decay of displacements (a). © Cartoon of a small red supercell containing a substitutional defect, in which the forces are computed with DFT. Outside this red supercell, the forces are set to zero because of their short-range decay, while the displacements are non-zero and are obtained with the IFCs computed in the large blue supercell.#
IFCs embedding#
In the computation of the partial Huang-Rhys factors, one would like to obtain the phonons in very large supercells at wave vector \(\mathbf{q}\rightarrow(0,0,0)\) (\(\Gamma\) point), so that the coupling with long-wavelength phonons is correctly captured. One would also like to include the coupling with the (localized) phonon modes introduced by the defect. A direct approach, either with DFPT or finite differences, is not computationally attractive. Indeed, DFPT is best suited to small pristine primitive cells with dense \(\mathbf{q}\)-meshes (which correspond to large pristine supercells after a folding procedure), while finite differences, which have no problem with the introduction of a defect, remain costly for large supercells. One way to include the effect of the defect (local modes) while converging the long-wavelength phonons is to employ the IFCs embedding approach [Alkauskas et al., 2014, Jin et al., 2021].
Fig. 4 Schematic view of the interatomic force constants (IFC) embedding approach to obtain defect phonons in large supercells. Each arrow represents an IFC between a pair of two atoms. See text for details.#
The procedure, illustrated in Fig. 4, is as follows:
First, the real-space interatomic force constants (IFCs) of the defect system, denoted as \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{defect}}\), are computed within a relatively small supercell (typically with a few hundred atoms). This is often done with finite differences on the supercell at the \(\Gamma\) point, although Density Functional Perturbation Theory (DFPT) could also be employed if the memory requirements are manageable and the level of theory required to describe the defect (e.g., DFT+U) is supported by the DFPT code.
Second, the IFCs of the pristine system, \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{pristine}}\), are obtained within a much larger supercell (typically with a few thousand atoms). A practical approach is to use DFPT on the unit cell to calculate the dynamical matrix on a \(\mathbf{q}\)-mesh, possibly followed by Fourier interpolation to obtain the matrix on a denser [\(N_x,N_y,N_z\)] \(\mathbf{q}\)-mesh. The information in the dynamical matrix computed on a unit cell for a dense [\(N_x,N_y,N_z\)] \(\mathbf{q}\)-mesh is equivalent to that in the dynamical matrix of an [\(N_x,N_y,N_z\)] supercell at \(\mathbf{q} = [0,0,0]\) (the \(\Gamma\) point). Mapping a \(\mathbf{q}\)-mesh to \(\mathbf{q} = [0,0,0]\) is referred to as a folding procedure, and yields the pristine IFCs \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{pristine}}\).
Third, an embedded IFC matrix \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{emb}}\) is constructed with the following rules: if both atoms \(\kappa\) and \(\kappa’\) are within a sphere of cut-off radius \(R_c\) centered on the defect, then \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{emb}}\) is set to \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{defect}}\). For all the other atomic pairs, the embedded IFC is set to the pristine value: \(C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{emb}} = C_{\kappa\alpha,\kappa’\alpha}^{\mathrm{pristine}}\). The cut-off radius \(R_c\) is typically determined by the size of the initial supercell used to compute the defect IFCs. Additionally, one could use a second cut-off radius \(R_b\): if two atoms are not within a sphere of radius \(R_b\), the IFC is set to zero. This second cut-off has been used in some studies [Jin et al., 2021, Razinkovas et al., 2021] to obtain large sparse matrices, for which dedicated techniques allow faster diagonalization (see, for instance, reference [Razinkovas et al., 2021]). In this work, we do not use this second cut-off.
Overall, the whole procedure can break the acoustic sum rule:
To correct for this, we follow the same approach as in [Jin et al., 2021, Razinkovas et al., 2021]:
As for the Born effective charges (BECs), if the BEC of the dopant atom was computed in the defect phonon calculation, it replaces the BEC of the substituted atom. If not, one can use the most common oxidation state of the dopant atom as its new BEC. The same approach applies to interstitials but, in this case, the dopant atom and its BEC are simply added. For vacancies, the atom and its BEC are removed.
Diagonalizing the embedded IFC matrix provides the phonon modes of a defective [\(N_x,N_y,N_z\)] supercell, in which the effect of the defect on the IFCs is included through the embedding procedure.
The localization of the phonon modes can be characterized by computing the inverse participation ratio (IPR), defined as [Alkauskas et al., 2014, Jin et al., 2021]
where \(\mathbf{e}_{\nu,\kappa}\) is the eigenvector of the phonon mode \(\nu\) associated with atom \(\kappa\). The IPR is a measure of the localization of the phonon mode \(\nu\), and is normalized to the number of atoms \(N\) in the supercell. Roughly speaking, the IPR can be interpreted as the number of atoms that participate in the phonon mode \(\nu\). For example, \(\mathrm{IPR}_{\nu}=1\) means that only one atom vibrates, and the phonon mode is therefore maximally localized, while \(\mathrm{IPR}_{\nu}=N\) means that all the N atoms in the supercell vibrate with the same amplitude.
The localization ratio \(\beta_{\nu}\) is defined as the ratio between the total number of atoms in the supercell, \(N\), and the \(\rm{IPR}\): \(\beta_{\nu} =N/\mathrm{IPR}_{\nu}\), where \(\beta_{\nu} \approx 1\) represents a bulk-like delocalized mode, while \(\beta_{\nu} \gg 1\) corresponds to a quasi-local or local mode.
References#
Brian Henderson and G Frank Imbusch. Optical spectroscopy of inorganic solids. Volume 44. Oxford University Press, 2006.
Audrius Alkauskas, Bob B Buckley, David D Awschalom, and Chris G Van de Walle. First-principles theory of the luminescence lineshape for the triplet transition in diamond NV centres. New J. Phys., 16(7):073026, 2014.
Yu Jin, Marco Govoni, Gary Wolfowicz, Sean E Sullivan, F Joseph Heremans, David D Awschalom, and Giulia Galli. Photoluminescence spectra of point defects in semiconductors: validation of first-principles calculations. Physical Review Materials, 5(8):084603, 2021.
Lukas Razinkovas, Marcus W Doherty, Neil B Manson, Chris G Van de Walle, and Audrius Alkauskas. Vibrational and vibronic structure of isolated point defects: the nitrogen-vacancy center in diamond. Physical Review B, 104(4):045303, 2021.
Julien Bouquiaux, Samuel Poncé, Yongchao Jia, Anna Miglio, Masayoshi Mikami, and Xavier Gonze. Importance of long-range channel sr displacements for the narrow emission in sr [li2al2o2n2]: eu2+ phosphor. Advanced Optical Materials, 9(20):2100649, 2021.
Julien Bouquiaux, Samuel Poncé, Yongchao Jia, Anna Miglio, Masayoshi Mikami, and Xavier Gonze. A first-principles explanation of the luminescent line shape of srlial3n4: eu2+ phosphor for light-emitting diode applications. Chemistry of Materials, 35(14):5353–5361, 2023.