ADAQ-SYM: Automated Symmetry Analysis of Defect Orbitals

Point defects in semiconductors can provide a platform for solid-state quantum technology, with applications such as qubits[1, 2], sensors [3, 4] and single photon emitters [5, 6]. One significant benefit of quantum applications made with solid-state point defects is room temperature operation [5, 2, 7]. Theoretical calculations have been proven useful for identification of potentially interesting defects in wide-band gap semiconductors and quantitative estimations of their properties [8]. Indeed, first-principles methods based on density functional theory (DFT) can simulate the electronic structures and predict multiple properties [9, 10, 11]. Each semiconductor material may host a multitude of intrinsic and extrinsic point defects. To probe the large combinatorically complex chemical space in an efficient manner, high-throughput workflows have been developed [12, 13, 14] to simulate thousands of defect combinations, calculate relevant properties and store the results into a searchable database [12, 15]. Automatic Defect Analysis and Qualification (ADAQ) [14] is one such high-throughput workflow which currently has been used to screen about 52000 defects in 4H-SiC [16] and 21000 defects in diamond [17].

Symmetry-based characterization is a well-established practice in analyzing point defect systems to understand the selection rules determining the polarization of absorbed and emitted light, to better understand features of the electronic structure such as degeneracies, and in general help understand and explain the physics of defects. With symmetry analysis on the theoretical side, polarization specific PL measurements can be more accurately matched with simulated defects [18] and orientation for single defects can be identified [19]. Before analyzing the orbitals, the point group symmetry of the crystal hosting the defect needs to be found. There are two broadly used codes for this, spglib [20] and AFLOW-SYM [21]. We use AFLOW-SYM because of its reported lowest mismatch when identifying space groups for known crystals from the Inorganic Crystal Structure Database (ICSD) [21]. In addition, there are several codes that calculate irreducible representations of bands but mostly with focus on topological insulators [22, 23, 24, 25]. Within quantum chemistry, one method to quantitatively analyse the symmetry of molecular orbitals is with continuous symmetry measures (CSM) [26, 27, 28, 29], which provide a numerical measure of how close molecular orbitals are to certain irreducible representations. Defect orbitals in the band gap are localized much like molecular orbitals, yet methods similar to CSM have not been applied to point defects in solid host materials. Presently, a common method of symmetry analysis of defect orbitals is visually inspecting an isosurface of the wave function and how it behaves under the symmetry transformations, this may be prone to human error especially for high symmetry structures, and is not applicable in high-throughput workflows. Another method of analyzing the symmetry is to describe the defect orbitals as a linear combination of atomic orbitals and performing a group theory analysis by making comparisons to basis functions of the character table [30, 31]. This method may be possible to automate, however, since the structure has already been relaxed with plane waves we focus on analyzing these directly without projecting to atomic orbitals. By omitting the projection we can keep all information present in the plane waves.

This paper presents a quantitative symmetry analysis method for defect orbitals in solid host materials simulated with a plane wave basis set and the selection rules of optical transitions between defect orbitals.

We introduce ADAQ-SYM, a Python implementation of this method. The tool is fast and automated, requiring little user input, making it applicable as an analysis tool for simulations of defects, also with applicability for high-throughput screening of defects. Section 2 presents an introduction to the group theory, specifically applied to defects. Section 3 describes the ADAQ-SYM algorithm that performs the symmetry analysis, and Appendix A deals with how the software is constructed and what approximations are used. Computational details of the simulations in this paper are described in Section 4. Section 5 presents the results from symmetry analyses of several known defects; nitrogen vacancy (NV) center and silicon vacancy (SiV) center in diamond and the silicon vacancy ($\mathrm{{V_{Si}}}$) and several divacancy ($\mathrm{V_{Si}V_{C}}$) configurations in 4H-SiC. Section 6 discusses these results and Appendix B presents troubleshooting advice when using the tool.

In this paper we consider point groups. For convenience we summarize basic concepts following Ref. [32]. We refer to symmetry transformations as unitary transformations in three dimensional space which have at least one fixed point, meaning no stretching or translation. In the Schönflies notation, these transformations are:

• 1.

  Identity, E.

• 2.

  Rotation of $2\pi/n$ or $2\pi m/n$, where $n$ and $m$ are integers, $C_{n}$ or $C_{n}^{(m)}$.

