Spin–orbit coupling and Jahn–Teller effect in Td symmetry: an ab initio study on the substitutional nickel defect in diamond

We analyse the spin–orbit and Jahn–Teller interactions in Td symmetry that are relevant for substitutional transition metal defects in semiconductors. We apply our theory to the substitutional nickel defect in diamond and compute the appropriate fine-level structure and magneto-optical parameters by means of hybrid density functional theory. Our calculations confirm the intepretations of previous experimental findings that the 2.56 and 2.51 eV optical centres are associated with this defect. Our analysis of the electronic structure unravels possible mechanisms behind the observed optical transitions and the optically detected magnetic resonance signal, too. This article is part of the Theo Murphy meeting issue ‘Diamond for quantum applications’.


I. INTRODUCTION
Nickel is a typical contaminant in high-temperature high-pressure synthesized (HPHT) diamonds that paved the way to observe various nickel related optical signals in diamond in the past decades [1,2].However, various assignments of optical signals to defect structures were only tentative assumptions.Recent advances in the development of ab initio methods made it feasible to revisit the optical centres with high prediction power and identify the microscopic structure of the optical centres.As a recent example, we have identified the nickel-vacancy complex by means of advanced density functional theory (DFT) methods [3] as the origin of the 1.4-eV photoluminescence (PL) centre [4][5][6] which was previously associated with the positively charged nickel interstitial [7][8][9][10].Nickel-vacancy has D 3d symmetry where group theory analysis combined with DFT calculations showed the electron-phonon coupling within Jahn-Teller (JT) effect formalism intertwins with spin-orbit coupling, and the full complexity of the system should be considered to obtain the key magneto-optical parameters for identification of nickel-related optical centres [3].In this paper, we apply this approach to the substitutional nickel defect (Ni s ) in diamond which has T d symmetry.Ni s introduces d orbitals to the electronic structure which presumably also a subject of spin-orbit coupling, and partial occupation of degenerate d states could lead to JT effect.As we shall show below, a very general group theory analysis of degenerate orbitals in T d symmetry is required for understanding this type of defect systems.
The Ni s configuration in diamond was observed in the electron paramagnetic resonance (EPR) spectra which is labelled as W8 [11][12][13].After the discovery of W8 EPR centre [11], the milestone experiments in the identification of Ni s were (i) the fingerprint of 61 Ni isotope in the EPR spectrum with the hyperfine interaction between the electron spin and I = 3/2 nuclear spin in isotopically 61 Ni enriched HPHT diamond samples [12] and (ii) observation of S = 3/2 electron spin with hyperfine signals of four identical 13 C I = 1/2 nuclear spins [13].The W8 EPR centre exhibits an isotropic g = 2.0310 ± 0.0005 tensor.These facts imply that the defect has T d symmetry.Furthmore, S = 3/2 spin state indicates that it is a substitional defect in the negative charge state, i.e.Ni − s .This model was later confirmed by additional measurements and ab initio calculations [14][15][16][17][18][19].The ground state of Ni − s is well understood: Ni s introduces a five times degenerate d orbital of which level splits due to the tetrahedral crystal field of diamond as a triple degenerate one (t 2 ) lying in the fundamental band gap whereas the double degenerate d orbital (e) is resonant with the valence band.The in-gap t 2 level is occupied by three electrons in the negative charge state which establishes the paramagnetic S = 3/2 state.This electronic configuration is the 4 A 2 orbital singlet multiplet state in T d symmetry.
The correlation between the Ni-related optical centres and the W8 EPR centre were studied either indirectly such as common appearance in the appropriate optical and EPR spectra with similar estimated concentrations or directly via optically detected magnetic resonance (ODMR).
In the former method, the correlation between the 1.883-eV and 2.51-eV Ni-related absorption centres and the W8 EPR centre was investigated.It was concluded that the 2.51-eV absorption centre is likely linked to the W8 EPR centre rather than the 1.883-eV absorption centre [20].Initial uniaxial stress measurements indicated that a triple degenerate T 2 state is involved in the 2.51-eV absorption spectrum [21] which was later revisited and concluded that 4 A 2 ↔ 4 T 2 optical transition is involved in the absorption process [22], being consistent with the ground state of Ni − s .Although, the 4 T 2 state is not optically allowed from 4 A 2 ground state but it was assumed that spin-orbit coupling between the 4 T 2 and 4 A 2 makes this transition partially allowed.The shape of the first peak exhibit a doublet feature separated by 1.5 meV [22].They interpreted this feature to the spinorbit splitting of the 4 T 2 state.Beside this feature a replica shows up at around 16.5 meV above the zerophonon-line (ZPL) energy which is broadened compared to the shape of the ZPL emission [22].The origin of the replica was not explained.
In the latter method, the 2.56-eV PL centre could be directly linked to the W8 EPR centre via ODMR measurements at cryogenic temperatures [23,24].Two individual features were observed in the ODMR spectrum excited by 325-nm laser in a 35-GHz microwave cavity [24]: (i) an electron spin resonance with the isotropic g = 2.0324 (5) associated with the W8 EPR centre with producing a PL spectrum in the ODMR contrast which agrees well with the 2.56-eV PL spectrum in terms of the zero-phonon-line (ZPL) energy and the features in its phonon sideband; (ii) the hyperfine splitting originating from the 14 N I = 1 nuclear spin in the P1 EPR centre, i.e. the neutral substitutional nitrogen donor, N 0 s with g ≈ 2.00 which produces a very broad gaussianshape PL spectrum in the ODMR contrast.The ODMR contrast of the S = 1/2 P1 centre was explained by a donor-acceptor pair (DAP) model, where the unknown acceptor A interacts with the N s donor as where hν photon is emitted for which the energy depends the distance between the N s and the acceptor A, so it results in a broad fluorescence spectrum for the ensemble of DAPs.If the unknown acceptor is boron (acceptor level is at 0.37 eV above the valence band maximum, VBM) then the fluorescence spectrum can be well explained by taking the huge reconstruction energy of the N s upon ionisation into account [24].It was speculated that electron spin resonance of the acceptor in the ODMR spectrum was not observed because if the acceptor state is an effective-mass-like, then the degeneracy of the top of the VBM causing its EPR signal to be strongly strain broadened and difficult to detect, as in all cubic semiconductors [24].Regarding the Ni-related ODMR feature, it is more difficult to come up with a unique model that explains the observation of the ground state EPR of Ni, in the luminescence.Nazaré et al. speculated [24] that one possibility is spin-dependent hole transfer from a distant acceptor to Ni − s , producing an excited luminescent state of Ni 0 s .We note that this may explain that the observed isotropic g-factor in ODMR spectrum (2.0324) slightly differs from the g-factor in the W8 EPR spectrum (2.0310) [24].It is worth to note too that the 2.56-eV PL centre can be activated though a broad excitation at 200 ± 20 nm range (∼ 6.2 ± 0.6 eV) efficiently as performed by photo-luminescence excitation (PLE) measurements [25].This suggests that carriers trapping is indeed involved behind the ODMR process of Ni s defect.We note that the reported negative ODMR contrast of the 2.56-eV PL in Ref. 23 was about 10 −5 in their experimental conditions (excitation wavelength of 365 nm applied in 36.2-GHzmicrowave cavity) which showed up as a positive contrast without reporting its magnitude in Ref. 24 in their experimental conditions (excitation wavelength of 325 nm applied in 35-GHz microwave cavity).This fact also implies that the ODMR process is not intrinsic to the Ni s defect alone.Photo-EPR measurements were carried out in nickel doped diamonds and found that the intensity of W8 EPR signal decreases by photo-excitation with ∼ 2.5 eV threshold energy which can be reinstated by photo-excitation with ∼ 3.0 eV [26].