• 3.

  Reflection in a plane, $\sigma_{x}$. $x=h$, $v$ or $d$, denoting reflection in a horizontal, vertical or diagonal plane.

• 4.

  Inversion, i.

• 5.

  Improper rotation of $2\pi/n$ or $2\pi m/n$, where $n$ and $m$ are integers, which is a rotation $2\pi/n$ or $2\pi m/n$ followed by a reflection in a horizontal plane, $S_{n}$ or $S_{n}^{(m)}$.

A set of these symmetry transformations, if they have a common fixed point and all leave the system or crystal structure invariant, constitutes the point group of that system or crystal structure. The axis around which the rotation with the largest $n$ occurs is called the principal axis of that point group. The point groups relevant to solid materials are the 32 crystallographic point groups, of which the following four are used in this paper, $\mathrm{C_{1h}}$, $\mathrm{C_{2h}}$, $\mathrm{C_{3v}}$ and $\mathrm{D_{3d}}$.

In brief, character describes how a physical object transforms under a symmetry transformation, (1 = symmetric, -1 = anti-symmetric, 0 = orthogonal), and representation $\Gamma$ describes how an object transforms under the set of symmetry transformations in a point group. Each point group has a character table which has classes of symmetry transformations on the columns and irreducible representation (IR) on the rows, with the entries in the table being characters. IRs can be seen as basis vectors for representations. Each point group has an identity representation, which is an IR that is symmetric with respect to all transformations of that point group. Character tables can have additional columns with rotations and polynomial functions, showing which IR they transform as. Appendix C contains the character tables of the point groups used in this paper, these character tables also show how the linear (x, y and z) and quadratic ($\mathrm{x^{2}}$, xy, …) polynomials transform.

For defects in solids, the point group is determined by the crystal structure, and the symmetry of the orbitals can be described by characters and IRs. Figure 1 shows a divacancy defect in silicon carbide with the point group $\mathrm{C_{3v}}$ as an example. Comparing with the character table for $\mathrm{C_{3v}}$, Table 4 in Appendix C, one sees that the orbital has the IR $a_{1}$.

When defects are simulated with DFT, one obtains (one-electron) orbital wave functions $\phi_{i}$ and corresponding eigenvalues $\epsilon_{i}$. Optical transitions, where an electron moves from an initial state with orbital i to a final state with orbital f, has an associated transition dipole moment (TDM) $\vec{{\mu}}$ , which is expressed as:

where e is the electron charge and $\vec{r}$ is the position operator. Selection rules can be formulated with group theory [32]. For TDM the following applies: for optical transitions to be allowed the representation of the TDM $\Gamma_{\mu}$ must contain the identity representation, where

with $\otimes$ being the direct product, and $\Gamma_{r}$ is the IR of the polarization direction of the light, corresponding to the linear functions in the character tables.

Figure 2 shows the symmetry analysis process of ADAQ-SYM. Here, we describe the steps in detail.

First, we perform a DFT simulation on a defect in a semiconductor host material. This produces a relaxed crystal structure and a set of orbital wave functions and their corresponding eigenvalues. These are the main inputs for ADAQ-SYM. The electron orbitals associated with defects are localized around the defect, and the inverse participation ratio (IPR) $\chi$ is a good measure of how localized an orbital is [33, 34, 35]. The discreet evaluation of IPR is

and can be used to identify defect orbitals in the band gap, since they have much higher IPR than the bulk orbitals. There are also defect orbitals in the bands which are hybridized with the delocalized orbitals, their IPR are lower than the ones in the band gap, but still higher than the other orbitals in the bands. We employ IPR as a tool for selecting the orbitals to be analyzed, and this method can also identify defect orbitals in the bands by spotting outliers. After the careful selection of ab initio data and inputs, ADAQ-SYM is able to perform the symmetry analysis.

Second, the ”center of mass” $\vec{c}$ of the each orbital is calculated according to

These centers are used as the fixed points for the symmetry transformations. Orbitals are considered degenerate if the difference in their eigenvalues are less than a threshold. When calculating $c$ for degenerate orbitals, they are considered together and the average center is used.

This method does not consider periodic boundary conditions and necessitates the defect be in the middle of the unit cell. To mitigate skew of the center of mass, the wave function is sampled in real space, and points with moduli under a certain percentage $p$ of the maximum are set to zero according to

Third, the point group and symmetry transformations of the crystal structure is found via existing codes. Each symmetry transformation has an operator $\hat{U}$. To get characters, the overlap of an orbital wave function and its symmetry transformed counterpart, the symmetry operator expectation value (SOEV), is calculated

for each orbital and symmetry transformation. The wave function is expanded in a plane wave basis set, with G-vectors within the energy cutoff radius. Therefore, Eq. 6 can be rewritten to be evaluated by summing over these G-vectors only once, the plane wave expansion is also truncated by reducing the cutoff radius when reading the wave function and renormalizing [36].

Fourth, the character of a conjugacy class is taken to be the mean of the overlaps of operators within that class, and the overlaps of degenerate orbitals are added. To find the representation of a set of characters, the row of characters is projected on each IR, resulting in how many of each IR the representation contains. Consider an IR $\Gamma$, where $\vec{W}_{\Gamma}$ is a vector with the characters of $\Gamma$ times their order. For example, for $\mathrm{C_{3v}}$ the vector for $a_{2}$ is $\vec{W}_{a_{2}}=(1\cdot 1,1\cdot 2,-1\cdot 3)=(1,2,-3)$. Let $\vec{V}$ be the vector with a row of characters and $h$ be the order of the point group, then $N_{\Gamma}$ is the number of times the IR $\Gamma$ occurs which is calculated as follows

For degenerate states the found representation should be an IR with dimension equal to the degeneracy, e.g. double degenerate orbital should have a two-dimensional $e$ state. If an IR is not found, the overlap calculation is rerun with the center of another orbital as the fixed point. The CSM $S$ for the IRs of molecular orbitals [28] is used for the defect orbitals and calculated with

which produces a number between 0 and 100. $S(\phi,\Gamma)=0$ means that the orbital is completely consistent with IR $\Gamma$, and $S(\phi,\Gamma)=100$ means that the orbital is completely inconsistent with the IR $\Gamma$.

Fifth, to calculate the IR of the TDM and find the allowed transitions, the characters of the TDM is calculated by taking the Hadamard (element-wise) product of the character vectors of each ’factor’

and Eq. 7 is used to calculate $\Gamma_{\mu}$. The representation of the resulting character vector is found in the same way as the IR of the orbital was found. As an example, consider the group $C_{3v}$ with the three IRs $a_{1}$, $a_{2}$ and $e$. If some TDM in this group has the character vector $\vec{V_{\mu}}=(4,1,0)$ calculating the representation would look like:

Since $\Gamma_{\mu}$ contains $a_{1}$, the transition is allowed. The software contains a function to convert a representation array of the above format to a string such as $"a_{1}+a_{2}+e"$.

Finally, the information produced by ADAQ-SYM is entered into a script to produce an energy level diagram which shows the position in the band gap, orbital occupation, IR and allowed transitions.

The DFT simulations are executed with VASP [37, 38], using the projector augmented-wave method [39, 40]. We apply the periodic boundary conditions, and the defects in the adjacent supercells cause a degree of self-interaction. To limit this, the supercell needs to be sufficiently large. In our case supercells containing more than 500 atoms are used. The defects are simulated with two different exchange-correlation functionals, the semi-local functional by Perdew, Burke and Ernzerhof (PBE) [41], and the hybrid functional by Heyd, Scuseria, and Ernzerhof (HSE06) [42, 43]. The PBE simulations are presented in the supplementary material, and the HSE simulations are presented in section 5. These simulations only include the gamma point, run with a plane-wave cutoff energy of 600 eV, with the energy convergence parameters $1\times 10^{-6}$ eV and $5\times 10^{-5}$ eV for the electronic and ionic relaxations, respectively. The simulations are done without symmetry constraints so symmetry breaking due to the Jahn-Teller effect can occur when relaxing the crystal structure. Excited states are simulated by constraining the electron occupation [44]. For cases where an orbital in the valance band is excited to a state in the band gap the setting LDIAG=.FALSE. is used.

To illustrate the capability of our method, we apply ADAQ-SYM to several defects in two different host materials, diamond and 4H-SiC, and analyze the symmetry properties for the defects orbitals. The symmetry analysis provides a coherent picture of the known defects, finds the allowed optical transitions between defect orbitals, and, specifically, explains the different ZPL polarization of the hk and kh divacancies in 4H-SiC. The following results are from HSE simulations, see the supplementary materials for results in PBE, as well as experimental comparisons of band gaps and ZPL energies.