II. METHODOLOGY
We characterise Ni s in diamond by plane wave supercell calculations within spin-polarised DFT as implemented in the vasp 5.4.1 code [27].We determine the electronic structure within the Born-Oppenheimer approximation where the ions are treated as classical particles.We relaxed the atomic positions until the forces acting on the ions fall below 10 −2 eV/ Å.We embed the defect in a 512-atom diamond supercell and we sample the Brillouin-zone at the Γ-point.We used the projector-augmentation-wave-method (PAW) [28,29] as implemented in vasp.We used the standard PAW projector for carbon and Ni pv PAW projector for nickel that includes the 3p atomic orbitals as valence.We applied a plane wave cut-off energy at 370 eV for the plane wave basis that was already proven to be convergent for SiV defect and many other defect system in diamond [3,[30][31][32][33][34][35].We calculate the excited states with the constrainedoccupation DFT method or ∆SCF method [36].We employed the Heyd-Scuseria-Ernzerhof (HSE06) hybrid functional [37,38] which reproduces the experimental band gap and the charge transition levels in diamond or other group-IV semiconductors within 0.1 eV accuracy [36,39].For the charged supercell we applied the Freysoldt-Neugebauer-van de Walle correction to the total energy [40].We use our home-built code to compute the electron-phonon spectrum of the Jahn-Teller modells that we successfully applied to other systems too [33,35].

III. RESULTS
We combine group theory description of the Ni s defect with DFT calculations.We start with the calculated electronic structure and charge transition levels of the defect.We identify the most common charge states under typical experimental conditions in this section.We continue with the detailed discussion of the ground state manifolds of the most relevant charged states of the defect and then we describe the optically accessible excited states.Finally, we discuss our findings in the light of the observed 2.51-eV absorption centre and the 2.56-eV PL centre including the possible mechanisms behind the ODMR signature.

A. Electronic structure of Nis
The Ni atom introduces 3d orbitals which splits to double degenerate e orbitals and triple degenerate t 2 orbitals in the tetrahedral crystal field of diamond.In the defect molecule diagram, these 3d orbitals and the 4s orbitals interact with the four dangling bonds of the neighbour carbon atoms that results in a 1 and t 2 molecular orbitals under T d symmetry.The 4s and a 1 orbitals can form bonding and anti-bonding molecular orbitals lying deep in the valence band and high in the conduction band, respectively.The e orbitals do not recombine with the dangling bonds and they remain atomic like.The level of e orbitals is resonant with the valence band in the neutral charge state of the defect as plotted in Fig. 1(a).Finally, the t 2 representations of the 3d orbitals and the dangling bonds again form bonding and antibonding molecular orbitals.The bonding combination falls deep in the valence band whereas the antibonding combination appears in the gap.In the context, we call this antibonding t 2 orbital simply the t 2 state in the gap which is responsible for the ground state manifolds of the Ni s defect.In the optical excitation, the atomic like e states may play a role.The calculated Kohn-Sham states are depicted in Fig. 1(b).
According to this defect molecule analysis, the defect has the a 2 1 t 6 2 e 4 t 2−q 2 electronic configuration, where q is the charge state of the defect.As the triple degenerate t 2 may accommodate six electrons, the possible charge states may go from the (2+) to (4−).We calculated the total energies of the corresponding charge states (E q tot ) and calculated the adiabatic charge transition levels as (q|q − 1) = E q tot − E q−1 tot + qE VBM , where E VBM is the calculated VBM energy and now it is aligned to zero in Fig. 1(a).
We find that the aniticipated (3−) and (4−) states are not stable at all, and the (2+) and (2−) charge states are marginally stable as a function of the Fermi-level in diamond.Therefore, we do not analyze the (2+) and (2−) charge states further.The calculated (0|−) level is at about E VBM + 2.6 eV, which means that it is about E CBM − 2.8 eV, where E CBM is the energy of the conduction band minimum (CBM).In other words, if the defect was in the neutral charge state then about 2.6 eV photo-excitation energy is required to ionise to (−) charge state, and if the defect was in the (−) charge state then about 2.8 eV photo-excitation energy is required to re-ionise it back to the neutral charge state.Ni + s can be readily ionised to the neutral charge state by the photo-excitation threshold energy at about 1.5 eV.It can be read out from Fig. 1(a) that Ni + s may be stable in p-type diamond and near-infrared photo-excitation is needed below the photo-ionisation threshold energy in order to observe its absorption or fluorescence signal.In practice, these conditions are not fulfilled, thus we will only briefly discuss the ground state of Ni + s .We rather focus our discussion to the negatively charged and neutral Ni s defects.
We briefly discuss now the 2.51-eV absorption and 2.56-eV PL centres in the light of the calculated photoionisation threshold energies.According to the calculated (0|−) level at around E VBM + 2.6 eV, both optical centres might correspond to the neutral excitation of both Ni 0 s and Ni − s defects.The 2.56-eV PL centre has the ZPL energy just below the calculated Ni 0 s →Ni − s ionisation threshold energy.Thus, the calculated ionisation threshold energies could not distinguish alone which optical centre belongs to which defect charge state.In order to identify these optical centres, we analyze the ground and excited states of Ni − s and Ni 0 s in detail in the next sections.
It is intriguing that the 2.56-eV emission can be observed above-band-gap excitations [25], and its ODMR signals were observed by 325-nm (3.81-eV) photo-excitation [24] and 365-nm (3.40-eV) photoexcitation [23] experiments.In the latter two experiments, the photo-excitation energies are both above the ionisation threshold energies for both the Ni 0 s →Ni − s and Ni − s →Ni 0 s processes, where the electron is promoted from the valence band to the empty defect level and from the occupied defect level to the conduction band, respectively.The calculated (0|−) level at around E VBM + 2.6 eV implies that the electron is promoted from a deeper energy region of the valence band in the first process rather than it is excited from the occupied defect level to little above the conduction band minimum in the second process.As the electron density of states rapidly increases going farther from the band edges this implies that Ni 0 s →Ni − s process is much more efficient than Ni − s →Ni 0 s process.As a conseuqence, the ultraviolet excitation will stabilise the Ni − s charge state, so the defect stays much longer time in the (−) charge state than in the (0) charge state.

B. Ground state of Ni − s
It is likely that ultraviolet photo-excitation stabilises the negative charge state of Ni s , and the totally filled t 2 electronic configuration of Ni − s can be readily analyzed as a basis for the more complex electronic configurations.Therefore, we start the analysis with the ground state electronic structure of Ni − s .According to first Hund's rule the three electrons in the triple degenerate t 2 orbitals will form a high-spin S = 3/2 electronic configuration.The wavefunction of such quartet state can be expressed as Fermi energy (eV) where we define the three particle (S) symmetrisation and (A) anti-symmetrisation operators for the sake of simplicity.In addition, we use shortened (x, y, z) notation to label (t 2x , t 2y , t 2z ) orbitals in the context from now on.Our results are in agreement with previous ab initio studies [17,19], the ground state is a | 4 A 2 ⟩ multiplet.The S = 3/2 spin of W8 EPR centre emerges from the three unpaired electrons fillling the t 2 orbitals in Ni − s .
C. Shelving doublet states of t electronic configuration can form doublet states too as the total number of possible combinations of muliplet states is 6  3 = 20 where 4 A 2 quartet state produces only four multiplet states from them.The other configurations form doublet manifolds that should be expressed as a combination of Slater-determinants. First, we determine the possible configurations with single Slater-determinants.There are 8 possible occupations where we pin electrons on individual x, y, z spin-orbitals as where S z = ± 3 2 combinations and the totally symmetric S z = ± 1 2 combinations are reserved for the | 4 A 2 ⟩ ground state as that is equivalent to Eq. ( 2).Therefore, the remaining four-dimensional subspace spans a | 2 E⟩ state as By constraint occupation of the t 2 Kohn-Sham orbitals, the total energy of the system can be calculated by spin-polarised HSE06 functional as given in Eq. ( 5).The exchange-correlation functional of HSE06 may not be able to capture the high correlation which can be described as combinations of Slater-determinants and the resulting solutions will be not the true eigenstate of the doublet spin state.
With this caveat, we carried out ab initio calculations for the single Slater-determinant configuration, e.g.A x ↑ y ↑ z ↓ of the Kohn-Sham t 2 states which yields 0.35 eV with respect to the calculated total energy of 4 A 2 state.This was achieved by restricting the symmetry to T d symmetry.By enabling the symmetry reduction for the A x ↑ y ↑ z ↓ electronic configuration, the energy gain has become negligible at 1.4 meV.We conclude that there will be no Jahn-Teller effect for the doublet configurations in the first order.The A x ↑ y ↑ z ↓ determinant is an 33% admixture of | 4 A 2 ⟩ and 66% of | 2 E⟩ multiplets, thus we expect the level of | 2 E⟩ lying above that of | 4 A 2 ⟩ by roughly 0.35 eV.
Other possible configurations are which reduces into a | 2 T 1 ⟩ spin doublet as and an another | 2 T 2 ⟩ spin doublet as We expect that the levels of these states lie above those of the maximally open shell configurations' doublets because the Coloumb repulsion is not compensated by exchange interaction in the former.In Kohn-Sham DFT method, we could calculate one of the single Slater-determinants electronic configuration as expressed in Eq. (8).Their total energy exceeds that of | 2 E⟩ by about 0.5 eV within T d symmetry.
Finally, we conclude based on the Coulombic repulsion principle that | 2 T 2 ⟩ has the highest total energy among the considered doublet manifolds as it exhibits the most symmetric polynomial, that is, The ground state manifolds of Ni 0 s exhibit t (2) 2 electronic configurations in T d symmetry which may be a subject of Jahn-Teller effect and spin-orbit coupling.The total number of manifolds is 6 2 = 15.According to direct product tables [41] of T d point group, there are four multiplets that can be read as where we kept the fermionic anti-symmetric property of the wavefunctions, thus do not occur.Next, we show the symmetrically adapted wavefunctions.According to our ab initio results the ground state of Ni 0 s is a | 3 T 1 ⟩ that reads as It is followed by the | 1 T 1 ⟩ singlet state which can be expressed as where we interchanged the two-particle A |ab⟩ = The triplet configuration can be determined by Kohn-Sham DFT as a |x ↑ y ↑ ⟩ configuration.We note that | 3 T 1 ⟩ triplet is Jahn-Teller unstable.We obtained 159 meV relaxation energy when we removed the T d symmetry constraint during the geometry optimization.We obtained 47 meV higher total energy for the |x ↑ y ↓ ⟩ with T d symmetry restriction that gained 125 meV relaxation energy upon removal the symmetry constraint.We note that the |x ↑ y ↓ ⟩ configuration is the 50%-50% admixture of the | 3 T 1 ⟩ and | 1 T 2 ⟩ multiplets, thus our results can be interpreted as that the singlets are above the triplet by about ∼47 meV.
Eqs. ( 12) and ( 13) span a 9+3 = 12 dimensional openshell subspace where the electrons are placed to different orbitals within t 2 orbitals.Next, we discuss the case of the closed shell configurations (xx, yy, zz) which leads to the | 1 A 1 ⟩ and | 1 E⟩ multiplets which may read as and respectively.We note that | 1 E⟩ and | 1 A 1 ⟩ are multiconfigurational states, thus Kohn-Sham DFT may not provide good estimates for their total energies.Nevertheless, we can determine the total energy within spinaveraged Kohn-Sham DFT calculation for the individual (xx, yy, zz) configurations that forces the same orbitals in each spin channel.According to our results, the total energy of the closed shell configurations lies above that of the ground state by 0.70 eV within T d symmetry.
After removing the symmetry constraint in the geometry optimization procedure, the calculated energy difference is reduced to 0.18 eV that shows a giant 0.52 eV Jahn-Teller energy.This effect might suppress the electronic correlation energy between | 1 E⟩ and | 1 A 1 ⟩, similarly to the product Jahn-Teller effect in SiV 0 defect of diamond [42,43].In the next section, we discuss the Jahn-Teller effect in detail for t (2) 2 electronic configuration in T d symmetry.The theory of Jahn-Teller interaction in T d symmetry was already discussed in the literature (see Ref. 44 and references therein).We apply this theory to the Ni 0 s defect.The T d symmetry is dynamically distorted by t and e vibration modes, i.e.T ⊗ (t ⊕ e) vibronic system.With assuming an effective t vibration and an effective e vibration, the vibronic system may read as where X, Ŷ , Ẑ operators depict the geometry distortion from the t vibration mode governed by ĤT JT while V , Ŵ depict to that of e vibration mode governed by ĤE JT .F E and F T electron-vibration coupling parameters can be determined [44] as where the E T,E JT is the Jahn-Teller energy and ℏω T,E is the respective effective phonon t and e frequencies.The global minima in the five-fold configurational space are either one of the 3 equivalent tetragonal extrema with D 2d symmetry or 4 equivalent trigonal extrema with C 3v symmetry.The actual global energy minima, whether tetragonal or trigonal, depend on the properties of the given system.As the t and e vibrations are orthogonal each other and they lead to trigonal and tetragonal distortions, there individual contributions can be determined by DFT calculations as where E DFT tot (Γ) total energies are obtained from HSE06 DFT calculations with geometry optimization under the Γ point group.We determined ℏω T,E by mapping the parabolicity of the adiabatic potential energy surface along the e and t normal coordinates.
We solved the JT problem up to 8-phonon limit by substituting our ab initio parameters E T JT = 160.5 meV, ℏω T = 45.0 meV, E E JT = 127.8meV, ℏω E = 69.5 meV into Eqs.( 17) and (16).In accordance to previous results [44], the lowermost eigenvalue of the electronvibration or vibronic system is a triple degenerate T level which is followed by an A level by 9.0 meV.| 3 T 1 ⟩ also experience spin-orbit effect as the maximum degeneracy in the T d double group is four-fold.In the next section, we study the spin-orbit interaction of this state.
2. Spin-orbit fine structure of Ni 0 s We introduce the spin-orbit operator as follows.The Lx,y,z operators depict an effective L = 1 orbital moment within the t 2 orbitals.One may note that there are two electrons in the t 2 orbital, thus where and Ŝx,y,z 's are regular spin operators for a given spin multiplicity.For the triplet spin of Ni 0 s , we use The complex i unit in Lx,y,z emphasizes that the t 2 orbitals lose their real-value character when spin-orbit interaction is considered unlike the case of electron-phonon interaction, e.g.Eq. ( 16)).For example, if the quantisation axis is along the "z" [001] axis then the eigenfunc- and |ml = 0⟩ = |t 2z ⟩ that we visualize in Fig. 1(b).
We determined the bare intrinsic spin-orbit strength at single particle level by ab initio calculations: .78 meV.The single particle l operator can be depicted by the same matrices that of Eq. ( 20) in contrast to a single electron's spin ŝ that can be depicted by two-times-two Pauli matrices: ŝx We note that λ (1) acts only between single particle levels but Ni 0 s is effectively a two-electron system, see Eq. (12).Therefore, the spin-orbit strength for the twoparticle system will be λ (1) ⟩.In order to solve the previous equation, we introduce an 1/2 Clebsch-Gordan coefficient for λ 0 = λ (1)   2 because only t ↑ 2z± posseses both angular and spin moments, so t ↑ 2z exhibiting ml = 0 remains intact from spin-orbit coupling.
However, the intrinsic value of spin-orbit coupling is reduced by the p T1 Ham reduction factor.That is, the L orbital operators are mediated by strong Jahn-Teller renormalization thus are subject to the so called Ham reduction factors (see section 5.6 about reduction factors in Ref. 44 for details).In other words, the triply degenerate t 2 phonons can entangle the three |t 2z± ⟩, |t 2z ⟩ orbital characters.That is the ground state that can be depicted in a Born-Oppenheimer basis by ladder operators acting on vibration modes: .., where we depicted only the first terms in the series expansion.Therefore the electron-phonon entangled | t2z+ ⟩ will exhibit a reduced spin-orbit strength as the ground state is not a pure |t 2z+ ⟩ electronic character anymore: |t 2z− ⟩ mixes in.The intrinsic orbital angular momentum will be partially quenched by the p T1 factor: Leff = p T1 ⟨ L⟩ phonons .Thus the observable spin-orbit coupling will be partially quenched too: λ eff = p T1 λ 0 .
We report the reduction factor as p T1 = 0.0349 by means of Eq. ( 16) when we capped the phonon quanta as n ≤ 9. Finally, the λ eff acting for | 4 T 1 ⟩'s L = 1 orbital and S = 1 spin moments will be where λ eff = 0.12 meV.Our positive λ eff is in agreement with Hund's rules: the lowermost state in a J = 0 (A 1 ) singlet followed by J = 1 (T 1 ) at λ eff energy and J = 2 (T 2 ⊕ E) at 3λ eff , where we depicted double group representations [45] for the T d symmetry in (. . . ) parentheses.However, higher order spin-orbit terms may dominate over the linear λ eff similarly to that of Ni − s 's excited state that we will describe in Eq. ( 29) in the next section.

E. Optically excited states of Ni − s
We turn to the interpretation of the optical spectra of Ni-related centres.First, we consider the optically allowed excited states of Ni − s .To this end, one has to go beyond t (3) 2 multiplets because no other quartet state can be formed than the 4 A 2 ground state.As can be seen in Fig. 1(a), a resonant e state occurs in the valence band from which an electron may be promoted to the empty t 2 state in the spin minority channel.This results in e (3) electron configuration which is equivalent to e (1) t (2) 2 hole configuration.We assume that the t 2 configuration reads as According to Hund's first rule, we expect that the highspin quartet states are more stable than the low-spin doublet states.We focus now on the discussion of the optically allowed quartet states.In the construction of the quartet states in Eq. ( 23) we may start with the doublet E state which can be expressed by (x 2 −y 2 ; 2z and It may be recognized that 4 T 1 states can be constructed by swapping |v⟩ and w⟩ in the multiplets of 4 T 2 states.Another important and surprising observation is that both 4 T 2 and 4 T 1 states may be described as a single Slater-determinant (see the last row in Eqs. ( 24) and ( 25)).As a consequence, the total energy of 4 T 2 and 4 T 1 states can be principally estimated from KS DFT ∆SCF method once the respective |xyv⟩ and |xyw⟩ electronic configurations could be converged in this procedure.We obtained 2.61 eV and 2.95 eV excitation energies w.r.t. the 4 A 2 ground state's energy within T d symmetry.
The electronic configurations and their energy levels may then be expressed by omitting the spin degrees of freedom as follows, where we define the P2 and P1 projectors as ) and ) However, we learnt from the 3 T 1 of Ni 0 s that is a subject to JT effect.Since 4 T 2 and 4 T 1 inherit this electron configuration one can suspect that JT effect is not negligible in 4 T 2 and 4 T 1 states of Ni − s .Indeed, ∼84 meV energy is released once we removed the T d symmetry constraint in our ab initio calculations during the ∆SCF procedure.In the next section, we the JT effect in the 4 T 2 and 4 T 1 states with using P2 and P1 projectors in Eqs.(27a) and (27b).

(T1 ⊕ T2) ⊗ (t ⊕ e) Jahn-Teller effect
We already analysed the JT effect for the T ⊗ (t ⊕ e) system in Sec.III D 1.We rely on this analysis which becomes more complex with taking both T 2 and T 1 states that might also interact with each other through electronphonon coupling that we describe by JT effect.This may be described as

+ ( ( ( ( ( ( ( ( ( ( h h h h h h h h h h
where we assume that ĤJT should not depend on |v⟩ and |w⟩ orbitals.That is, ĤJT in Eq. ( 28) is the dot product of ĤJT from Eq. ( 16) with a 2×2 unit matrix to span the six dimensional degrees of freedom from the orbitals.In the derivation, we used the lemma Î = P1 + P2 which leads to diagonal and offdiagonal terms where we can distribute the JT interaction acting on inside the 4 T 1,2 manifolds (diagonal terms) or induce orbital excitation  26).It can be seen that E T JT /4 minima at negative distortion is independent on the choice of Λ. Theoretically, it is possible to determine E T JT directly by ab initio calculations; however, we were not able obtain convergent results.Therefore, we assumed E T JT = 138 meV from Ni 2 s − charge state that exhibits the same electron occupancy in the t2 orbital.We note that yellow filled-in data points are ab inito results.(b) Polaronic eigenspectrum by means of Eq. ( 28) for Ni − s without the 3 2 ℏωT + 2 2 ℏωE zero point energy terms.(c) APES for E vibration modes towards the V distortion.Quantization axis and shape for t2x,y,z = x, y, z orbitals are shown in the first row of Fig. 1(b).
between the 4 T 1,2 multiplets.As a consequence, JT interaction induces mixing of electronic characters.
First, one may check that the "E" vibration modes do not lead to this kind of orbital mixing that we depicted in Eq. (28).That is, Eqs.(27a) and (27b) show that the projectors are diagonal at their bottom-right edges, i.e. 0 4 , 4 0 , respectively.This bottom-right edge is influenced only by the −F E V diagonal matrix element from Eq. ( 16) and thus does not inflict any interaction between A|xyv⟩ and A|xyw⟩ configurations of | 4 T 2 ⟩ and | 4 T 1 ⟩ multiplets, respectively.One may argue that the upper two 2×2 matrix blocks in Eq. ( 27) are not diagonal.However, these matrix blocks can be diagonalized with A|xy( w 2 w)⟩ orbital rotations.Therefore, offdiagonal terms in Eq. ( 28) do not appear for E-type distortions.Therefore, APES remains the same for | 4 T 1,2 ⟩ states.They are just shifted by Ξ electronic splitting of Eq. ( 26) and thus Eq. (18b) can be used to determine F E strength.According to our ab initio calculations, the JT distortion energy is E E JT = 84.2meV mediated by ℏω E = 68.7 meV vibrations.
However, F T in Eq. ( 18a) is effective in this regard.X, Ŷ , Ẑ in Eq. ( 16) are purely offdiagonals that can create cross-talks between 2×2 blocks of the projectors in Eq. ( 27) and thus mix the | 4 T 1,2 ⟩ electronic characters since the three 2×2 blocks are not diagonal in the same basis.As we saw in the previous paragraph, the third block is diagonal in {v; w} basis but the first two blocks are diagonal only in 2 v}, thus ĤT JT will cause orbital excitation between | 4 T 1,2 ⟩ multiplets.In other words, the APES will not be just a shifted replicas that of Ni 0 s 's T ⊗ (t ⊕ e) JT APES depicted in Sec.III D 1. Indeed, we report that the dis-tortion energy is severely quenched for "T "-type distortions: E T JT = 11.6 meV (that was 160 meV for Ni 0 s ) that demonstrates that Eq. ( 17) should not be directly employed for Ni − s 's excited states.Therefore we used E T JT = 138.4meV and ℏω T = 46.8meV that of Ni 2− s that resembles this excited state.Both configurations confine four electrons in t 2 orbitals thus we expect similar JT instability for these two cases approximately.We also assume that the pure atomic d orbitals in e states do not affect the Jahn-Teller instability too much.Therefore, we take an exact diagonalisation of Eq. ( 28) in order to approximate the vibronic levels up to n ≤ 8 phonon limit.By considering the offdiagonal terms, the | 4 T 2 ⟩ and | 4 T 1 ⟩ can be entangled by phonons which will ultimately lead to a vibronically strongly coupled lowermost | 4 T2 ⟩ polaronic state quickly followed by an additional | 4 T1 ⟩ state at Λ eff = 18.0 meV.We observe that at any choice for Λ below ≤ 342 meV the effective energy spacing will be strongly quenched as Λ eff = p E Λ by the p E Ham reduction factor.We note that p E affects orbital operators transforming as E representation of T d symmetry (see section 5.6 about reduction factors in Ref. 44 for details).In our present case, Λ energy exhibits the same order of magnitude that of E T JT or E E JT thus the two effects will compete with each other.As a consequence, p E will be dependent on Λ thus it is p E (Λ = 342 meV) = 18.02 meV 342 meV = 0.0527 in our case.

Spin-orbit fine structure
We learnt from 3 T 1 ground state of Ni 0 s that spin-orbit interaction occurs.Nazaré et al. was able to fit the following spin-orbit Hamiltonian (see Eq. (2) in Ref. 22) to the fine structure of the 2.51-eV absorption centre, where the L orbital operators are that of Eq. ( 20) and Ŝ spin operators are extended for quartet (S = 32 ) multiplicity to depict the spin-orbit interaction similar way to that of Eq. ( 19) for Ni 0 s .However, in contrast to that of Eq. ( 19), Eq. ( 29) is an effective Hamiltonian in which the L orbital operators are mediated by strong Jahn-Teller renormalization, i.e. those are subject to the so-called Ham reduction factors.In this case, the reduction factor is p T1 because our three L orbital operators transform as the T 1 representation.Therefore, the intrinsic orbital angular momentum will be partially quenched by the p T1 factor: Leff = p T1 ⟨ L⟩ phonons thus the observable spin-orbit coupling will be partially quenched too, i.e. λ eff = p T1 λ 0 .We report the reduction factor as p T1 = 0.0604 by means of Eq. ( 28) with using n ≤ 7 for the expression of polaronic states.
The spin-orbit coupling parameter should be also determined.We were not able to converge the excited state together with spin-orbit coupling, thus we determined λ 0 by taking the value of Ni2− s that bears the same number of electrons on its t 2 orbital: 18 meV.We note that λ (1) , strictly speaking, acts only on single particle KS levels.
However, our current level of theory was not able to determine the additional higher order κ and ρ parameters in Eq. ( 29)).That is, determining those coefficients requires second order effects that would utilize 2 ⟩⟨e ↑ | flipping terms.We report λ ± ≈ 50 meV from ab initio calculations that would scale the ϱ, κ parameters by means of second order perturbation theory as with assuming a heuristic ∼1 eV electronic splitting.Our result in Eq. ( 30) is in the same order of magnitude that was reported in Ref. 22 as κ = −0.326meV, ρ = +0.532meV further supporting the association of 2.51-eV absorption centre to the 4 A 2 → 4 T 1 optical transition of Ni − s .

F. Excited states of Ni 0 s
We continue the investigation with the excited states of e (1) 2 electronic configuration which could be either t ↑↑↑ 2 or t ↑↑↓ 2 spin configurations coupled to 2 E which is a hole left on the d orbital with e symmetry with either ↑ or ↓ spin state.We start the analysis with the t ↑↑↑ 2 case which reads as The | 3 E⟩ is the case when the S = 3/2 spin from | 4 A 2 ⟩ and the S = 1/2 spin from the |e⟩ hole is coupled by the opposite spin states.The maximally spin-polarised wavefunctions are the followings, where the |ξ⟩ can be either |v⟩ or |w⟩.We note that the A|x ↑ y ↑ z ↑ ξ ↓ ⟩ single Slater-determinant dominates (75%) the | 3 E⟩ ms=+1 wavefunction, thus we expect that its total energy can be well estimated by ∆SCF DFT method which results in 2.61 eV excitation energy within T d symmetry.The JT energy is negligibly small (3 meV).The total energy of the 5 E state should be accurate by ∆SCF DFT method as a single Slater-determinant, and it yields 2.08 eV excitation energy.The JT energy is again tiny (5 meV) as expected.
We now turn to the other manifolds based on t ↑↑↓ 2 electronic configurations which result in 2 E, 2 T 1 and 2 T 2 multiplets.We already proved during the discussion of Ni − s that 2 E multiplet exhibits the lowest total energy.Therefore, we only focus to that manifold as it is combined with a hole left on the d orbital with e symmetry which reads as follows where the other spin states can be readily obtained and explicitly given here.It can be observed that these states are highly correlated.The largest determinant contribution is 2 √ 12 2 = 33% prefactors to where the ∆SCF DFT results can converge to As a consequence, the total energies of these individual triplet states cannot be well determined by ∆SCF calculations.Nevertheless, we calculated the ∆SCF total energy of Eq. (35).We obtain 2.45 eV excitation energy within T d symmetry.The calculated JT energy is 43 meV which is expected due to the partially filled t 2 orbitals.
As the direct calculation of these states is not feasible with DFT methods, we do not analyze these states further but we give a very rough estimate for their ZPL energies at around 2.4 eV.

G. Ground state of Ni + s
We briefly discuss paramagnetic Ni + s which may be observed in special p-type doped diamond sample.Ni + s exhibits t (1) 2 electronic configuration which yields | 2 T 2 ⟩ multiplet.This is a subject to JT effect and spin-orbit interaction.
We can apply JT theory as given in Eqs. ( 16) and ( 17) which results in E T JT = 177.4meV, ℏω T = 46.8meV and E E JT = 115.2meV, ℏω E = 66.7 meV.These results imply that the Ham reduction factor is significant that may act on the spin-orbit energy gaps.We determined p T1 = 0.0209 up to the 8-phonon limit.We determined the spin-orbit coupling λ 0 as 8.6 meV.We note here that no Clebsch-Gordan coefficient is needed in this case since the | 2 T 2 ⟩ multiplet can be described as a single t ↑/↓ 2(±/0) ⟩ orbital.Therefore, the spin-orbit coupling will lead to an J = 1/2 followed by an high moment J = 3/2 state with λ eff = p T1 λ 0 = 0.179 meV.We again note that higher order spin-orbit terms may appear similarly to that of Eq. (29).

IV. DISCUSSION
We summarise the calculated excitation energies and photo-ionisation threshold energies of Ni s defects in diamond in Fig. 4. As we can see the neutral and negatively charged Ni s has photo-ionisation energies and excitation energies are close to 2.5-2.7 eV, thus one has to study the nature of optical excitations in great detail.This includes the vibronic states and levels as well as the spin-orbit couplings and fine structure levels (see Fig. 4).These results are the base for the discussion of the origin of 2.51-eV absorption centre and 2.56-eV ODMR centre.We note that the ∼2.5 eV and ∼3.0 eV threshold energies in the photo-EPR measurements [26] as the de-activation and re-activation of the W8 EPR centre associated with the ground state of Ni − s cannot be explained by simple ionisation and re-ionisation processes of Ni − s , according to our results.In particular, we show below multiple evidences that the 2.51-eV absorption centre is associated with Ni − s , thus the photo-ionisation threshold energy towards the conduction band minimum cannot be smaller than 2.51 eV.It is rather likely that carriers trapping is also involved in the process where the carriers might be generated by photo-ionisation of other defects in that diamond sample.
A. The 2.51-eV absorption centre and Ni − s We reiterate here that the 2.51-eV absorption centre is associated with a | 4 A 2 ⟩ → | 4 T 2 ⟩ optical transition according to the last observations which shows a replica at 16.5 meV above the ZPL feature [22].Both the ZPL and the replica show a fine structure with about ∼1.5 meV splitting where the ZPL's fine structure is depicted by Eq. ( 29)'s effective spin-orbital Hamiltonian.The Ni 0 s has T 1 orbital ground state whereas Ni − s has the appropriate A 2 ground state, thus we analyse the excited state of Ni − s as the possible candidate for the 2.51-eV absorption center.Indeed, the calculated ZPL energy is at 2.53 eV, close to the experimental one.Additionally, we were able to predict the experimentally observed SOC parameter λ expt = −0.163meV as λ eff = −0.145meV suggesting that the model in Ref. 22 is correct.
We found that the lowest energy quartet excited state is indeed | 4 T 2 ⟩ for Ni − s which is optically forbidden in the first order, see Fig. 3(b).The next quartet state is a 4 | T 1 ⟩ separated by about 18.02 meV from the | 4 T 2 ⟩ state that is very close to the observed replica at 16.5 meV.We emphasize that this effect is very similar that we found for the neutral silicon-vacancy center in diamond [42] where two polaronically entangled states lie close to each other in the Jahn-Teller solution.An additional analog to the properties of the neutral silicon-vacancy optical centre in diamond is that the first state is an optically for- FIG. 4. Ground and excited states of Nis defect in diamond in various charge states.The estimated zero-phonon-line energies of neutral excitation and photo-ionisation threshold energies are also given.bidden 3 | A 2u ⟩ which is quickly followed by the optically allowed 3 | E u ⟩ by ∼ 7 meV that is an order of magnitude smaller than that of the vibration mode's energy (ℏω = 76 eV).A similar situation is apparent for Ni − s : Λ expt = 16.8 meV < 47 meV = ℏω T .Therefore, the vibration replica is significantly broader (∼ 5 meV, see Fig. 1 in Ref. 22) than that of the ZPL (∼ 2 meV) but not as broad as that of a real first phonon vibronic replica (∼ 30 meV) of Huang-Rhys theory [33,46].Therefore, we argue that the 16.8-meV splitting is a fingerprint of the (T 1 ⊕ T 2 ) ⊗ (t ⊕ e) Jahn-Teller instability.
The optical activity of the dark | 4 T 2 ⟩ state may be explained by the JT coupling to the bright | 4 T 1 ⟩ multiplet which brings a | 4 T 1 ⟩ character to the | 4 T2 ⟩ state of about 5% due to offdiagonal terms in Eq. ( 28).The remaining 95% of | 4 T 1 ⟩'s optical transition dipole moment would be still active around Ξ + Λ = 2.51 eV + 0.34 eV = 2.85 eV (435 nm).However, this peak may be not recognized because it is simply obscured by the phonon sideband of the 2.51-eV peak and overlapping with the photo-ionisation threshold at about 2.72 eV converting Ni − s to Ni 0 s .Unfortunately, absorption data beyond > 2.54 eV (< 488 nm) is not reported for the 2.51-eV centre, to our best knowledge, to gain insight about this issue.
Next, we discuss the strength of the ZPL and the replica.First, we note that even the ZPL of the lowermost | 4 T 2 ⟩ only contains 7% | 4 T 1 ⟩ contribution and 1% Debye-Waller factor, thus only 1% optical transition dipole moment is visible in the replica of 2.51 eV.Additionally, the λ ± |t ↑ 2 ⟩⟨e ↓ | spin-orbit coupling may turn The missing fluorescence of the 2.51-eV band may be explained by the fact that the optical transition dipole moment is weak and | 2 T 2 ⟩, | 2 T 1 ⟩ and | 2 E⟩ doublets exist between the 2.51-eV optical transition that could lead to an efficient non-radiative channel for the 4 T 2 excited state.The estimated spin-orbit coupling between the 2 E and 4 A 2 states is about λ ± ∼ 50 meV which mediates the intersystem crossing in the non-radiative process.
B. The 2.56-eV ODMR centre and Ni 0 s The 2.56-eV ODMR centre is associated with the ground state of Ni − s but observed in luminescence which may be explained by the optical emission from Ni 0 s as speculated in Ref. 24.We discuss this scenario in the light of our results.One of the 3 E excited states could be calculated with relatively high accuracy at 2.61 eV which is resonant with the calculated photo-ionisation threshold energy at 2.62 eV.The excitation energy of the other three triplet excited states could not well calculated by ∆SCF method and it might fall below the photo-ionisation threshold.However, the strength of optical transition associated with |t 2 ⟩ → |e⟩ at Kohn-Sham orbital level is very weak because the t 2 state is well localised on Ni 3d orbitals.We think that it is unlikely that those states are involved in the ODMR signal.We here sketch another scenario.As we discussed in Sec.III A, Ni s spends in the negative charge state upon ultraviolet illumination according to our calculations, and ultraviolet illumination also releases holes in diamond samples [24].Ni − s can trap holes with very high hole capture rates that results in bound exciton states of Ni 0 s .We note that ODMR signal via bound exciton excited states have been observed for the neutral silicon-vacancy defect [47], and it was argued for nitrogen-vacancy centre that bound exciton states have giant hole capture cross section [48].We analysed the bound exciton states for silicon-vacancy centre in D 3d symmetry [47].For Ni − s plus hole system, we should consider T d symmetry as the Ni − s ground state electronic configuration is stable in T d symmetry and the hole split from valence band maximum should have negligible JT effect.For the hole state, we may use the effective mass theory as explained in Ref. 47.The lowest energy state transforms by A 1 symmetry which has an 1s type envelope function.The 1s effective mass state is relatively well localised around the core of the defect for which the optical transition dipole moment is relatively strong.The phonon sideband of the fluorescence spectrum may be well estimated by the optimized geometries of the Ni − s representing the excited state and Ni 0 s for the ground state within Franck-Condon approximation.Indeed, the calculated and observed PL spectra agree well (see Fig. 5).We conclude that the 2.56-eV PL and ODMR signals are associated with the bound exciton emission of Ni 0 s .Understanding the ODMR contrast requires the analysis of the fine structure of the ground state and the excited state.The ground state is 3 T 1 which splits to J = 0 (A 1 ) state, J = 1 (T 1 ) state and J = 2 (T 2 ⊕ E) states in ascending energy order where we label here the double group representations.The splitting is caused by the spin-orbit interaction which is roughly estimated to be at 0.12 meV based on our ab initio calculations that do not contain quadratic spin-orbit contributions that might be also significant.Nevertheless, we assume that the order of these levels retains.In this case, the ground state has no effective coupling to the external magnetic field in the ground state, thus it may be not observable in EPR measurements.The Jahn-Teller nature of the ground state could also make challenging to observe the triplet state in EPR experiments.In the excited state, the 4 A 2 electronic configuration has S 1 = 3/2 spin state which is combined with a hole in the valence band maximum with S 2 = 1/2 state (T 2 ).The composed system transforms at T 1 which again results in S = 0 A 1 state, S = 1 T 1 state and S = 2 T 2 ⊕E states split by spin-orbit interaction, where S = S 1 + S 2 here.The initial population of the excited state branch may depend on the applied external magnetic field, temperature and strain in the sample.The final strength and sign of the ODMR contrast will depend on the population of these states in the excited state manifold.Further analysis of this process and the g-factor of the state is out of the scope of the present study.

C. Interaction of nickel defects with other defects in diamond
In the previous sections we already outlined that defect-defect interactions play a crucial role in the manifestation of the ODMR signal associated with the substitutional nickel defect in diamond.Since the ionisation energies of defects are much larger than the operation temperatures (typically, room temperature) the exchange of carriers are mediated by photo-excitation.The nitrogen donor can be efficiently activated by photo-excitation energy at 2.1 eV and above.The optical signals of substitional nickel defects have ZPL energies at 2.51 eV and above, thus surrounding nitrogen donors around the nickel defects provide electrons towards the nickel defects upon illumination which may stabilise their negative charge state.On the other hand, if nickel-vacancy (NiV) defects are also present in the sample then this illumination results in ejecting holes by converting the 1.4-eV PL centre (negatively charged NiV defect) to doubly negative charged NiV defects (see Ref. [3]).Thus, the ratio and location of the nitrogen and NiV defects with respect to the substitutional nickel defects can be decisive about the magneto-optical fingerprints of the substitional nickel defects.If vacancies or vacancy aggregates are also present in diamond, e.g. after implantion, then they also enter as traps for holes (in their negative charge states) and as sources of holes upon illumination of neutral vacancies or vacancy clusters.The latter can indeed occur as the calculated acceptor level of the vacancy or divacancy with respect to VBM is lower than 2.5 eV [49].These defects and possibly boron acceptors are the major players in the photodynamics of the substitional nickel defect in diamond.

V. SUMMARY
We studied the electronic structure of the Ni 3d transition metal substitutional defect with T d symmetry by means of group theory analysis and plane wave supercell density functional theory calculations.We find that the negative and neutral charge states are the most relevant configurations for the substitutional Ni defect in diamond.We observed strong electron-phonon coupling and spin-orbit coupling in certain states of the defect.We compared the ab initio results with the previously reported 2.51-eV absorption centre and the 2.56-eV ODMR centre in diamond, and we associate these centres to the negatively charged defect and emission from the bound exciton excited state of the neutral defect in diamond, respectively.

FIG. 1 .
FIG.1.(a) Single particle levels and charge transition levels for Nis.Note that the spin-polarisation in the applied method splits the occupied (▲,▼) and unoccupied (△,▽) defect levels, and the Fock-exchange potential in the hybrid DFT further lifts the degeneracy within the same spin channel for the occupied and unoccupied orbitals.Additionally, we show the lowest possible optical excitation between t2 and e orbitals for Ni 0 s and Ni − s .We note that the e orbitals Ni 0 s are smeared into multiple valence band states, thus its position is only approximate.However, upon ∆SCF procedure the e hole is getting localized and elevated into the band gap.(b) Kohn-Sham orbitals of Ni − s in diamond.In the real wavefunctions, the red and cyan isovalent surfaces correspond to positive and negative values, respectively.We plot the |t2x,2y,2z⟩ orbitals from Ni − s 's triply degenerate (▲ ▲ ▲) occupied t2 levels.The |eu,w⟩ are taken from Ni 2− s 's nearly doubly degenerate (▼ ▼) occupied e orbitals.

FIG. 2 . 3 2
FIG. 2. (a) Lowest layer of APES for T vibration mode showing the four equivalent minima.(b) Geometry distortions for (i) E and (ii) T type modes towards configuration coordinates V and (X + Y + Z)/ (3) respectively.(c) APES for E vibration mode depicting the three equivalent minima points towards directions +V and − 1 2 V ± √ 3 2 W .(d) APES for T vibration mode towards the [111] distortion.Quantization axis and shape for |t [111] 2... ⟩ orbitals are shown in the second row of Fig. 1(b).(e) Polaronic eigenspectrum by means of Eq. (16) for Ni 0 s without the 3 2 ℏωT + 2 2 ℏωE zero point energy terms.(f) APES for E vibration mode towards the V distortion.Quantization axis and shape for |t2...⟩ orbitals are shown in the first row of Fig. 1(b).We note that yellow filled-in data points depict ab inito results.

FIG. 3 .
FIG. 3. (a) APES for T vibration modes towards the [111] distortion for Ni −s 's excited state.The dotted lines represent the APES where Λ = 0 is assumed in Eq.(26).It can be seen that E T JT /4 minima at negative distortion is independent on the choice of Λ. Theoretically, it is possible to determine E T JT directly by ab initio calculations; however, we were not able obtain convergent results.Therefore, we assumed E T JT = 138 meV from Ni 2 s − charge state that exhibits the same electron occupancy in the t2 orbital.We note that yellow filled-in data points are ab inito results.(b) Polaronic eigenspectrum by means of Eq. (28) for Ni − s without the 3 2 ℏωT + 2 2 ℏωE zero point energy terms.(c) APES for E vibration modes towards the V distortion.Quantization axis and shape for t2x,y,z = x, y, z orbitals are shown in the first row of Fig. 1(b).