From the symmetry analysis by ADAQ-SYM, one can attribute the differing polarization behavior of the $hk$ and $kh$ configurations to the symmetry of the lowest excited state, anti-symmetric and asymmetric (possibly symmetric) respectively. Due to the principal axis laying in-plane, different defect orientations will also emit out-of-plane with different polarization. Thus, it may be possible to experimentally determine the orientation of individual defects measuring the in-plane polarization angle of the PL detected along the c-axis, in an experiment similar to Alegre et al. [19]. In such an experiment, the $hk$ divacancy will exhibit a luminescence intensity maxima when the polarization is parallel to the principal axis, and a minima when the polarization in perpendicular. The opposite would be true for the $kh$ divacancy, and the two configurations could be distinguished by the approximately 30 meV difference in ZPL [18], or by the 30 degree polarization differences between the respective maxima, see supplementary material for polar plots of this polarization dependence.

Having a loose tolerance parameter for the AFLOW-SYM crystal symmetry finder can be useful in ambiguous cases since ADAQ-SYM will then run for a larger set of symmetry operators, which gives an overview and can provide insight to what extent the orbitals are asymmetric with regards to each operator. It is also recommended to do this when multiple gradual distortions of the same defect are examined. The initial PBE excited state calculation of the silicon vacancy in 4H-SiC seemed to show a case of the pseudo Jahn-Teller effect where the symmetry was reduced and the degenerate states split despite not being partially occupied in either spin channel. Upon running a simulation with more accurate tolerance parameters the point group remained $\mathrm{C_{3v}}$ and the splitting reduced to less than the threshold of 10 meV. This effect did not appear in the HSE simulation. For cases where the occupied state is very close to empty states convergence becomes more important and looser high-throughput simulations may exaggerate these effects. To resolve this one can either run more accurate calculations, or have a higher degeneracy tolerance parameter which will cause more states to be grouped together as degenerate.

We have presented a method of determining the symmetry of defect orbitals, and implemented this method in the software ADAQ-SYM. The implementation calculates the characters and irreducible representations of defect orbitals, the continuous symmetry measure is also calculated to get a numerical measure of how close the orbitals are described by the irreducible representations. Finally, ADAQ-SYM applies selection rules to the optical transitions between the orbitals. The tool is applicable to efficient analysis of defects. We have applied the tool to a variety of known defects with different point groups and host materials, and it reliably reproduces their symmetry properties. In addition to the reproduction of prior results, we find that the polarization of the allowed transition for $hk$ is parallel to an in-plane axis, and while $kh$ exhibits slightly broken symmetry, the polarization is perpendicular to an in-plane axis. A method to determine the orientation of individual $hk$ and $kh$ divacancies is also proposed. In summary, ADAQ-SYM is an automated defect symmetry analysis tool which is useful for both manual and high-throughput calculations.

For availability of ADAQ-SYM and instructions, see https://httk.org/adaq/.

William Stenlund: Conceptualization, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft. Joel Davidsson: Conceptualization, Methodology, Supervision, Writing - review & editing. Rickard Armiento: Conceptualization, Supervision, Writing - review & editing. Viktor Ivády: Conceptualization, Methodology, Supervision, Writing - review & editing. Igor A. Abrikosov: Conceptualization, Supervision, Writing - review & editing.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This work was partially supported by the Knut and Alice Wallenberg Foundation through the Wallenberg Centre for Quantum Technology (WACQT). We acknowledge support from the Knut and Alice Wallenberg Foundation (Grant No. 2018.0071). Support from the Swedish Government Strategic Research Area Swedish e-science Research Centre (SeRC) and the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Linköping University (Faculty Grant SFO-Mat-LiU No. 2009 00971) are gratefully acknowledged. JD and RA acknowledge support from the Swedish Research Council (VR) Grant No. 2022-00276 and 2020-05402, respectively. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at NSC, partially funded by the Swedish Research Council through grant agreements no. 2022-06725 and no. 2018-05973. This research was supported by the National Research, Development, and Innovation Office of Hungary within the Quantum Information National Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004) and within grant FK 145395.