The Group-IV-Vacancy Color Center in Diamond

TL;DR

This paper reviews and refines theoretical models of Group-IV vacancy centers in diamond, predicting their properties and guiding experimental efforts for quantum information applications.

Contribution

It provides a comprehensive theoretical framework for XV centers, integrating intrinsic and external interactions, and offers predictions to support future experimental research.

Findings

Refined models explain experimental data.

Predicted properties of XV centers.

Guidance for quantum information applications.

Abstract

Group-IV vacancy (G4V, or XV, where X = Si, Ge, Sn, Pb) color centers constitute a novel and promising class of defects in diamonds. This chapter reviews and refines the theoretical models for the XV systems, encompassing the intrinsic interactions, including spin-orbit coupling and electron-phonon interactions, and the external interactions involving strain, electric, light, and magnetic fields. Based on the refined model, we predict their properties, explain the experimental data, and suggest follow-up experiments. This article established a solid foundation for controlling the XV system, thus paving the way for quantum information processing.

Tables

Table 1. Table 1: Parameters for C–C and X–C bonds. The C–C bond length ( d 0 , C d_{\rm 0,C} ) and dissociation energy ( E 0 , C E_{\rm 0,C} ) are taken from Ref.  [ 25 ] , as are the X–C bond parameters ( d 0 , X d_{\rm 0,X} , E 0 , X E_{\rm 0,X} ). The potential-width parameter α C \alpha_{\rm C} is determined from both experiment  [ 26 ] and DFT calculations  [ 27 ] , while α X \alpha_{\rm X} values are estimated via Eq.  10 . Note: These values correspond to molecular systems and should be interpreted accordingly in crystalline contexts.

Bond	$d_{0,C/X}$ (Å)	$E_{0,C/X}$ (kJ/mol)	${\alpha }_{C/X}$ (Å-1)
C–C	1.54	346	2.2
Si–C	1.85	318	1.7
Ge–C	1.95	238	1.2
Sn–C	2.16	192	0.87
Pb–C	2.30	130	0.55

Table 2. Table 2: Dimensions of the triangular prism cantilever and electrodes, extracted from Fig.  17 .

Dimension	$x$	$l$	$y$	$z$	$h$
Size ($\mu$m)	29.0	22.6	2.0	1.6	2.4

Table 3. Table 3: Comparison of FEM results for | ϵ x ​ x M − ϵ y ​ y M | |\epsilon_{xx}^{\rm M}-\epsilon_{yy}^{\rm M}| extracted from Ref.  [ 14 ] and results for | ϵ x ​ x R − ϵ y ​ y R | |\epsilon_{xx}^{\rm R}-\epsilon_{yy}^{\rm R}| calculated using Eq. ( 110 ) with dimensions from Table  2 as inputs.

Voltage (V)	25	50	75	100	125	150
$\left|{ϵ}_{xx}^M-{ϵ}_{yy}^M\right|$ (${10}^{-4}$)	0.0200	0.0500	0.110	0.200	0.300	0.430
$\left|{ϵ}_{xx}^R-{ϵ}_{yy}^R\right|$ (${10}^{-4}$)	0.0154	0.0615	0.138	0.246	0.384	0.553

Table 4. Table 4: Inner products between the normalized vectors 𝐕 j \mathbf{V}_{j} ( j ∈ { 0 , x , y } j\in\{0,x,y\} ) and vectors 𝐃 k / ( a 0 ​ ϵ k ) \mathbf{D}_{k}/(a_{0}\epsilon_{k}) ( k ∈ { z ​ z , x ​ x , y ​ y , x ​ y , y ​ x , x ​ z , z ​ x , y ​ z , z ​ y } k\in\{zz,xx,yy,xy,yx,xz,zx,yz,zy\} ).

$\frac{{𝐕}_j\cdot {𝐃}_k}{a_0{ϵ}_k}$	${𝐃}_{zz}$	${𝐃}_{xx}$	${𝐃}_{yy}$	${𝐃}_{xy}$	${𝐃}_{yx}$	${𝐃}_{xz}$	${𝐃}_{zx}$	${𝐃}_{yz}$	${𝐃}_{zy}$
${𝐕}_{0+}$	$-\frac{25}{4\sqrt{114}}$	$-\frac{4}{\sqrt{114}}$	$-\frac{4}{\sqrt{114}}$	0	0	0	0	0	0
${𝐕}_{x+}$	0	$-\frac{2}{\sqrt{57}}$	$\frac{2}{\sqrt{57}}$	0	0	$-\frac{5}{\sqrt{114}}$	$-\frac{5}{\sqrt{114}}$	0	0
${𝐕}_{y+}$	0	0	0	$\frac{2}{\sqrt{57}}$	$\frac{2}{\sqrt{57}}$	0	0	$-\frac{5}{\sqrt{114}}$	$-\frac{5}{\sqrt{114}}$

Table 5. Table 5: Relative energy Δ ​ E \Delta E and (unnormalized) Boltzmann distribution n ​ ( Δ ​ E ) n(\Delta E) of the population in the ground (a) and excited manifolds (b). Calculations are based on the SnV parameters obtained from ab initio calculations  [ 13 ] , as listed in Table  6 . The notation ” A A e B B ” denotes number A × 10 B A\times 10^{B} .

(a)	$|1⟩$	$|3⟩$	$|5⟩$	$|7⟩$	$|9⟩$	$|11⟩$	$|13⟩$	$|15⟩$
$\Delta E$ (meV)	0	3.82	57.9	60.9	100	101	120	123
$n$ (5 K)	1	1e-4	0	0	0	0	0	0
$n$ (50 K)	1	0.412	0	0	0	0	0	0
$n$ (100 K)	1	0.642	1.2e-3	9e-4	0	0	0	0

Table 6. Table 6: XV system parameters ℏ ​ ω \hbar\omega , E JT E_{\text{JT}} , and λ \lambda reported in DFT calculations  [ 13 ] . The suffix ”-G”(”-E”) denotes the ground and first-excited manifold. The variable Δ 5 , 4 \Delta_{5,4} represents the energy difference between the fifth and fourth lowest eigenenergies.

	$\hslash \omega$	$E_{\text{JT}}$	$\lambda$	${\Delta }_{5,4}$
	(meV)	(meV)	(meV)	(meV)
SiV-G	85.2	42.3	0.82	52.0
GeV-G	82.2	30.1	2.20	54.5
SnV-G	79.4	21.6	8.28	54.1
PbV-G	74.9	15.6	34.6	41.6
SiV-E	73.5	78.5	6.96	31.4
GeV-E	73.0	85.7	36.1	26.7
SnV-E	75.6	83.1	96.8	23.0
PbV-E	78.6	91.6	245	14.5

Table 7. Table 7: Precisions of the quench factors based on the original and refined definitions for various XV systems. The suffix ”-G”(”-E”) denotes the ground (first-excited) manifold. The notion A ​ e ​ B AeB means A × 10 B A\times 10^{B} .

	$p$	$p^{\prime }$	${ϵ}_{\mathrm{qch},z}$	${ϵ}_{\mathrm{qch},z}^{\prime }$	$q$	$q^{\prime }$	${ϵ}_{\mathrm{qch},x}$	${ϵ}_{\mathrm{qch},x}^{\prime }$
SiV-G	0.30	0.30	-0.029	-0.029	0.65	0.65	1.7e-4	1.5e-4
GeV-G	0.38	0.38	-0.029	-0.029	0.69	0.69	3.7e-4	2.1e-4
SnV-G	0.46	0.46	-0.089	-0.092	0.73	0.73	4.5e-3	7.0e-4
PbV-G	0.53	0.49	-0.24	-0.030	0.77	0.73	0.10	9.8e-3
SiV-E	0.12	0.12	-0.16	-0.15	0.56	0.56	1.1e-3	5.6e-4
GeV-E	0.11	0.11	-0.50	-0.50	0.55	0.55	1.9e-2	3.2e-3
SnV-E	0.12	0.12	-0.71	-0.71	0.56	0.53	0.15	2.4e-2
PbV-E	0.11	0.10	-0.84	-0.85	0.55	0.42	0.99	0.16

Table 8. Table 8: XV system parameter estimations from Secs.  2.5 and  2.6 .

	$\hslash \omega$	$F_0$	$F$	$\lambda$
	(meV)	(eV/Å)	(eV/Å)	(meV)
SiV	136	39${\alpha }_{\mathrm{Si}}$	27${\alpha }_{\mathrm{Si}}$	52${\beta }_{\mathrm{Si}}$
GeV	120	35${\alpha }_{\mathrm{Ge}}$	25${\alpha }_{\mathrm{Ge}}$	603${\beta }_{\mathrm{Ge}}$
SnV	114	28${\alpha }_{\mathrm{Sn}}$	20${\alpha }_{\mathrm{Sn}}$	1841${\beta }_{\mathrm{Sn}}$
PbV	111	25${\alpha }_{\mathrm{Pb}}$	18${\alpha }_{\mathrm{Pb}}$	7706${\beta }_{\mathrm{Pb}}$

Table 9. Table 9: Partition factors ( α X \alpha_{\rm X} ) and screening factors ( β X \beta_{\rm X} ) of the XV systems in their ground (-G) and first-excited (-E) manifolds. The calculation is based on Thiering’s DFT simulations  [ 13 ] as shown in Tab.  6 and relation ( 175 ).

	${\alpha }_X$	${\beta }_X$
SiV-G	0.049	0.016
GeV-G	0.043	0.0036
SnV-G	0.044	0.0045
PbV-G	0.039	0.0045
SiV-E	0.058	0.13
GeV-E	0.065	0.060
SnV-E	0.082	0.053
PbV-E	0.10	0.032

Table 10. Table 10: Lowest splittings in the ground and first-excited manifolds for all the XV species measured in experiments ( Δ 0 , X , mea g / e \Delta_{0,\rm X,mea}^{g/e} ) along with the ones calculated by Thiering et al. ( Δ 0 , X , Thi g / e \Delta_{0,\rm X,Thi}^{g/e} )  [ 13 ] . One can derive these values with the Hamiltonian ( 55 ) along with Thiering’s DFT simulated parameters as shown in Tab.  6 . Note that the unit is in GHz, and 1 meV = 242 GHz.

	${\Delta }_{0,X,\mathrm{mea}}^{g/e}$ (GHz)	${\Delta }_{0,X,\mathrm{Thi}}^{g/e}$ (GHz) [13]
SiV-G	50 [2]	61
GeV-G	170 [3]	207
SnV-G	850 [6]	946
PbV-G	5700 [10]	4514
SiV-E	260 [2]	215
GeV-E	1120 [3]	987
SnV-E	3000 [6]	2897
PbV-E	no data	7051

Table 11. Table 11: Strain susceptibility parameters from 1) T-conv: Thiering’s DFT calculation for the SiV and SnV system parameters  [ 13 ] (shown in Tab.  6 ) converted with relations ( 175 ) and ( 174 ); 2) M-res: restored parameters based on the Meesala’s experimental measurement for the quenched (See Sec.  7.2 ) SiV strain susceptibility  [ 14 ] ; Based on Thiering’s DFT calculation as shown in Tab.  7 , the quench factors q g ′ = 0.65 q_{\rm g}^{\prime}=0.65 , q e ′ = 0.56 q_{\rm e}^{\prime}=0.56 ; 3) G-PBE/G-SCAN: Guo’s strain susceptibility DFT calculation for the SnV strain susceptibilities with the PBE/SCAN functionals  [ 15 ] . The unit for the strain susceptibility parameters d d and f f is PHz/strain.

		Ground Manifold			first-excited Manifold
	data set	$d$	$f$	$d/f$	$d$	$f$	$d/f$
SiV	T-conv	0.30	1.07	0.28	0.36	1.26	0.28
	M-res	1.3/$q_g^{\prime }$	-1.7/$q_g^{\prime }$	-0.76	1.8/$q_e^{\prime }$	-3.4/$q_e^{\prime }$	-0.53
SnV	T-conv	0.20	0.71	0.28	0.38	1.33	0.28
	G-PBE	1.57	-1.12	-1.40	1.91	-5.11	-0.37
	G-SCAN	1.67	-1.12	-1.49	1.84	-5.18	-0.35

Table 12. Table 12: Strain susceptibility parameters from 1) T-conv: Thiering’s DFT calculation for the SiV system parameters  [ 13 ] (shown in Tab.  6 ) converted with relations ( 175 ), ( 174 ), and ( 173 ); 2) M-res: Meesala’s experimental measurement for the quenched (See Sec.  7.2 ) SiV strain susceptibility  [ 14 ] . Unlike d d or f f , no quench factors are associated with t ∥ t_{\parallel} and t ⟂ t_{\perp} ; hence, the restoration to the original values is exempted. The unit for the strain susceptibility parameters t ∥ t_{\parallel} and t ⟂ t_{\perp} is PHz/strain.

	data set	$t_{\parallel ,e}-t_{\parallel ,g}$	$t_{⟂,e}-t_{⟂,g}$
SiV	T-conv	-0.16	-0.11
	M-res	-1.7	0.078

Equations

$$\begin{array}{rl}|{\psi }_{+,s}⟩=& \frac{1}{\sqrt{2}}\left(|{\psi }_s^{\mathrm{upp}}⟩+|{\psi }_s^{\mathrm{low}}⟩\right),\\ |{\psi }_{-,s}⟩=& \frac{1}{\sqrt{2}}\left(|{\psi }_s^{\mathrm{upp}}⟩-|{\psi }_s^{\mathrm{low}}⟩\right).\end{array}$$

$$\begin{array}{rl}|{\psi }_0^{\mathrm{upp}}⟩& =\frac{1}{\sqrt{3}}{\left(1,1,1\right)}^{\text{T}},\\ |{\psi }_{+1}^{\mathrm{upp}}⟩& =\frac{1}{\sqrt{3}}{\left(1,{\text{e}}^{\text{i}2\pi /3},{\text{e}}^{\text{i}4\pi /3}\right)}^{\text{T}},\\ |{\psi }_{-1}^{\mathrm{upp}}⟩& =\frac{1}{\sqrt{3}}{\left(1,{\text{e}}^{-\text{i}2\pi /3},{\text{e}}^{-\text{i}4\pi /3}\right)}^{\text{T}}.\end{array}$$

$$\begin{array}{rl}\frac{1}{2}\left({\mathrm{e}}^{\mathrm{i}{\varphi }_{+}^{\mathrm{upp}}}+{\mathrm{e}}^{\mathrm{i}{\varphi }_{-}^{\mathrm{upp}}}\right)& =x,\\ \frac{-\mathrm{i}}{2}\left({\mathrm{e}}^{\mathrm{i}{\varphi }_{+}^{\mathrm{upp}}}-{\mathrm{e}}^{\mathrm{i}{\varphi }_{-}^{\mathrm{upp}}}\right)& =y.\end{array}$$

$$\begin{array}{rl}|{\psi }_x^{\mathrm{upp}}⟩& =\frac{1}{\sqrt{2}}\left(|{\psi }_{+1}^{\mathrm{upp}}⟩+|{\psi }_{-1}^{\mathrm{upp}}⟩\right)=\frac{1}{\sqrt{6}}{\left(2,-1,-1\right)}^{\text{T}},\\ |{\psi }_y^{\mathrm{upp}}⟩& =\frac{-\mathrm{i}}{\sqrt{2}}\left(|{\psi }_{+1}^{\mathrm{upp}}⟩-|{\psi }_{-1}^{\mathrm{upp}}⟩\right)=\frac{1}{\sqrt{2}}{\left(0,1,-1\right)}^{\text{T}}.\end{array}$$

$$\begin{array}{rl}|{\psi }_{+,0}⟩& =\frac{1}{\sqrt{6}}(1,1,1,1,1,1{)}^{\text{T}},\\ |{\psi }_{-,0}⟩& =\frac{1}{\sqrt{6}}(1,1,1,-1,-1,-1{)}^{\text{T}}.\end{array}$$

$$\begin{array}{rl}|{\psi }_{+,+1}⟩& =\frac{1}{\sqrt{6}}(1,{\text{e}}^{\text{i}2\pi /3},{\text{e}}^{\text{i}4\pi /3},1,{\text{e}}^{\text{i}2\pi /3},{\text{e}}^{\text{i}4\pi /3}{)}^{\text{T}},\\ |{\psi }_{+,-1}⟩& =\frac{1}{\sqrt{6}}(1,{\text{e}}^{-\text{i}2\pi /3},{\text{e}}^{-\text{i}4\pi /3},1,{\text{e}}^{-\text{i}2\pi /3},{\text{e}}^{-\text{i}4\pi /3}{)}^{\text{T}},\\ |{\psi }_{-,+1}⟩& =\frac{1}{\sqrt{6}}(1,{\text{e}}^{\text{i}2\pi /3},{\text{e}}^{\text{i}4\pi /3},-1,-{\text{e}}^{\text{i}2\pi /3},-{\text{e}}^{\text{i}4\pi /3}{)}^{\text{T}},\\ |{\psi }_{-,-1}⟩& =\frac{1}{\sqrt{6}}(1,{\text{e}}^{-\text{i}2\pi /3},{\text{e}}^{-\text{i}4\pi /3},-1,-{\text{e}}^{-\text{i}2\pi /3},-{\text{e}}^{-\text{i}4\pi /3}{)}^{\text{T}}.\end{array}$$

$$\begin{array}{rl}|{\psi }_{+,x}⟩& =\frac{1}{\sqrt{12}}(2,-1,-1,2,-1,-1{)}^{\text{T}},\\ |{\psi }_{+,y}⟩& =\frac{1}{2}(0,1,-1,0,1,-1{)}^{\text{T}},\\ |{\psi }_{-,x}⟩& =\frac{1}{\sqrt{12}}(2,-1,-1,-2,1,1{)}^{\text{T}},\\ |{\psi }_{-,y}⟩& =\frac{1}{2}(0,1,-1,0,-1,1{)}^{\text{T}}.\end{array}$$

$$\begin{array}{rl}|{\psi }_x⟩& =\frac{1}{\sqrt{2}}\left(|{\psi }_{+1}⟩+|{\psi }_{-1}⟩\right),\\ |{\psi }_y⟩& =\frac{-\text{i}}{\sqrt{2}}\left(|{\psi }_{+1}⟩-|{\psi }_{-1}⟩\right).\end{array}$$

$$V(r_{\mathrm{a}})=E_0{\left[1-{\mathrm{e}}^{-\alpha (r_{\mathrm{a}}-d_0)}\right]}^2-E_0,$$

$${\alpha }_{\mathrm{X}}\approx \frac{E_{0,\mathrm{X}}{d}_{0,\mathrm{C}}}{E_{0,\mathrm{C}}{d}_{0,\mathrm{X}}}{\alpha }_{\mathrm{C}},$$

$$\Delta D\ll 1/\alpha ,$$

$$V(r_{\mathrm{a}})\approx -E_0+\frac{1}{2}k(r_{\mathrm{a}}-d_0{)}^2,$$

$$k=2E_0{\alpha }^2.$$

$$K_{\mathrm{C},x/y}$$

$$K_{\mathrm{C},\mathrm{z}}$$

$$K_{\mathrm{X}}$$

$$\hslash {\omega }_{\mathrm{X}}=\{\begin{array}{ll}74.7\text{meV}& \text{(SiV)}\\ 28.5\text{meV}& \text{(GeV)}\\ 14.5\text{meV}& \text{(SnV)}\\ 5.8\text{meV}& \text{(PbV)}\end{array},$$

$$\hslash {\omega }_{\mathrm{C},x/y}$$

$$\hslash {\omega }_{\mathrm{C},\mathrm{z}}$$

$$\Delta D=\sqrt{\frac{\hslash }{2\mu \omega }},$$

$$\Delta D_{\mathrm{X}}=\{\begin{array}{ll}0.032\text{A˚}& \text{(SiV)}\\ 0.032\text{A˚}& \text{(GeV)}\\ 0.035\text{A˚}& \text{(SnV)}\\ 0.042\text{A˚}& \text{(PbV)}\end{array}.$$

$$\Delta D_{\mathrm{C},x/y}$$

$$\Delta D_{\mathrm{C},\mathrm{z}}$$

$$\Delta d=\frac{1}{2}{z}_i.$$

$$\frac{\Delta E}{\Delta d}\approx -2F_{d0}:=-\frac{8\alpha e^2}{\pi {ϵ}_0{d}_{0,\mathrm{X}}^2}.$$

$$H_{\mathrm{vo}}=-F_{d0}\sum _{i=1}^6{z}_i|{\sigma }_i⟩⟨{\sigma }_i|+H_{\mathrm{h}-\mathrm{x}},$$

$${\mathbf{z}}_i^{\mathrm{T}}{\mathbf{V}}_j:=-⟨{\sigma }_i||{\psi }_j⟩,$$

$$H_{\mathrm{vo}}=-F_{d0}\sum _j{Q}_j\mathrm{diag}({\mathbf{V}}_j)+H_{\mathrm{h}-\mathrm{x}},$$

$$Q_j=\sum _{i=1}^6{\mathbf{V}}_j[i]\cdot z_i,$$

$$H_{\text{vo}}=F_0{Q}_0{I}^{\text{o}}+F\left(Q_x{\sigma }_z^{\text{o}}-Q_y{\sigma }_x^{\text{o}}\right),$$

$$\begin{array}{rl}{I}^{\text{o}}& =|{\psi }_x⟩⟨{\psi }_x|+|{\psi }_y⟩⟨{\psi }_y|,\\ {\sigma }_z^{\text{o}}& =|{\psi }_x⟩⟨{\psi }_x|-|{\psi }_y⟩⟨{\psi }_y|,\\ {\sigma }_x^{\text{o}}& =|{\psi }_y⟩⟨{\psi }_x|+|{\psi }_x⟩⟨{\psi }_y|.\end{array}$$

$$\begin{array}{rl}{F}_0& =\frac{1}{\sqrt{3}}{F}_{d0}=\{\begin{array}{ll}39{\alpha }_{\mathrm{Si}}\text{eV/A˚}& \text{(SiV)},\\ 35{\alpha }_{\mathrm{Ge}}\text{eV/A˚}& \text{(GeV)},\\ 28{\alpha }_{\mathrm{Sn}}\text{eV/A˚}& \text{(SnV)},\\ 25{\alpha }_{\mathrm{Pb}}\text{eV/A˚}& \text{(PbV)},\end{array}\\ F& =\frac{1}{\sqrt{6}}{F}_{d0}=\{\begin{array}{ll}27{\alpha }_{\mathrm{Si}}\text{eV/A˚}& \text{(SiV)},\\ 25{\alpha }_{\mathrm{Ge}}\text{eV/A˚}& \text{(GeV)},\\ 20{\alpha }_{\mathrm{Sn}}\text{eV/A˚}& \text{(SnV)},\\ 18{\alpha }_{\mathrm{Pb}}\text{eV/A˚}& \text{(PbV)}.\end{array}\end{array}$$

Full text

1 Introduction

Group-IV vacancy (XV, where X = Si, Ge, Sn, Pb) color centers in diamond [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] have recently emerged as promising platforms for quantum information processing. Their appeal stems from their high zero-phonon line (ZPL) proportion—the spectral feature corresponding to phonon-free electronic transitions—as well as their robustness against charge noise and resistance to thermal excitation. These properties position XV centers as highly attractive candidates for applications in quantum technologies. The motivation for the work presented in this chapter is fourfold.

Firstly, despite the enormous potential of XV centers, these systems consist of complex interactions between the orbital and spin movements of a hole (the absence of an electron in a fully occupied system), phonons (vibrations of atomic nuclei), and their coupling with strain, as well as external electric and magnetic fields. These intricate interactions make understanding them difficult. Several models have been developed and used for numerical calculations and experimental data analysis [11, 12, 13, 14, 15]; however, as detailed in Sec. 8, mutual discrepancies and missing interaction terms occasionally appear in these models. This emphasizes the necessity for a structured and comprehensible framework to elucidate the definitions of established concepts and to identify critical prerequisite conditions for using a specific approximate model in the characterization of XV centers. Starting from the foundational level, this article tackles these challenges by offering precise and comprehensive explanations of essential concepts and methodologies. In particular, this article elucidates the interactions arising from the fundamental force between a solitary electron and the atomic level, offering a detailed physical depiction of the XV systems. Additionally, it avoids employing sophisticated mathematical techniques like group theory, but ultimately arrives at the same conclusions regarding the configurations of electronic orbital eigenstates. This approach is intended to promote an accurate comprehension of XV centers, thereby minimizing errors in the application of concepts and models.

Secondly, this article addresses two key issues found in the previous modeling of the XV systems. To begin with, it contends that phonon modes should not be classified without considering the electronic orbital eigenstate configurations, as was done in Hepp’s work [11]. Instead, they should be classified with these configurations taken into account, as explained in Sec. 2.5 and Sec. 8. This approach results in a simpler expression for the electron-phonon interaction Hamiltonian and validates the assumption of two-mode coupling employed in Thiering’s calculations [13]. Moreover, this article highlights the discrepancy between the model formulated by Hepp [11] and the one implemented in Thiering’s calculations [13] regarding the Jahn-Teller effect (Sec. 8). Through meticulous analysis, this article endorses the model utilized in Thiering’s calculations and introduces a refined model that builds upon it, achieving enhanced precision for the XV systems with higher spin-orbit coupling strengths (Sec. 7.2.2).

Thirdly, despite abundant experimental measurements, those data lack sufficient theoretical explanation. This article explores possible theoretical frameworks behind these observations and obtains four main achievements. Firstly, given the fact that the electron-phonon interaction and strain interaction arise from the same mechanism, the interaction between the hole and the nuclear displacements, we explored the connection between the two interactions as presented in the quantitative relationship between their interaction strengths(Sec. 4.5). Besides, we explained the magnitude of the Stark shift (Sec. 21) and the intensity branching ratio in the light emission/absorption spectra (Sec. 5.4). Finally, we anticipate that an unstrained XV system will preserve its degeneracy when subjected to a magnetic field perpendicular to the principal axis, a phenomenon that Kramer’s theorem does not account for (Sec. 3.3). Our research indicates that this degeneracy originates from a novel symmetry associated with a joint reflection operation acting on the phonon, hole orbital, and hole spin movements (Sec. 3.2.2).

Finally, with the relationship between electron-phonon interaction strength and strain susceptibilities now determined, this article identifies a discrepancy between the XV system parameters derived by Thiering et al. and the strain susceptibilities reported by Meesala et al. [14] as well as those calculated by Guo et al. [15]. These differences are observed in both signs and magnitude. This article offers possible explanations for these contradictions, attributing them to the intrinsic strain of the XV sample. Additionally, we introduce a technique for measuring the intrinsic strain in the XV sample (Sec. 9). This approach allows for mitigating any potential effects from the intrinsic strain.

2 Intrinsic interactions

2.1 XV Center Formation

A Group-IV vacancy (XV, where X=Si, Ge, Sn, Pb) center in diamond is a lattice defect formed by two stages. Firstly, two neighboring carbon atoms are removed and replaced with a single X (X=Si, Ge, Sn, Pb) atom positioned centrally between the two vacancies, as illustrated in Fig. 1. Secondly, an electron from a nearby donor is captured by the X atom and the six surrounding carbon atoms [16].

The X atom, together with six surrounding carbon atoms and the captured electron, constitutes an XV system—the central subject of this chapter. When two adjacent carbon atoms are removed from the diamond lattice, it leaves the six remaining neighboring carbon atoms each with an unsaturated covalent bond, commonly referred to as a dangling bond. Each dangling bond hosts an unpaired electron occupying an $sp^{3}$ hybrid orbital.

The introduction of the X atom, which has four valence electrons, alters this configuration. In the presence of two additional external electrons, the X atom could form covalent bonds with all six neighboring carbon atoms through its $sp^{3}d^{2}$ hybrid orbitals, fully saturating the bonding structure. However, in the actual XV center system, only a single external electron is added instead of two. As a result, one of the six X-carbon (X-C) covalent bonds remains incomplete, leaving behind a hole that exists in a quantum superposition over the six bonds. This delocalized hole is fundamental to the quantum properties of the XV center, as it enables the emergence of well-defined and controllable quantum states, which are central to its quantum applications. In the following, we will determine the eigenstates governing the hole’s orbital motion and derive the corresponding interaction Hamiltonians that describe its coupling to both nuclear vibrations (phonons) and the hole’s spin.

2.2 Electronic Orbital Eigenstates

In order to determine the orbital eigenstates of the hole, we must first establish the orbital eigenstates of the 11 valence band electrons in the XV center system. To achieve this, we begin by determining the orbital eigenstates of a single electron in the Coulomb potential formed by the X atomic core and its six surrounding carbon atomic cores. Here, an atomic core refers to an atom with its valence-shell electrons removed, leaving only the nucleus and the fully occupied inner electron shells.

Fig. 2(a) shows the six covalent bonds in the XV center system, each bond resembles a $\sigma$-bond [17] consisting of one $sp^{3}d^{2}$ hybrid orbital from the X atom and one $sp^{3}$ hybrid orbital from the nearby carbon atom as shown in Fig. 2(b). These six covalent bonds form a simplified set of basis states, nearly orthogonal and complete to present the eigenstate of the electron. We define the six covalent bonds as quantum states $\left\{\ket{\sigma_{i}}\right\}$ ($i=1,2,\dots,6$), ordered according to a Cartesian coordinate system as illustrated in Fig. 2(a). This coordinate system is established in the reference frame of the six covalent bonds, which are assumed to remain fixed relative to one another. Within this frame, the X atomic core and the six surrounding carbon atomic cores are free to vibrate around their equilibrium positions, which themselves are considered stationary.

The origin of the coordinate system is set at the equilibrium position of the X atomic core, effectively defining the center of the XV system. The $z$-axis is oriented along the line connecting the two vacant carbon sites, while the $x$-axis is chosen such that the $x$-$z$ plane contains two of the X–C bonds—one located on the $+z$ side and the other on the $-z$ side. The $y$-axis is then uniquely determined to ensure a right-handed coordinate system.

Based on this coordinate system, we establish the spatial arrangement of the six basis states. The states $\ket{\sigma_{1}}$, $\ket{\sigma_{2}}$, and $\ket{\sigma_{3}}$ reside on the $+z$ side, with $\ket{\sigma_{1}}$ lying in the $x$-$z$ plane, while $\ket{\sigma_{2}}$ and $\ket{\sigma_{3}}$ are positioned counterclockwise when viewed from the $-z$ direction. Similarly, the states $\ket{\sigma_{4}}$, $\ket{\sigma_{5}}$, and $\ket{\sigma_{6}}$ are symmetrically positioned on the $-z$ side, directly opposite to $\ket{\sigma_{1}}$, $\ket{\sigma_{2}}$, and $\ket{\sigma_{3}}$, respectively, as depicted in Fig. 2(a).

Although group theory provides a rigorous framework for determining the eigenstates of an electron [11], the standard procedure is extensive and difficult to master quickly. Instead, we adopt a simplified approach based on a single foundational theorem: the symmetry of the eigenstates reflects the symmetry of the potential in a quantum system (detail shown below). By decomposing the XV center into two simpler subsystems with well-understood symmetries, this method enables a more intuitive and accessible derivation of the electronic eigenstates within the Coulomb potential of the XV center, referred to here as the XV potential.

The XV potential exhibits inversion symmetry, indicating that the $-z$ side of the potential mirrors the $+z$ side under an inversion operation. Consequently, the eigenstates of the XV system can be categorized into two groups, $\{\ket{\psi_{+,s}}\}$ and $\{\ket{\psi_{-,s}}\}$, with the subscript $s$ denoting some indices. They take the following forms:

[TABLE]

Here, $\left\{\ket{\psi^{\rm upp}_{s}}\right\}$ and $\left\{\ket{\psi^{\rm low}_{s}}\right\}$ represent wavefunctions in the subspace spanned by $\left\{\ket{\sigma_{1}},\ket{\sigma_{2}},\ket{\sigma_{3}}\right\}$ and $\left\{\ket{\sigma_{4}},\ket{\sigma_{5}},\ket{\sigma_{6}}\right\}$, respectively. The states $\ket{\psi^{\rm upp}_{s}}$ and $\ket{\psi^{\rm low}_{s}}$ share the same coefficients according to the $\ket{\sigma_{i}}$–$\ket{\sigma_{i+3}}$ correspondence ($i=1,2,3$). Consequently, determining the eigenstates of the XV system is reduced to identifying the three eigenstates of either subspace.

Considering the subspace spanned by $\left\{\ket{\sigma_{1}},\ket{\sigma_{2}},\ket{\sigma_{3}}\right\}$, the $\frac{2\pi}{3}$ rotational symmetry of the $+z$ side of the XV system mandates that eigenstates must be invariant under cyclic permutation of the coefficients within the basis, except for a global phase factor. Only three distinct wavefunctions fulfill this criterion. They are:

[TABLE]

Notice that $\ket{\psi_{+1}^{\rm upp}}$ and $\ket{\psi_{-1}^{\rm upp}}$ correspond to clockwise and counterclockwise permutations, respectively, both acquiring the same global phase factor $\text{e}^{\text{i}2\pi/3}$. Due to the symmetry between these two rotational directions, these states possess the same energy eigenvalue and are thus degenerate.

As dictated by rotational symmetry, these states are closely related to the eigenstates of the $z$-component of the orbital angular momentum operator, $L^{\rm o}_{z}$, making them particularly convenient for analyzing spin-orbit coupling interactions, as we will explore in Sec. 2.6. However, since these states exhibit uniform populations across the three basis orbitals, they interact identically with nuclear vibrational modes, regardless of the spatial direction of motion. For problems involving such interactions, it is advantageous to work with eigenstates that exhibit spatially polarized populations.

Fortunately, since $\ket{\psi_{+1}^{\rm upp}}$ and $\ket{\psi_{-1}^{\rm upp}}$ are degenerate in energy, any linear combination of these states is also an eigenstate of the Hamiltonian describing the upper subsystem of the XV center. This allows us to construct an alternative set of eigenstates with spatially polarized charge distributions.

The X–C bond orbitals are approximately linear in the radial directions (see Fig. 2), and the phase factors $\{\phi_{\pm}^{\rm upp}\}=\left\{0,\pm\frac{2\pi}{3},\pm\frac{4\pi}{3}\right\}$ in $\ket{\psi_{\pm}^{\rm upp}}$ are spatially dependent, with $\phi^{\rm upp}=\arctan(y/x)$, where $x$ and $y$ are the Cartesian coordinates in the XV system. This spatial dependence motivates the construction of the following linear combinations:

[TABLE]

Using these relations, we define the spatially polarized eigenstates:

[TABLE]

These states can be interpreted as the state $\ket{\psi_{0}^{\rm upp}}$ modulated by the spatial functions $x$ and $y$, respectively, as illustrated in Fig. 3. This polarization renders them particularly useful for describing interactions with nuclear vibrations.

Substituting the states (2) or their spatially polarized alternatives (4) into Eq. (1), we obtain the six eigenstates of the XV center’s bare Hamiltonian. The first two take the form:

[TABLE]

The remaining four are given by:

[TABLE]

Alternatively, using spatially polarized bases, the four states can be rewritten as:

[TABLE]

For problems involving spin-orbit interactions, we employ the eigenstates in Eq. (6), as they are expressed in terms of eigenstates of the orbital angular momentum. For electron-phonon interactions, we use the spatially polarized states in Eq. (7), which facilitate analysis of interactions with nuclear vibrations.

2.3 Eigenenergy Ordering

Evidently, $\ket{\psi_{+,x}}$ and $\ket{\psi_{+,y}}$ are degenerate in energy, as are $\ket{\psi_{-,x}}$ and $\ket{\psi_{-,y}}$. Consequently, the XV bare Hamiltonian has only four distinct eigenenergies, denoted as $E_{+,0}$, $E_{-,0}$, $E_{+,1}$, and $E_{-,1}$, corresponding to four energetic manifolds: $\left\{\ket{\psi_{+,0}}\right\}$, $\left\{\ket{\psi_{-,0}}\right\}$, $\left\{\ket{\psi_{+,x}},\ket{\psi_{+,y}}\right\}$, and $\left\{\ket{\psi_{-,x}},\ket{\psi_{-,y}}\right\}$, respectively.

Notice that there are 11 outer-shell electrons in an XV system. They tend to occupy the six eigenstates in the order of their eigenenergies from low to high. Therefore, it is crucial to determine the ordering of these four eigenenergies.

In the case of a single electron, we can employ a heuristic method to determine the order of the four eigenenergies: eigenstates with more nodes tend to have higher eigenenergies. Here, a node refers to a region where the wavefunction—and hence the probability density—vanishes, such as a plane in three dimensions, a line in two, or a point in one. This rule originates from the uncertainty principle: the presence of nodes confines the spatial extent of the wavefunction, reducing position uncertainty and increasing momentum uncertainty, which typically raises the energy [18].

According to this principle, $\ket{\psi_{+,0}}$, having no nodes, is the ground state. The states $\ket{\psi_{-,0}}$, $\ket{\psi_{+,x}}$, and $\ket{\psi_{+,y}}$ each contain a single nodal plane: $z=0$, $x=0$, and $y=0$, respectively [see Eq. (1), Fig. 5, and Eq. (3)]. However, as shown in Fig. 4, since the XV center’s spatial extent is smaller along $x$ and $y$ than along $z$, reducing position uncertainty, $\ket{\psi_{-,0}}$ with node plane $z=0$ has lower energy corresponding to the first excited state, while $\ket{\psi_{+,x}}$ and $\ket{\psi_{+,y}}$ have node planes $x=0$ and $y=0$, respectively, have higher energies forming the degenerate second excited states. Finally, $\ket{\psi_{-,x}}$ and $\ket{\psi_{-,y}}$, each possessing two nodal planes ($z=0$ and either $x=0$ or $y=0$), correspond to the highest-energy degenerate states. This argument is initially introduced in the work of Hepp [11].

Nonetheless, this rule may not apply to systems with multiple electrons as Coulomb interaction and/or exchange interaction take effect. The smaller the spatial extent of confinement, the larger the interaction energy. In the large interaction energy case, when occupied by two electrons, the orbitals with more nodes may have lower energies as allowing for larger spacing between the two electrons. Therefore, a more accurate determination of the order of the eigenenergies relies on the DFT (Density Functional theory) simulation. Coincidentally, the DFT result [13, 19] aligns with the ordering of the single-electron case.

The presence/absence of the central ($z=0$) node affects the partition factor as introduced in Sec. 2.5 and the screening factor as introduced in Sec. 2.6. Therefore, by calculating these factors based on the DFT calculations, one can infer the node configurations. We present these calculations based on Thiering’s DFT simulation [13] in in Sec. 8.

2.4 Hole Orbital Eigenstate

With the eigenstates and energy ordering of the XV system established, we now determine the orbital states of the 11 free electrons. Due to the low electron density, Coulomb interactions have a negligible effect on the electronic orbital eigenstates. Following the Pauli exclusion principle, electrons sequentially occupy the six eigenstates from lowest to highest energy, with each level accommodating up to two electrons of opposite spin. Consequently, the lowest-energy configuration leaves a single hole in the $\{\ket{\psi_{-,x}},\ket{\psi_{-,y}}\}$ manifold. The first excited configuration arises when an electron is promoted from the $\{\ket{\psi_{+,x}},\ket{\psi_{+,y}}\}$ manifold to the $\{\ket{\psi_{-,x}},\ket{\psi_{-,y}}\}$ manifold, as illustrated in Fig. 6(a).

From the hole’s perspective, the $\{\ket{\psi_{-,x}},\ket{\psi_{-,y}}\}$ manifold constitutes the ground states, while $\{\ket{\psi_{+,x}},\ket{\psi_{+,y}}\}$ forms the first excited manifold. Higher-energy manifolds can be constructed similarly but lie beyond the scope of this thesis.

Including spin, each manifold is four-fold degenerate in the absence of a magnetic field, as shown in Fig. 6(b). The energy gap between the ground and first-excited manifolds is significantly larger than that of other interactions, such as electron-phonon, static electric and magnetic field, and strain interactions. Thus, without influences (e.g., light) of high enough energy, inter-manifold couplings are suppressed, allowing us to treat these manifolds as independent subspaces when analyzing additional interactions.

For convenience, when specifying the manifold is unnecessary, we may use $\ket{\psi_{x}}$ and $\ket{\psi_{y}}$, or equivalently $\ket{\psi_{+1}}$ and $\ket{\psi_{-1}}$, to denote states within either manifold. These eigenbases are related by

[TABLE]

Thus, the Pauli-$Z$ operator in the $\{\ket{\psi_{+1}},\ket{\psi_{-1}}\}$ basis corresponds to the Pauli-$Y$ operator in the $\{\ket{\psi_{x}},\ket{\psi_{y}}\}$ basis. Since these two sets of bases are convertible, in the following sections, when introducing various interaction Hamiltonians, we will adopt the convention of expressing all the interaction Hamiltonians in the $\{\ket{\psi_{x}},\ket{\psi_{y}}\}$ basis.

2.5 Electron-phonon interaction

2.5.1 Localized phonon

As studied in solid-state physics, phonons, the quantized vibrations of a crystal lattice, come in various types. In the context of group-IV vacancy color centers in diamond, the most relevant phonons are localized phonons—vibrations of nuclei near the defect that are largely decoupled from the rest of the crystal lattice.

In contrast, bulk phonons, which propagate throughout the entire crystal, are macroscopically numerous but have limited influence on the electronic (or hole) orbitals localized at the defect. This is because bulk phonons have long wavelengths, especially at low temperatures, and primarily induce a global translation of the defect rather than energy shifts.

Localized phonons arise due to the unique properties of the defect site. The nuclei near the defect have different masses and bond stiffness compared to those in the bulk crystal, causing them to vibrate at distinct frequencies. This frequency mismatch limits their coupling to the bulk phonons, effectively confining the vibrations to the defect region. Localized or quasi-localized photons have been theoretically predicted and/or experimentally observed in GeV [20], SiV [21], and nitrogen-vacancy (NV) color centers in diamond [22].

In this section, we assume that these vibrations are effectively localized, so we confine the involved nuclei to the scope of the central X nucleus and its surrounding six carbon nuclei.

2.5.2 nuclear Vibrations

To begin with, we choose the rest frame of the equilibrium positions of the seven nuclei as the reference frame for analyzing the nuclear movement. Furthermore, for the free movement of the seven nuclei, we neglect the hole’s influence and treat all bonds as saturated. Following the concept of localized photons, we model each nucleus in the XV system as moving independently within a potential created by the covalent bonds surrounding it.

For the six carbon nuclei, each is attached to one X–C bond on one side and three C–C bonds on the other. Fig. 1(b) shows the 3D bonding structure, while Fig. 7(a) details the bond orientations. Due to the distinct X–C bond, the potential profile along the $z_{i}$ direction differs from those along the perpendicular $x_{i}$ and $y_{i}$ directions.

For simplicity, we model the potential along $x_{i}$ or $y_{i}$ as a fictitious 1D structure where the nucleus is confined between two C–C bonds (one on each side), as shown in Fig. 7(b). Conversely, the potential along $z_{i}$ is modeled as a 1D structure with the carbon nucleus confined between an X–C bond on one side and a C–C bond on the other (Fig. 7(c)).

For the X nucleus, the six surrounding X–C bonds are approximately evenly spaced (Fig. 2(a)). We thus model its potential along all directions as a 1D structure with the X nucleus confined between two X–C bonds (Fig. 7(d)).

For each covalent bond, we model the bond potential energy $V(r_{\rm a})$ as a function of internuclear distance $r_{\rm a}$ using the Morse potential [23, 24]:

[TABLE]

where $E_{\rm 0}$ is the bond dissociation energy, $d_{\rm 0}$ the equilibrium bond length, and $\alpha$ the potential-width parameter, representing the potential’s depth, minimum position, and curvature, respectively.

The potential-width parameter $\alpha_{\rm C}$ for the C–C bond has been well-characterized through experimental measurements [26] and theoretical calculations [27]. For X–C bonds, we estimate $\alpha_{\rm X}$ using the empirical relation:

[TABLE]

where subscripts C and X denote C–C and X–C bonds, respectively.

The Morse potential has a characteristic concave shape that can be well approximated by a harmonic oscillator potential when the vibration amplitude satisfies

[TABLE]

where $\Delta D$ represents the typical nuclear displacement. The approximating harmonic oscillator potential is expressed as

[TABLE]

where the stiffness $k$ is given by

[TABLE]

Based on the simplified models in Fig. 7(b)–(d), where each carbon or X nucleus is confined between two covalent bonds, the effective potential can be approximated as a superposition of two Morse potentials. The resulting effective stiffnesses are:

[TABLE]

where $K_{\rm C,\,x/y}$, $K_{\rm C,\,z}$, and $K_{\rm X}$ represent the stiffnesses for carbon nuclei along the $x/y$ and $z$ directions, and for X nuclei, respectively.

The phonon energy of the X nucleus, derived from the quantum harmonic oscillator model, is:

[TABLE]

where $\omega_{\rm X}=\sqrt{K_{\rm X}/\mu_{\rm X}}$ is the vibrational angular frequency and $\mu_{\rm X}\in\{28.0,72.6,118.7,207.2\}\,\text{Da}$ (1 Da = $1.66\times 10^{-27}$ kg) is the mass of the X nucleus.

Similarly, the carbon nucleus phonon energies in the XV system are:

[TABLE]

with $\omega_{\rm{C},x/y}=\sqrt{K_{\rm{C},x/y}/\mu_{\rm C}}$ and $\omega_{\rm{C},z}=\sqrt{K_{\rm{C},z}/\mu_{\rm C}}$.

To verify the condition (11) for XV systems, we calculate the root-mean-square displacement of nuclei using the ground state variance of a 1D quantum harmonic oscillator:

[TABLE]

where $\mu$ is the reduced mass and $\omega$ the vibrational angular frequency.

For different species of the X nuclei, we obtain the following typical displacements:

[TABLE]

For carbon nuclei, the displacements are:

[TABLE]

Comparing these displacements with typical $\alpha$ values (2.2 Å*-1* for C–C and 0.55-1.7 Å*-1* for X–C bonds from Table 1), we find $1/\alpha$ ranges from 0.59 Å  to 1.8 Å. Since all calculated $\Delta D$ values are $\ll 1/\alpha$, the harmonic approximation is justified for all cases considered.

2.5.3 Interaction Hamiltonian

As discussed, by employing the picture of a hole, the system originally comprising 11 covalent electrons is simplified to one single hole within a structure of six fully filled X–C bonds (totaling 12 electrons). The electron-phonon interaction problem thus results in the interaction between a single hole and the seven nuclei within the XV system. In this context, the interaction species under investigation is the Coulomb interaction.

Within this model, the nearest-neighbor approximation is applicable: When a hole is located in the $i$-th X–C covalent bond orbital ($i=1,2,\dots,6$), it primarily interacts with (1) the central X nucleus and (2) the $i$-th carbon nucleus through the Coulomb potential. The other five carbon nuclei are disregarded because their longer distances make their influence negligible. The symmetry of the six X–C bonds permits the interaction analysis to be focused on a single X–C bond unit, as depicted in Fig. 2(b).

We adopt the rest frame defined by the equilibrium positions of the nuclei. In this frame, the nuclei are allowed to move, and the X–C bond can deform accordingly. To first-order approximation, the hole-nucleus interaction is considered separately for each nucleus. That is, when evaluating the interaction between the hole and the motion of the carbon nucleus, the X nucleus is assumed to remain fixed, and vice versa.

The interaction Hamiltonian between a hole and a nucleus is characterized by the potential energy of the hole-nucleus setup, which depends on the positions of the hole and the nuclei. However, similar to the idea of an electronic cloud, the hole is spread in space across the X–C bonding orbital. Additionally, its configuration and density alter in response to nuclear movement. This non-local and adaptable feature of the hole makes evaluating the Coulomb energy between the hole and the nucleus challenging. For a rough magnitude estimate, we employ a simplified model. To begin with, the portion of the hole population adjacent to the X and carbon nuclei can be considered stationary relative to their respective nuclei and thus has minimal impact on the Coulomb potential variation. It is primarily the segment of the population farther from the nuclei, specifically the middle portion of the X–C bond, that alters the distance to the two nuclei and consequently the Coulomb potential energy. Consequently, we define a partition factor $\alpha$ ($0<\alpha<1$), such that the hole’s effective charge is $\alpha e$, with $e$ representing the elementary charge.

Since the middle portion of the hole cloud is distant from the X and carbon nuclei, we can approximate it as a point charge. In this context, the effective position of the point charge is the average position of the middle portion of the hole cloud. In Sec. 2.2, we showed that there is a node present (absent) in the ground (first-excited) manifold of the XV system as shown in Fig. 5. The presence or absence of the node modifies the shape of the hole cloud, consequently affecting the average position of the effective point charge. In the ground (first-excited) manifold, the average position is deviated from the average position as in an isolated X–C bond toward the carbon (X) nucleus. That is the reason why the electron-phonon interaction strengths differ between the two manifolds. Here, to give an order-of-magnitude estimation, we neglect this difference, examining an isolated X–C bond and assuming that the average position of the middle portion of the hole cloud is at the midpoint of the X–C bond.

The X and carbon nuclei can also be approximated as point charges. Regarding the effective charge of the nuclei, we assume that the middle portion of the hole cloud is distant enough from the two nuclei. Consequently, the inner-shell electrons provide full screening for the charges of the two nuclei, resulting in effective charges of $4e$ for each nucleus.

The Coulomb energy between the hole and a nucleus depends on their relative distance. Being independent of gauge choices, only changes in potential energy are physically meaningful. Therefore, we investigate only the energy variation induced by the displacements of the nuclei.

We first examine the displacement of the carbon nucleus, assuming a stationary X nucleus. For small displacements, as justified in Sec. 2.5.2, the displacements in the $x_{i}$ and/or $y_{i}$ directions approximately result in rotations of the X–C bond, thus no energy variations. In contrast, the displacements in the $z_{i}$ direction cause changes in the relative distance between the hole and the two nuclei. Assuming that the effective position of the hole remains at the midpoint of the X–C bond, the variation in distance between the point charge of the hole and either the X or the carbon nucleus is expressed by

[TABLE]

Additionally, the derivative of the sum of the hole-X and hole-carbon potential energies with respect to $\Delta d$ is given by

[TABLE]

In classical terms, this corresponds to the Coulomb force between the hole and the screened X or carbon nucleus. Here, the minus sign comes from the fact that the hole is of positive charge.

For the time being, we consider the interaction term related to the displacement of the X nucleus, $H_{\rm h-x}$, as an undefined entity.

In conclusion, the expression of the hole-nucleus interaction Hamiltonian as

[TABLE]

As the final step, we convert the above interaction Hamiltonian (26) into the hole basis $\{\ket{\psi_{x}},\ket{\psi_{y}}\}$ as specified in Eq. (7). Recognizing that $\{z_{i}\}$ forms a linear space, choosing an optimal set of basis vectors can simplify the expression of $H_{\rm vo}$. Also, noticing that the coefficients in the spatially polarized eigenstates of the hole, $\ket{\psi_{j}}$ ($j\in\{+,0;\;-,0;\;$ $+,x;\;-,x;\;$ $+,y;\;-,y\}$) in the bases of $\ket{\sigma_{i}}$ are all real (see Sec. 2.2), we may choose the set of basis vectors ($\{\mathbf{V}_{j}\}$) of the displacements of the carbon nuclei following the coefficients of $\ket{\psi_{j}}$, i.e.,

[TABLE]

where $\mathbf{z}_{i}$ is the unit vector in the $z_{i}$ direction and the minus sign is intentionally included.

In this displacement basis, $\{\mathbf{V}_{j}\}$, the interaction Hamiltonian in Eq. (26) becomes

[TABLE]

where variables $\{Q_{j}\}$ and $\{z_{i}\}$ are related via

[TABLE]

and $\mathrm{diag}(\cdot)$ denotes a function that maps a vector to a diagonal matrix with its components placed on the diagonal in order.

Upon transforming the interaction Hamiltonian into the hole orbital basis $\{\ket{\psi_{x}},\ket{\psi_{y}}\}$, we observe that all terms with $j\in\{-,0;-,x;-,y\}$ vanish. These correspond to collective displacement modes with odd inversion symmetry—i.e., the nuclear displacements reverse sign under spatial inversion (see Fig. 8). Likewise, the term $H_{\rm h-x}$ vanishes, as the displacement of the central X nucleus also exhibits odd inversion symmetry. As a result, the interaction Hamiltonian simplifies to

[TABLE]

where the $+$ subscript on $Q_{j}$ is dropped for brevity, and

[TABLE]

The coefficients $F_{0}$ and $F$ are given by

[TABLE]

where $\alpha_{\rm X}$ is the partition factor $\alpha$ of the X nucleus. And, we have

[TABLE]

We have thus derived the electron-phonon interaction Hamiltonian in the hole basis $\left\{\ket{\psi_{x}},\ket{\psi_{y}}\right\}$. Note that the full Hamiltonian of the XV system also includes the free Hamiltonian for the vibrational variables $Q_{x}$ and $Q_{y}$, while the one for $Q_{0}$ may be neglected, as it induces only a global energy shift. The free Hamiltonian takes the form of a two-dimensional harmonic oscillator:

[TABLE]

where $\mu=\mu_{\rm C}$ and $\omega=\omega_{\rm{C},z}$. This Hamiltonian relates to those of the individual carbon displacements $z_{i}$ through

[TABLE]

with

[TABLE]

In summary, we proposed an approximate model characterizing the electron-phonon interaction in the XV systems with undetermined parameters $\alpha_{\rm X}$ ($0<\alpha_{\rm X}<1$). We list the approximations used in App. A. We compare this model with previous theoretical works [11], reported experimental measurements [14], and DFT simulations [13, 15] in Sec. 8.

2.6 Spin-Orbit Interaction

The origin of spin-orbit interaction describes the interaction between the spatial (orbital) movement of the hole and its own spin. It can be explained by the principles of classical physics. According to classical electromagnetism, even when no magnetic field is present in the laboratory frame, an object moving in an electric field experiences a magnetic field in its own rest frame. In the non-relativistic limit ($|\mathbf{v}|\ll c$), this induced magnetic field is given by

[TABLE]

where $\mathbf{v}$ is the object’s velocity, $\mathbf{E}$ is the electric field at its position, and $c$ is the speed of light. If the object possesses spin, the motion of the spin couples to this effective magnetic field generated by its orbital motion, giving rise to the spin-orbit interaction.

In quantum mechanics, the velocity $\mathbf{v}$ is not a well-defined operator and is replaced by the momentum $\mathbf{p}$. The association between them is given by $\mathbf{p}=m_{\rm e}\mathbf{v}$, where $m_{\rm e}$ is the object’s mass. Additionally, the electric field can be expressed as the gradient of the electrostatic potential, $\mathbf{E}=-\nabla V$. Substituting these relations into Eq. (37), the spin-orbit interaction Hamiltonian takes the form [28]

[TABLE]

where $\mathbf{S}=\frac{\hbar}{2}\left[\sigma_{x}^{\rm s},\sigma_{y}^{\rm s},\sigma_{z}^{\rm s}\right]$ is the spin operator with $\{\sigma_{i}^{\rm s}\}$ ($i=x,y,z$) being the Pauli matrices, $g_{\rm s}$ is the particle’s $g$-factor, and $\mu_{\rm B}=\frac{e\hbar}{2m_{\rm e}}$ is the Bohr magneton.

Return to the problem of the spin-orbit coupling of the hole in the XV system, the hole is delocalized over the six X–-C bonds, which surround the central X atom and are approximately confined by the outer six carbon nuclei of the diamond lattice (see Fig. 2). This particular arrangement leads to the consequence that the majority of the spin-orbit interaction is predominantly determined by the electric field generated by the central X nucleus and its inner-shell electrons, while the effect from the electric field generated by the carbon nuclei and their inner-shell electrons can be neglected. To explain, we take the classical picture. For a state with non-zero orbital momentum, the hole can be seen as orbiting around the central X nucleus [29]. Only the electric field generated by the X nucleus and its inner-shell electrons can have non-zero divergence within the area enclosed by the hole’s trajectory, giving rise to non-zero spin-orbit coupling. In contrast, the electric field stemming from any carbon nucleus and its inner-shell electrons has zero divergence within that area. The corresponding magnetic field derived from Eq. (37) reverses with the hole’s motion, resulting in a vanished net spin-orbit coupling effect. Consequently, the spin-orbit interaction is primarily governed by the electric field of the central X nucleus, as in the case of an isolated X atom.

According to the central field approximation [30], the electric field of the central X atom is spherically symmetric and can be expressed as

[TABLE]

where $Z_{\rm X}$ is the atomic number of the X atom, and $\beta_{\rm X}$ ($0<\beta_{\rm X}<1$) is the screening factor (unspecified) accounting for the reduction of the nuclear electric field by inner-shell electrons.

Substituting this form of $\nabla V$ into the spin-orbit Hamiltonian $H_{\rm os}$ yields

[TABLE]

where $\mathbf{L}=\mathbf{r}\times\mathbf{p}$ is the orbital angular momentum of the hole, and $\mu_{0}=\frac{1}{\epsilon_{0}c^{2}}$ is the vacuum permeability. Here we have used $g_{\rm s}\approx 2$ and included the $\frac{1}{2}$ factor arising from Thomas precession [29, 31, 32, 33].

The hole’s orbital angular momentum is defined in three-dimensional space, e.g., $\mathbf{L}=\left[L_{x}^{\rm o},L_{y}^{\rm o},L_{z}^{\rm o}\right]$, where each component $L_{i}^{\rm o}$ (with $i\in\{x,y,z\}$) corresponds to rotation about the respective axis in the XV system [see Fig. 2(a)].

As previously discussed, each X–C bond consists of an $sp^{3}d^{2}$ hybridized orbital from the X atom and an $sp^{3}$ hybridized orbital from the carbon atom. These bonds may be considered as extended $sp^{3}d^{2}$ orbitals oriented along the radial directions. As a result, the orbital angular momentum of the electronic state (see Sec. 2.4) is determined by the orbital angular momentum associated with the $sp^{3}d^{2}$ component.

Since $d$-orbitals have the largest angular momentum quantum number $l=2$, the magnetic quantum number $m_{z}$ of the $sp^{3}d^{2}$ orbital takes on values in $\{-2,-1,0,1,2\}$. Consequently, the $z$-component of the orbital angular momentum can be written as

[TABLE]

where $m_{z}$ is defined with respect to the $z$-axis of the XV system, taken as the quantization axis. The corresponding non-normalized eigenstate of $L_{z}^{\rm o}$ is

[TABLE]

with $\phi^{\rm o}=\arctan\left(\frac{y}{x}\right)$ denoting the azimuthal angle of the hole’s position.

Analyzing the electronic orbital eigenstates as represented by Eqs. (6), the state $\ket{\psi_{+1}}$ ($\ket{\psi_{-1}}$) exhibits a non-zero overlap exclusively with $\ket{m_{z}=+1}$ and $\ket{m_{z}=-2}$ ($\ket{m_{z}=-1}$ and $\ket{m_{z}=+2}$). Additionally, since the longitudinal extent of the X–C bond is small compared to the bond length, the magnitudes of the two overlaps are approximately equal. Thus,

[TABLE]

A visualization of these composition of the $\ket{\psi_{\pm 1}}$ states is shown in Fig. 9.

Utilizing Eqs. (44) and (8), the operator $L_{z}^{\rm o}$ can be represented in the basis $\{\ket{\psi_{x}},\ket{\psi_{y}}\}$ as

[TABLE]

where

[TABLE]

Regarding the $y$-component of the orbital angular momentum, $L_{y}^{\rm o}$, the XV system has two X–C bonds ($\ket{\sigma_{2}}$ and $\ket{\sigma_{6}}$) on the $+y$ side, two ($\ket{\sigma_{3}}$ and $\ket{\sigma_{5}}$) on the $-y$ side, and the remaining two ($\ket{\sigma_{1}}$ and $\ket{\sigma_{4}}$) on the $y=0$ plane. Thus, all states within the XV system exhibit no preference towards either the positive or negative orientation of $L_{y}^{\rm o}$, leading to the conclusion that

[TABLE]

Concerning the $x$-component of the orbital angular momentum, $L_{x}^{\rm o}$, similar to the case of $L_{z}^{\rm o}$, the XV system has three X–C bonds ($\ket{\sigma_{1}}$, $\ket{\sigma_{5}}$, and $\ket{\sigma_{6}}$) on the $+x$ side and the other three ($\ket{\sigma_{4}}$, $\ket{\sigma_{2}}$, and $\ket{\sigma_{3}}$) on the $-x$ side [see Fig. 4]. The six X–C bonds can be regarded as presenting a distorted $\frac{2\pi}{3}$ symmetry about the $x$ axis. Therefore, $L_{x}^{\rm o}$ is smaller than $L_{z}^{\rm o}$ in magnitude but not completely vanish, i.e.,

[TABLE]

The magnitude and form of the operator $L_{x}^{\rm o}$ are left for future examination. Due to its smaller magnitude, this term is often neglected in most of the thesis; nevertheless, it affects the energy degeneracy as discussed in Sec. 3.3.

Therefore, considering only the $z$-component of orbital angular momentum, $L_{z}^{\rm o}$, Eq. (41) is simplified to

[TABLE]

where $\frac{\hbar}{2}\sigma_{z}^{\rm s}$ is the $z$-component of the spin operator of the hole.

Below, we estimate the expectation value, $\left\langle r^{-3}\right\rangle$. The $\ket{\psi_{+1}}$ states can be regarded as comprising half $sp^{3}$ orbitals from the six adjacent carbon atoms, along with half $sp^{3}d^{2}$ hybrid orbitals from the X atom. We neglect the contribution from the $sp^{3}$ half, i.e.,

[TABLE]

The reason is twofold: (1) the spin-orbit coupling strength is reciprocal cubic decay in distance from the X nucleus; (2) the population in the $sp^{3}$ orbitals suffers a more substantial screening effect on the electric field of the X nucleus than that in the $sp^{3}d^{2}$ orbitals.

For the $sp^{3}d^{2}$ half, we estimate the expectation value $\left\langle r^{-3}\right\rangle$ based on a single atom as given by [29]

[TABLE]

where $n_{\rm X}=3,4,5,6$ for Si, Ge, Sn, and Pb, respectively, is the principle quantum number, $l$ is the angular momentum quantum number and $a_{\rm B}=5.29\times 10^{-11}\text{m}$ is the Bohr radius. As shown in Eq. (44), for the $\ket{\psi_{+1}}$ states, the $sp^{3}d^{2}$ half has only two equal-weight components of $\ket{m_{z}=+1}$ and $\ket{m_{z}=-2}$. According to the Aufbau principle, the $\ket{m_{z}=+1}$ component arises from the $l=1$ (i.e., $p$) state, while the $\ket{m_{z}=-2}$ component comes from the $l=2$ (i.e., $d$) state. As such, the expectation value of $r^{-3}$ for the $\ket{\psi_{\pm 1}}$ state can be approximated as

[TABLE]

Hence, the spin-orbit interaction Hamiltonian becomes

[TABLE]

where

[TABLE]

In summary, we proposed an approximate model characterizing the spin-orbit interaction in the XV systems with undetermined parameters $\beta_{\rm X}$ ($0<\beta_{\rm X}<1$). We list the approximations used in App. A. We compare this model with DFT simulations [13] and measured ground splittings [2, 1, 3, 4, 6, 10] in Sec. 8.

3 Intrinsic Properties

In the previous section, we have shown the modeling of the intrinsic interactions of the XV system and derived the intrinsic Hamiltonian. It includes terms of free motions of hole orbitals, atomic vibrations, hole spin, and the interaction between them. The ground and first-excited energy manifolds are analyzed separately with their respective Hamiltonians sharing the same form but with different parameters. Ignoring the term $F_{0}Q_{0}I^{\text{o}}$ in Eq. (30), as it is irrelevant to the splitting of the hole orbital states, the intrinsic Hamiltonian for each manifold is given by

[TABLE]

where $Q_{x}$ and $Q_{y}$ are linear combinations of the movements on the six carbon atomic core in the $\{z_{i}\}$ directions and can be interpreted as the displacements of the center of the mass of the upper three carbon atomic core in the $x$ and $y$ directions, respectively. In this section, we analyze the properties of $H_{0}$.

Before we proceed, to simplify the analysis, we first make a unitary transformation $U_{c}$ for the hole orbital operators $\sigma^{\text{o}}_{i}$ $i=x,y,z$ in Hamiltonian $H_{0}$. The unitary operator $U_{c}$ is defined as

[TABLE]

where

[TABLE]

It transforms the hole-orbital Pauli matrices as

[TABLE]

Consequently, in the new hole orbital bases $\{U_{\rm c}\ket{\psi_{x}},U_{\rm c}\ket{\psi_{y}}\}$, the intrinsic Hamiltonian $H_{0}$ (55) is given by

[TABLE]

where the notations of $\tau_{x}^{\text{o}}$ and $\tau_{y}^{\text{o}}$ have been made uniform with $Q_{x}$ and $Q_{y}$ in terms of subscripts and signs.

3.1 Conserved Quantities

In addition to the total energy, we identify two other conserved quantities in the intrinsic Hamiltonian $H_{0}$. The first is the hole spin, as evidenced by the following commutation relation:

[TABLE]

The second conserved quantity is the $z$-component of the orbital-phonon angular momentum, or vibronic angular momentum, defined as

[TABLE]

[see Eqs. (42), (43), (45), and (58)]. Here, $\phi^{\rm v}=\arctan\frac{Q_{y}}{Q_{x}}$ represents the azimuthal angle in the $Q_{x}$-$Q_{y}$ plane.

The conservation of $L_{z}^{\rm vo}$ is confirmed by the commutation relation

[TABLE]

which can be verified using the following auxiliary relations:

[TABLE]

where $Q=\sqrt{Q_{x}^{2}+Q_{y}^{2}}$, $M_{\pm}=\mathrm{e}^{\pm\mathrm{i}\phi^{\mathrm{v}}}$, and

[TABLE]

These conserved quantities partition the Hilbert space of the XV center system into multiple independent subspaces. As a result, states in different subspaces do not interact. A state initially in one subspace remains confined to that subspace throughout its evolution. Furthermore, the intrinsic Hamiltonian can be block-diagonalized, meaning that each of its eigenstates can be constructed exclusively from basis states within a single subspace.

The conservation of the spin ($\tau_{z}^{\rm s}$) implies that states with spin-up (eigenstates of $\tau_{z}^{\rm s}$ with eigenvalue $+1$) form one subspace, while states with spin-down (eigenvalue $-1$) form another distinct subspace.

Similarly, the conservation of the $z$-component of the vibronic angular momentum ($L_{z}^{\rm vo}$) further divides the spin-up and spin-down subspaces into smaller subspaces, each characterized by a unique value of the $z$-component vibronic angular momentum. A detailed presentation and visualization of these subspaces will be provided in Sec. 3.4.

3.2 Symmetries

The XV center system exhibits both temporal and spatial symmetries: the well-known time-reversal symmetry (associated with motion directions) and a newly defined spatial symmetry termed ”joint-reflection symmetry”.

3.2.1 Time-Reversal Symmetry

The XV center system possesses half-integer spin, and thus the Kramers’ degeneracy theorem [34] applies. This theorem states that the system is at least doubly degenerate due to time-reversal symmetry.

Time-reversal symmetry, or intuitively, motion-reversal symmetry, is based on the time-reversal operation. This operation transforms the momentum and all spins of an object at a given time into their respective opposite directions while keeping the spatial positions unchanged. The time-reversal operator takes different forms in different representations. In the spatial representation (not the momentum representation), the time-reversal operator can be expressed as

[TABLE]

where $K_{0}$ represents the complex conjugation operation, which transforms all imaginary numbers in the quantum states or operators it acts on into their complex conjugates. This operation reverses the momentum operators and the Pauli-$y$ spin operators, as they are imaginary in the spatial representation. The $\sigma_{y}^{\rm s}$ operator reverses the Pauli-$x$ and Pauli-$z$ spin operators. The factor $\mathrm{i}$ is included to make the operator $\mathrm{i}\sigma_{y}^{\rm s}$ real, ensuring it commutes with $K_{0}$.

With the mathematical form of the time-reversal operator, the time-reversal symmetry of Hamiltonian (55) can also be verified by the commutation relation

[TABLE]

or equivalently $\mathcal{T}H_{0}\mathcal{T}^{-1}=H_{0}$. This can be easily proved, as the time-reversal operator transforms only the two operators $\sigma_{y}^{\rm o}$ and $\sigma_{z}^{\rm s}$ in $H_{0}$ by simply flipping their signs, i.e.,

[TABLE]

3.2.2 Joint-Reflection Symmetry

The joint-reflection symmetry is based on the operation of spatial reflections on both the phonon and orbital states, which flips both their positions and momenta, combined with a $\pi$-rotation that flips the $z$-component of the hole spin.

The $Q_{x}-Q_{y}$ plane and the $\tau_{x}^{\text{o}}-\tau_{y}^{\text{o}}$ ”plane” are reflected about the axes $Q_{y}=Q_{x}\tan\theta$ and $\cos\theta\tau^{\text{o}}_{x}+\sin\theta\tau^{\text{o}}_{y}$, respectively, with the same angle $\theta$. The transformations can be presented as

[TABLE]

with a reflection matrix

[TABLE]

The transformation is not unique with valid angle $\theta\in[0,2\pi)$.

Or alternatively, these transformations can be presented as operators $F^{\mathrm{v}}_{\theta}$ and $F^{\mathrm{o}}_{\theta}$, where and operator $F^{\mathrm{v}}_{\theta}$ applies the following transformation to a quantum state $\ket{\psi(Q_{x},Q_{y})}$:

[TABLE]

and

[TABLE]

The joint-reflection operation also involves a $\pi$-rotation on the hole spin, given by the operator

[TABLE]

whose definition is also not unique, with valid $\phi\in[0,2\pi)$.

The joint-reflection operation is then represented by the operator

[TABLE]

One can verify the commutation relation of the intrinsic Hamiltonian with the joint-reflection operator:

[TABLE]

This follows from the symmetry properties of the individual terms in $H_{0}$. The phonon bare Hamiltonian $H_{\rm v}$ possesses reflection symmetry, ensuring its invariance under $F^{\mathrm{vos}}_{\theta\phi}$. In the electron-phonon interaction Hamiltonian $H_{\rm vo}$, the relevant operators appear in pairs and transform in the same manner under the symmetry operation. Finally, in the spin-orbit interaction Hamiltonian $H_{\rm os}$, both contributing operators undergo sign changes, preserving the overall commutation relation. Consequently, the intrinsic Hamiltonian exhibits joint-reflection symmetry.

3.3 Degeneracy

Above, we see that the intrinsic Hamiltonian $H_{0}$ commutes with four operators:

• •

the $z$-component vibronic angular momentum operator $L_{z}^{\rm vo}$,

• •

the spin-$z$ operator $\tau_{z}^{\rm s}$,

• •

the time-reversal operator $\mathcal{T}$, and

• •

the flip-all operator $F_{\theta\phi}^{\rm vos}$.

However, operators among the four do not always commute with each other. Specifically,

[TABLE]

These non-commuting relationships imply two-fold degeneracy in the Hamiltonians $H_{0}$. To show this, we use the following lemma.

Lemma. Let $H$, $A$, and $B$ be operators acting on a Hilbert space $\mathcal{H}$, satisfying:

$[H,A]=[H,B]=0$, 2. 2.

The commutator $[A,B]$ acts as a full-rank operator on $\mathcal{H}$.

Then every eigenvalue of $H$ has at least two linearly independent eigenfunctions.

Proof. Let $\lambda$ be an eigenvalue of $H$, and denote its corresponding eigenspace by

[TABLE]

Since $[H,A]=[H,B]=0$, both $A$ and $B$ preserve $\mathcal{H}_{\lambda}$, i.e., $A\psi\in\mathcal{H}_{\lambda}$ and $B\psi\in\mathcal{H}_{\lambda}$ for all $\psi\in\mathcal{H}_{\lambda}$.

Suppose, for contradiction, that $\mathcal{H}_{\lambda}$ is one-dimensional. Then any vector $\psi\in\mathcal{H}_{\lambda}$ must satisfy $A\psi=a\psi$ and $B\psi=b\psi$ for some scalars $a$ and $b$. It follows that

[TABLE]

which contradicts the assumption that $[A,B]$ acts as a full-rank operator on $\mathcal{H}_{\lambda}$. Therefore, $\dim\mathcal{H}_{\lambda}\geq 2$.

Once an eigenstate of $H$ is given, the proof above also provides a method to find its degenerate counterpart.

Based on this lemma, any one of the relationships (77), (78), and (79) can demonstrate the double degeneracy of the intrinsic Hamiltonian $H_{0}$. It also reveals that, for an eigenstate with a $z$-component vibronic angular momentum expectation value $l$ and spin expectation value $s$, there must be a degenerate state with $-l$ and $-s$. The interaction diagram in Fig. 10 visualizes this symmetry and degeneracy.

In the latter sections, we show that the XV system may also couple to lattice strain (Sec. 4.2) and magnetic fields (Sec. 6). Below, we examine the degeneracy situations when the XV system is exposed to strain and/or the $x$- and/or $y$-components of magnetic fields.

The strain interaction Hamiltonian takes the form

[TABLE]

with a parameter $\theta_{0}$ depending on the stress direction, while the interaction Hamiltonian with the $x$- and/or $y$-components of the magnetic field takes the form

[TABLE]

with a parameter $\phi_{0}$ depending on the magnetic field direction.

The XV center system remains doubly degenerate in energy when strain is applied because the Hamiltonian $H_{0}+H^{\rm str}$ commutes with both $\sigma^{\text{s}}_{z}$ and $\mathcal{T}$, but the two of them do not commute.

Similarly, the XV center system remains doubly degenerate in energy when the the magnetic field in the $x$- and/or $y$-direction are applied because the Hamiltonian $H_{0}+H_{xy}^{\rm mag}$ commutes with both $L_{z}^{\rm vo}$ and $F_{\theta\phi}^{\rm vos}(\phi=\phi_{0})$, but the two of them do not commute. Please note that this result is based on the assumption of a vanished $x$-component of the orbital angular momentum. However, as discussed in Sec. 2.6, it does not strictly vanish [Eq. (48)], which may lead to observable splittings.

However, when both the strain and the magnetic field in the $x$- and/or $y$-direction are present, this double degeneracy will be lifted. This observation suggests that when a constant $x$- and/or $y$-direction magnetic field is applied, the XV system’s splitting depends on the strain’s magnitude. This phenomenon leads to a potential application as follows.

Currently, the in situ strain of the XV center is estimated by numerical calculation of the finite element analysis of the bulk sample, which is computationally costly and not always precise. Since the splitting under the magnetic field in the $x$- and/or $y$-direction depends on strain, it may serve as a means to probe local strain with improved accuracy. We provide the relationship between the splitting, magnetic field, and strain based on the improved quench model (see Sec. 7.2) in Sec. 9.

3.4 Interaction Diagram

The interaction diagram provides a visual representation of the interactions between quantum eigenstates in a certain representation. To construct this diagram, we first select a set of complete basis states for the atomic motions, hole orbitals, and hole spin.

For atomic motions, we use the basis set $\{\ket{n}\ket{m}\}$, obtained by solving the 2D harmonic oscillator in polar coordinates [35]. Specifically, $\ket{n}$ represents the eigenstate of $H_{\rm vib}$ [Eq. (34)]:

[TABLE]

and $\ket{m}=e^{\mathrm{i}m\phi^{\rm v}}$ is the eigenstate of $L_{z}^{\rm v}$:

[TABLE]

Here, $n\in\{0,1,2,\ldots\}$ denotes the total phonon number, and $m\in\{-n,-n+1,\ldots,n-1,n\}$ is the quantum number corresponding to the $z$-component of angular momentum for the atomic motion.

For the hole orbitals, we select the eigenstates of the operator $\tau_{z}^{\rm o}$, denoted as $\ket{+}_{z}^{\rm o}$ and $\ket{-}_{z}^{\rm o}$, as the basis set, i.e.,

[TABLE]

Similarly, for the hole spin, we use the eigenstates of the operator $\sigma_{z}^{\rm s}$, labeled as $\ket{\uparrow}_{z}^{\rm s}$ and $\ket{\downarrow}_{z}^{\rm s}$, i.e.,

[TABLE]

In this basis set, the states $\ket{n}$ and $\ket{n+1}$, corresponding to the same hole orbital and spin, are split by an energy $\hbar\omega$. Similarly, the states $\ket{+}_{z}^{\rm o}$ and $\ket{-}_{z}^{\rm o}$, associated with the same atomic motion and hole spin, are split by an energy $\lambda$. Furthermore, as analyzed in Sec. 3.1, the interactions in $H_{0}$ occur only between states $\ket{m}\ket{+}_{z}^{\rm o}$ and $\ket{m+1}\ket{-}_{z}^{\rm o}$ with the same hole spin because of the conservation of the $z$-component of the total angular momentum.

The interaction strengths can be obtained by expressing the operators

[TABLE]

[see Eq. (63)] in the $\{\ket{n}\ket{m}\}$ basis [36]:

[TABLE]

[TABLE]

Based on this, the minimum interaction strength, which occurs between the two states,

[TABLE]

is given by

[TABLE]

Using these results, we derive the interaction diagram of $H_{0}$, as shown in Fig. 10. From the interaction diagram, one sees that the whole Hilbert space of the XV system can be decomposed into independent subspaces, each characterized by a distinct $z$-component vibronic angular momentum number $l$. The Hamiltonian $H_{0}$ has non-zero matrix elements only between states belonging to the same subspace.

In the latter sections, we show that the XV system may also couple to lattice strain (Sec. 4.2) and magnetic fields (Sec. 6). The strain interaction Hamiltonian takes the form

[TABLE]

with a parameter $\theta_{0}$ depending on the stress direction. The Hamiltonian commutes with the operator $F^{\mathrm{vos}}_{\theta\phi}(\theta=\theta_{0})$ and $\tau_{z}^{\rm s}$, but not $L_{z}^{\rm vo}$.

Additionally, the interaction Hamiltonian with a magnetic field in the $x$- and/or $y$-directions takes the form

[TABLE]

with a parameter $\phi_{0}$ depending on the magnetic field direction. The Hamiltonian commutes with the operator $F^{\mathrm{vos}}_{\theta\phi}(\phi=\phi_{0})$ and $L_{z}^{\rm vo}$, but not $\tau_{z}^{\rm s}$.

The strain interaction Hamiltonian $H^{\rm str}$ has non-zero matrix elements only between the hole orbital states $\ket{+}_{z}^{o}$ and $\ket{+}_{z}^{o}$ with the same phonon number $n$, phononic angular momentum number $m$, and hole spin. One may easily draw the interaction diagram of Hamiltonian $H_{0}+H^{\rm str}$ based on Fig. 10.

On the other hand, the interaction Hamiltonian $H_{xy}^{\rm mag}$ has non-zero matrix elements only between the spin states $\ket{\uparrow}_{z}^{s}$ and $\ket{\downarrow}_{z}^{s}$ of the same phonon number $n$, phononic angular momentum number $m$, and hole orbital state. The cross-spin interaction in $H_{xy}^{\rm mag}$ makes the construction of the interaction diagram for the Hamiltonian $H_{0}+H_{xy}^{\rm mag}$ based on Fig. 10 difficult. A simple way to present it is by a rearrangement of the independent subspaces of $H_{0}$, grouping the states with the same hole orbital movement and different spins close to each other, as shown in Fig. 11.

From the interaction diagrams Fig. 10 and 11, each with two panels, we observed symmetric interaction configurations between the two panels, i.e., for every bare state (eigenstate of Hamiltonian $H_{\rm vib}+H_{\rm os}$, represented by the red and blue solid lines) in the upper panel, there exists a corresponding bare state in the lower panel with the same energy and the same corresponding interaction strengths. As a result, all eigenstates of the Hamiltonians $H_{0}$, $H_{0}+H^{\rm str}$, and $H_{0}+H_{xy}^{\rm mag}$ are two-fold degenerate. In particular, we also observe that the two degenerate states, $\ket{\psi_{l/s}}$ and $\ket{\psi^{\prime}_{l/s}}$, of certain $z$-direction vibronic angular momentum numbers ($l$) or certain $z$-direction hole spin ($s$) are related to each other as

[TABLE]

3.5 Numerical diagonalization

The Hamiltonian $H_{0}$ can be diagonalized numerically. However, this approach is limited to finite photon numbers $n$, requiring a cut-off such that $n\in[0,n_{\rm cut}]$. The cut-off value $n_{\rm cut}$ must be sufficiently large to ensure the accuracy of the derived eigenstates and eigenenergies. To assess this accuracy, we examine their precision as functions of $n_{\rm cut}$. The precision of an eigenenergy $E$ and an eigenstate $\ket{\psi}$ is defined as:

[TABLE]

Here, the reference value $n_{\rm cut}=40$ can be replaced with any sufficiently large number to ensure convergence.

For the lowest 10 eigenenergies (each corresponding to two degenerate eigenstates), when $\{J,\lambda\}<2\hbar\omega$, we find that the precision satisfies $\{\epsilon_{E},\epsilon_{\psi}\}<10^{-4}\;(10^{-15})$ for $n_{\rm cut}>10\;(20)$.

Figure 12 shows the spectrum for the parameters $\hbar\omega=75.6~\mathrm{meV}$ and $\lambda=96.8~\mathrm{meV}$ as $J$ is varied from $0.03\hbar\omega$ to $3\hbar\omega$. Figure 13 illustrates the two spin-up eigenstates with the lowest energies for $J=79.3~\mathrm{meV}$, while using the same parameter set.

4 Stress Effect

4.1 Stress

The stress experienced by an infinitesimal rectangular prism within a solid material is represented by the stress tensor, a second-order tensor expressed as

[TABLE]

where each tensor component, such as $\sigma_{xx}$, is defined as

[TABLE]

Here, $F_{x}$ is the force in the $x$-direction acting on the $y$-$z$ surface of the prism, and $\Delta_{y}$ and $\Delta_{z}$ are the dimensions of this surface. The other components of the tensor are defined analogously. For static stress, the stress tensor is symmetric, satisfying $\sigma_{ij}=\sigma_{ji}$.

Uniaxial stress is a type of stress where the force acts along a single direction, resulting in either stretching or compression along that axis. For uniaxial stress with a magnitude $S$ applied along the direction $\vec{e}=\alpha_{x}\vec{e}_{x}+\alpha_{y}\vec{e}_{y}+\alpha_{z}\vec{e}_{z}$, where $\vec{e}_{i}$ are unit vectors in the $i$-direction and $\alpha_{x}^{2}+\alpha_{y}^{2}+\alpha_{z}^{2}=1$, the stress tensor elements are given by [37]

[TABLE]

This expression follows from the definition of stress as internal pressure, where the stress $\sigma$ is the force divided by the cross-sectional area: $\sigma=\frac{F}{A}$. Both the force and the cross-sectional area must be properly rescaled, as illustrated in Fig. 14, which provides an example for a 2D material.

4.2 Strain

The strain in an infinitesimal rectangular prism within a solid material is described by the strain tensor, a second-order tensor expressed as

[TABLE]

Each tensor component is defined as

[TABLE]

where $\delta_{i}$ is the $i$-th component of the displacement of a point accumulated along the $j$-th direction relative to an arbitrary reference point, and $\Delta_{j}$ is the $j$-th component of the distance from the examined point to the reference point. Fig. 15 illustrates the deformation of a 2D square subjected to strains $\epsilon_{xx}$ and $\epsilon_{xy}$, respectively. Notably, for uniaxial strain, the strain tensor is symmetric, i.e., $\epsilon_{ij}=\epsilon_{ji}$.

Based on the definition of the strain tensor elements, consider a material point initially located at $\left(\Delta_{x},\Delta_{y},\Delta_{z}\right)^{\text{T}}$ that is displaced to $\left(\Delta_{x}^{\prime},\Delta_{y}^{\prime},\Delta_{z}^{\prime}\right)^{\text{T}}$ due to the application of a strain described by the tensor $\left(\epsilon_{ij}\right)$. The relationship between the initial and strained positions is given by

[TABLE]

4.3 Generalized Hooke’s Law

The generalized Hooke’s law describes the relationship between the stress tensor and the strain tensor:

[TABLE]

where $\left(E_{ijkl}\right)$ is the fourth-order elasticity tensor, and the expression involves a tensor product. While $\left(E_{ijkl}\right)$ contains 81 components, both $\left(\epsilon_{ij}\right)$ and $\left(\sigma_{kl}\right)$ are symmetric, i.e., $\epsilon_{ij}=\epsilon_{ji}$ and $\sigma_{kl}=\sigma_{lk}$, reducing each tensor to six independent terms. Consequently, $\left(E_{ijkl}\right)$ can be represented by a $6\times 6$ elastic matrix $C$.

For cubic lattices, such as diamond, the generalized Hooke’s law simplifies to:

[TABLE]

where the matrix $C$ is characterized by three independent parameters: $C_{11}$, $C_{12}$, and $C_{44}$. For diamond, these parameters are $C_{11}=1075\,\mathrm{GPa}$, $C_{12}=139\,\mathrm{GPa}$, and $C_{44}=567\,\mathrm{GPa}$.

4.4 Strain Engineering

Natural strains can arise from various lattice defects near the XV center within the diamond bulk. Artificial strain, on the other hand, can be induced by applying engineered external forces to the diamond. One method to implement artificial strain in an integrated optical system involves using electric forces by applying a voltage to electrodes positioned on a cantilever etched from a substrate, as shown in Fig. 16(a).

Given the size of the cantilever and electrodes shown in Fig. 16(a), along with the voltage $V$ applied between the electrodes on the cantilever and the substrate, we derive a theoretical formula to estimate the magnitude of the induced strain.

First, we perform a force analysis for the cantilever. As shown in Fig. 16(b), the cantilever experiences two forces: the stress from the substrate and the electric force $F$. We assume that the stress has a maximum value of $\mathrm{c}_{11}\epsilon_{xx}^{\rm R}$ at the upper surface of the cantilever, decreases linearly to zero at the middle plane, reverses direction, and increases linearly to $-\mathrm{c}_{11}\epsilon_{xx}^{\rm R}$ at the bottom surface. For static equilibrium, the torques caused by these two forces must balance, yielding:

[TABLE]

where $z^{\prime}$ is the coordinate along the $z$-direction with the origin at the middle plane of the cantilever, and $F$ is the electric force, assumed to act at the center of the upper electrode.

The electric force $F$ is given by the product of the electric field strength $E_{1}$ (assumed uniform near the upper electrode) and the charge $Q$ on the upper electrode:

[TABLE]

Due to the high relative permittivity of diamond ($5.5$–$10$), the electric field within the diamond bulk is nearly zero. For simplicity, we consider an equivalent case where the upper electrode is placed at the bottom of the cantilever. Next, we approximate the electric field induced by the upper electrode. As shown in Fig. 16(c), we assume the electric field inside the semi-cylinder covering the electrode is uniform, given by:

[TABLE]

where $\varepsilon_{0}=8.85\times 10^{-12}\,\mathrm{F/m}$ is the vacuum permittivity.

For the case $h>\frac{y}{2}$, we approximate the electric field $E_{2}$ outside the semi-cylinder as that induced by a charged cylinder:

[TABLE]

Note that $E_{1}$ and $E_{2}$ are the electric fields induced by the upper electrode alone, without accounting for the presence of the electrode placed on the lower substrate [see Fig. 16(a)].

Given that the electrode on the lower substrate is a large plane, the method of images can be used to calculate the voltage between the two electrodes. The main idea is that the electric field induced by the upper electrode and the substrate electrode is equivalent to the field induced by the upper electrode and an image electrode, located on the opposite side of the substrate plane and carrying an opposite charge, as shown in Fig. 16(d). Using this approach, the voltage between the upper and substrate electrodes can be calculated as

[TABLE]

From Eqs. (105)–(109), the strain-voltage relationship is derived as

[TABLE]

The above derivation assumes that the cantilever remains stationary. In reality, the cantilever bends toward the substrate electrode. For this case, $h$ can be replaced with $h^{\prime}=h-\frac{1}{2}\epsilon_{xx}(l+x)$ in Eq. (110) to obtain a more accurate relationship. However, this substitution introduces additional complexities in solving the equation.

Meesala et al. [14] demonstrated artificial strain experimentally using a triangular-prism cantilever, as shown in Fig. 17.

From Fig. 17(a), the triangular prism is estimated to have the dimensions listed in Tab. 2.

Using the finite element method (FEM), Meesala et al. [14] calculated the strain at the position of the SiV center near the joint of the cantilever, as shown in Fig. 17(b), under different electric voltages. Tab. 3 presents their calculated results for the strain $|\epsilon_{xx}^{\rm M}-\epsilon_{yy}^{\rm M}|$, where

[TABLE]

in the SiV coordinate system.

Separately, the setup of Meesala et al. can be approximated using the rectangular-prism cantilever model discussed above. Using this model and Eq. (110), Table 3 also lists the calculated strain $|\epsilon_{xx}^{\rm R}-\epsilon_{yy}^{\rm R}|$, where $\epsilon_{yy}^{\rm R}\approx-0.11\epsilon_{xx}^{\rm R}$. Since the $x$-$y$ plane of the cantilever coordinate system aligns with that of the SiV coordinate system, the values of $|\epsilon_{xx}^{\rm R}-\epsilon_{yy}^{\rm R}|$ and $|\epsilon_{xx}^{\rm M}-\epsilon_{yy}^{\rm M}|$ are expected to match. As shown in Table 3, the close agreement between the FEM results and the rectangular-prism cantilever model demonstrates the effectiveness of Eq. (110) as a practical tool for guiding cantilever design.

Finally, it is worth noting that mechanical forces have also been employed to induce strain in diamonds. Examples include compressing a diamond sample with two anvils [38], bending a diamond needle with an indenter [39], and stretching a diamond bridge with a gripper [40], as illustrated in Fig. 18(a), (b), and (c), respectively.

4.5 Strain Interaction Hamiltonian

In Sec. 2.5, we derived the electron-phonon interaction Hamiltonian, which describes the interaction between the hole orbital and the atomic displacements. The interaction between the hole orbitals and strain arises from the same mechanism, as strain can also induce atomic displacements, similar to thermal vibrations. Therefore, the strain interaction Hamiltonian can be derived by expressing the atomic displacement magnitudes $Q_{0}$, $Q_{x}$, and $Q_{y}$ (as in Eq. 30) in terms of the strain tensor elements $\epsilon_{k}$ ($k\in\{zz,xx,yy,xy,yx,xz,zx,yz,zy\}$) using the dimensions of the XV system.

As defined in Sec. 4.2, each of the nine elements in the strain tensor represents a distinct type of lattice distortion. Given the dimensions of the XV system (Fig.4), one can determine the atomic displacements induced by each strain component. These displacements, corresponding to the nine strain tensor elements, are illustrated in Fig. 19.

With the displacements of each atom under strain known for each strain tensor element, we can derive the generalized displacement vectors – concatenations of the displacement vectors of multiple atoms – for each strain tensor element. Since the displacement of the central X atom decouples from the hole orbitals, and the displacements of the lower ($-z$ side) three carbon atoms are symmetric to the upper ($+z$ side) three, we decompose the generalized displacement vector of $\epsilon_{k}$ as

[TABLE]

where $\mathbf{D}_{k}^{\rm upp}$ and $\mathbf{D}_{k}^{\rm low}$ are the generalized displacement vectors of the upper and lower three carbon atoms, respectively. The expressions for the set $\{\mathbf{D}_{k}^{\rm upp}\}$ are given by

[TABLE]

expressed in the $x^{1},y^{1},z^{1},x^{2},y^{2},z^{2},x^{3},y^{3},z^{3}$ basis, where $x^{i}$, $y^{i}$, and $z^{i}$ are the coordinates of the $i$-th carbon atom in the XV Cartesian coordinate system (see Fig. 2).

On the other hand, in the same basis, the unit vectors in the three generalized directions $\mathbf{V}_{j}$ ($j\in\{0+,x+,y+\}$) as in Eq. (27) can also be decomposed as

[TABLE]

with $\{\mathbf{V}_{j}^{\rm upp}\}$ given by:

[TABLE]

where

[TABLE]

are the unit vectors in the $z_{i}$ ($i=1,2,3$) directions, respectively.

Due to the symmetry between the displacements of the upper and lower three carbon atoms, when projecting the generalized displacement vector $\mathbf{V}_{j}$ onto the generalized direction $\mathbf{V}_{j}$, we can use the relationship

[TABLE]

The projections (inner products) are shown in Tab. 4.

With the projection coefficients, and noting that $\epsilon_{xy}=\epsilon_{yx}$, $\epsilon_{xz}=\epsilon_{zx}$, and $\epsilon_{yz}=\epsilon_{zy}$, we can express the atomic displacement magnitudes $Q_{0}$, $Q_{x}$, and $Q_{y}$ in terms of the strain tensor elements $\epsilon_{k}$ as

[TABLE]

Substituting these equations into Eq. (30), we obtain the strain interaction Hamiltonian as

[TABLE]

where

[TABLE]

are the electron-strain coupling strengths or strain susceptibilities.

The ratio between $d$ and $f$ is given by

[TABLE]

4.6 Stress-Strain-Spectrum Converter

With the Hamiltonian (119) describing the strain effect, the energy shifts for the states in both the ground and excited manifolds can be calculated. Furthermore, a transition from a state in the excited manifold to a state in the ground manifold is accompanied by an emitted photon. The average energy difference of the emitted photons, corresponding to the zero-phonon line (ZPL), is given by:

[TABLE]

where $E_{i}$ is the eigenenergy of the $\ket{i}$ state, $\ket{A}$ and $\ket{C}$ are the first and second spin-up states in the excited manifold, and $\ket{1}$ and $\ket{3}$ are the first and second spin-up states in the ground manifold. The ZPL of an unstrained sample can thus be measured to determine this energy difference. Using this information, along with the formulas in Secs. 4.1–4.5, the spectrum of the XV emitter under arbitrary stress can be computed.

To facilitate these calculations, we developed a user-friendly calculator, S-cubed, which performs stress-strain-spectrum conversions, as shown in Fig. 20. This MATLAB-based tool features a convenient user interface (UI) that allows users to input the ZPL wavelength, the splittings of the lowest two spin-up states in the ground and excited manifolds under zero stress, and the magnitude and direction of uniaxial stress. The tool provides the following outputs:

Visualization of the deformation of the XV atomic system. 2. 2.

Strain tensor in both the diamond and XV coordinate systems. 3. 3.

Variation of the splitting of the two lowest spin-up states in the ground and first-excited manifolds with respect to strain magnitude. 4. 4.

Change in the ZPL wavelength as a function of strain magnitude. 5. 5.

Evolution of the spectrum as a function of stress magnitude.

5 Stark Effect and Light Emission

5.1 Electric Dipole Moment

In classical physics, the electric potential energy of a particle with charge $Q_{\rm p}$ in a uniform electric field $\mathbf{E}$ is given by

[TABLE]

where $\mathbf{r}$ is the displacement of the particle from a reference point with zero potential energy. This leads to the definition of the electric dipole moment as

[TABLE]

In quantum mechanics, the physical quantity $\mathbf{r}$ is represented as an operator $\hat{\mathbf{r}}$. In the basis set $\{\ket{\psi_{i}}\}$, the dipole moment operator $\hat{\mathbf{d}}$ is represented by a matrix

[TABLE]

where $\hat{I}=\sum_{i}\ket{\psi_{i}}\bra{\psi_{i}}$ is the identity operator in the basis $\{\ket{\psi_{i}}\}$ and $\left[\bra{\psi_{i}}\hat{\mathbf{r}}\ket{\psi_{j}}\right]$ is a matrix with elements $\bra{\psi_{i}}\hat{\mathbf{r}}\ket{\psi_{j}}$. When the basis states $\{\ket{\psi_{i}}\}$ possess inversion symmetry, the diagonal elements $\{\bra{\psi_{i}}\hat{\mathbf{r}}\ket{\psi_{i}}\}$ are all zero since the position operator is odd under inversion symmetry. Consequently, the external electric field $\mathbf{E}$ only induces transitions between different states of $\{\ket{\psi_{i}}\}$.

5.2 Robustness Against Charge Noise

We can estimate the electric dipole moment based on the eigenstates of the hole orbitals in the XV Coulomb potential, Eq. (7), and the size of the XV system as shown in Fig. 4.

Given the estimated geometry of the X-C bond, the expected position of the hole in an X-C bond state $\ket{\sigma_{i}}$ ($i=1,2,\ldots,6$) is approximately at the center of the X-C bond. Based on the dimensions of the XV system (see Fig. 4), the operators $\hat{x}^{\rm o}$, $\hat{y}^{\rm o}$, and $\hat{z}^{\rm o}$ in the X-C bond basis $\{\ket{\sigma_{i}}\}$ are given by

[TABLE]

Using the expressions in Eq. (7), we calculate the electric dipole moments in the subspace $\{\psi_{x},\psi_{y}\}$. However, all matrix elements vanish:

[TABLE]

The vanishing inner products arise from the parity properties of $\ket{\psi_{x}}$, $\ket{\psi_{y}}$, and the operators $\hat{x}^{\rm o}$, $\hat{y}^{\rm o}$, and $\hat{z}^{\rm o}$ as functions of spatial coordinates. The state $\ket{\psi_{x}}$ ($\ket{\psi_{y}}$) is odd under $\hat{x}^{\rm o}$ ($\hat{y}^{\rm o}$) and even under $\hat{y}^{\rm o}$ ($\hat{x}^{\rm o}$). Both $\ket{\psi_{x}}$ and $\ket{\psi_{y}}$ share the same parity with respect to $\hat{z}^{\rm o}$, being odd (even) in the ground (first-excited) subspace. The operators $\hat{x}^{\rm o}$, $\hat{y}^{\rm o}$, and $\hat{z}^{\rm o}$ are odd functions of their respective coordinates. For example, in the first integral above, $\ket{\psi_{x}}\hat{x}^{\rm o}\ket{\psi_{y}}$ is even in $\hat{x}^{\rm o}$ but odd in both $\hat{y}^{\rm o}$ and $\hat{z}^{\rm o}$, making the integral zero. Similar reasoning applies to the other two integrals. Thus, the electric dipole moment in this subspace is zero.

Experimental observations confirm that the electric dipole moment is at least four orders of magnitude smaller than that of the nitrogen-vacancy (NV) center in diamond [8]. This characteristic has both advantages and disadvantages. On the one hand, the reduced dipole moment narrows the linewidths of energy levels, making them less susceptible to charge noise. On the other hand, it limits the feasibility of using electric fields to tune energy level positions.

5.3 Stark Shift

Although the electric dipole moments in both the ground and the first-excited subspaces are zero, the electric dipole moment in the space combined with these two is not zero. The off-diagonal block matrix of the three coordinates, mapping from the space spanned by $\{\psi_{x-},\psi_{y-}\}$ to the space spanned by $\{\psi_{x+},\psi_{y+}\}$, is given by

[TABLE]

Hence the full matrix in the combined space spanned by $\{\psi_{x-},\psi_{y-},\psi_{x+},\psi_{y+}\}$ is given by

[TABLE]

where 0 is the zero matrix and $I^{\rm 0}$ is the identity matrix.

Now we define

[TABLE]

then the there components of the electric dipole moment can be presented as

[TABLE]

where $\tau_{+}^{\rm o}$ are the raise and lower operators for the electronic orbital states.

Lastly, we must calculate the electric dipole moment in the bases of the eigenstates of the intrinsic Hamiltonian $H_{0}$. We do this numerically by adopting the SnV-G and SnV-E parameters based on the ab initio calculation [13] (see Tab. 6). The block matrix mapping from the space spanned by the lowest eight spin-up states $\{\ket{\rm A},\ket{\rm C},\ket{\rm E},\cdots,\ket{\rm O}\}$ in the first-excited subspace to the space spanned by the lowest eight spin-up states $\{\ket{1},\ket{3},\ket{5},\cdots,\ket{15}\}$ in the ground subspace is given by

[TABLE]

[TABLE]

[TABLE]

Here, the symbol D denotes the unit of the Debye (1 D $=3.336\times 10^{-30}$ C$\cdot$m). The fraction in parentheses beside the basis states on the border of the matrix represents the conserved quantity $\left\langle L_{z}^{\rm vo}\right\rangle/\hbar$ (see Sec. 3.1). It is evident that the matrix element $\hat{d}_{\pm}[i,j]$ ($\hat{d}_{z}[i,j]$) is nonzero only when the transition satisfies $\Delta\left\langle L_{z}^{\rm vo}\right\rangle/\hbar=\pm 1$ (0).

Using the nonzero matrix elements of the electric dipole moment, we estimate the Stark shifts induced by a static electric field $\mathbf{E}$. According to perturbation theory [29], the second-order energy shift (the first-order contribution, corresponding to the diagonal elements, vanishes) for a SnV system is given by

[TABLE]

where $E_{i}$ and $E_{f}$ are the energies of the initial and final states, respectively.

The change in the splitting between $E_{\rm A}$ and $E_{1}$ is expressed as

[TABLE]

For both $\Delta E_{\rm A}$ and $\Delta E_{1}$, we consider the two largest terms in their respective expansions, yielding

[TABLE]

Here, $\theta_{z}$ is the angle between the applied electric field $\mathbf{E}$ and the $z$ axis of the XV system coordinate, and $E=|\mathbf{E}|$. Substituting the values $E_{\rm A}-E_{1}=484.1\,\mathrm{THz}$ (corresponding to 619.3 nm), $E_{\rm A}-E_{3}=483.2\,\mathrm{THz}$, and $E_{\rm B}-E_{1}=486.9\,\mathrm{THz}$ into the equations, we obtain

[TABLE]

where the minimum (maximum) value corresponds to $\theta_{z}=\frac{\pi}{2}$ ([math]). This result agrees with experimental observations [8, 9], as shown in Fig. 21.

5.4 Light Emission

The nonzero elements of the electric dipole moment indicate that the XVs can couple to photons, making them photon-addressable quantum devices. Different spatial components of the electric dipole moment couple to photons with specific polarizations. Specifically, the $\hat{d}_{z}$ component couples to photons propagating in the $x$-$y$ plane with an electric field component along the $z$ axis. In contrast, the $\hat{d}_{+(-)}$ components couple to right (left) circularly polarized photons propagating along the $z$ axis, with their electric fields oscillating in the $x$-$y$ plane.

In experiments, the photon emission and absorption profiles of XV emitters can be characterized using Photoluminescence (PL) and Photoluminescence Excitation (PLE) techniques, respectively.

In a PL experiment, a laser with a fixed wavelength excites the emitter to an energy level above the first-excited manifold. The emitter first relaxes to states within the first-excited manifold, populating them according to a specific distribution, and subsequently transitions from the first-excited manifold to the ground manifold, accompanied by photon emission.

In contrast, a PLE experiment begins by initializing the XV emitters in the ground manifold with a defined population distribution. A wavelength-tunable laser is then used to scan the spectrum, identifying frequencies that can excite the emitter. The excited emitter subsequently decays, emitting a photon at the same frequency as the excitation laser, which is detected by a photon detector.

Thus, the PL spectrum provides information about the emission characteristics of the emitter, while the PLE spectrum reveals its absorption characteristics.

To simulate the PL and PLE spectra, it is necessary to calculate the relaxed population distribution in the first-excited and ground manifolds, respectively. Both distributions follow the Boltzmann distribution, which depends on the temperature:

[TABLE]

where $n_{i}$ represents the population of the $i$-th energy level of energy $E_{i}$, $E_{1}$ is the the lowest-energy within the respective manifold, $k_{B}$ is the Boltzmann constant, and $T$ is the absolute temperature. The normalization factor for the distribution is omitted here, as we are concerned only with the relative strengths of the photon peaks in the PL or PLE spectrums.

Table 5 provides the calculated energy differences $\Delta E_{i}$ and the corresponding populations $n$ for the lowest eight energy levels in the ground and excited manifolds at various absolute temperatures. These calculations are based on the SnV parameters obtained from ab initio calculations [13], as listed in Table 6.

Both the population distribution and the magnitude of the electric dipole matrix elements determine the intensity of the collected light in PL and PLE experiments. According to Fermi’s golden rule [42], the transition intensity from an initial state $\ket{i}$ to a final state $\ket{j}$, associated with the $k$-component ($k=+,-,z$) of the electric dipole moment, $\hat{d}_{k}$, is given by

[TABLE]

where $n_{i}$ is the population of the initial state $\ket{i}$.

Figure 22 presents the calculated magnitudes for transitions in the PL and PLE spectra for a SnV center at 100 K.

From Fig. 22(a), the ratio between the intensities of the $A1$ and $A3$ lines, $\frac{I_{+}^{A3}}{I_{z}^{A1}}$, is approximately 0.16. This can be understood as follows: due to the similarity between the eigenstates in the ground and first-excited manifolds,

[TABLE]

where $q_{\rm SnV-G}$ is the quench factor for the $\tau_{x}^{o}$ and $\tau_{y}^{o}$ operators in the ground manifold. Thus,

[TABLE]

which closely matches the simulated result of 0.16.

Next, we consider the proportion of light emitted into zero-phonon lines (ZPLs). Since the states $\ket{A},\ket{C},\ket{3},$ and $\ket{1}$ all reside in the phononic ground states of their respective manifolds, the transitions $\ket{A},\ket{C}\to\ket{1},\ket{3}$ do not involve phonon excitation. Hence, the four lines A1, A3, C1, and C3 are all ZPLs.

For a temperature of 50 K, the calculated proportions of light emitted in the ZPLs for SiV, GeV, SnV, and PbV are 0.84, 0.79, 0.83, and 0.85, respectively.

6 Zeeman Effect and Addressability Manipulation

6.1 Zeeman Hamiltonian

The interaction between the external magnetic field and the XV system, known as the Zeeman effect, affects both the hole spin and the hole’s orbital angular momentum. In Sec. 2.6, we established that the hole in the XV system has only a $z$-component of angular momentum. Consequently, the Zeeman Hamiltonian is given by [11]:

[TABLE]

where $\mu_{\rm B}=14.1\,\text{GHz/T}$ is the Bohr magneton, and $B_{i}$ ($i=x,y,z$) is the $i$-component of the magnetic field. Here the electron $g$-factor $g_{L}=1$ and we set the spin $g$-factor $g_{s}=2$. The hole orbital angular momentum is defined as $L_{z}^{\rm o}=-\frac{\hbar}{2}\sigma_{y}^{\rm o}$ [see Eqs. (42), (43), and (45)]. Note that the factor $-\frac{1}{2}$ is not apparently shown in Hepp’s thesis [11].

The Zeeman effect induces energy splitting; however, as discussed in Sec. 3.2, the $x$- and $y$-components of the magnetic field do not break the degeneracy of the XV system. Thus, only the $z$-component of the magnetic field causes energy splitting. Fig. 23 illustrates the energy levels within an XV manifold under various magnetic field strengths applied along the $z$- and $x$-directions, respectively.

6.2 Addressability Manipulation

The Zeeman effect is a crucial mechanism for selectively addressing transitions between specific quantum states in the XV system. Without an external magnetic field, the spin-up and spin-down states in the ground manifold ($\ket{\uparrow}_{1}$ and $\ket{\downarrow}_{2}$) are degenerate, as are the spin-up and spin-down states in the first-excited manifold ($\ket{\uparrow}_{\rm A}$ and $\ket{\downarrow}_{\rm B}$). Consequently, coupling light to the $\ket{\uparrow}_{1}\leftrightarrow\ket{\uparrow}_{\rm A}$ transition unavoidably also couples the $\ket{\downarrow}_{2}\leftrightarrow\ket{\downarrow}_{\rm B}$ transition with equal strength, as illustrated in Fig. 24(a).

This degeneracy presents challenges for utilizing the XV system in quantum information processing. For example, several protocols require asymmetric driving of the transitions between $\ket{\uparrow}_{1}\leftrightarrow\ket{\uparrow}_{\rm A}$ and $\ket{\downarrow}_{2}\leftrightarrow\ket{\downarrow}_{\rm B}$, such as:

• •

cnot gates between two XV qubits driven optically [43].

• •

cnot gates between a photon and an XV qubit [44].

• •

Single-photon generation [45].

• •

Single-qubit gates with optical driving [46, 47, 48].

To break the symmetry between the two transition energies, both a $z$-component and a combination of $x$- and $y$-components of the external magnetic field are required. The $x$- and/or $y$-components create spin orientations that differ between the ground and first-excited manifolds. When the $z$-component magnetic field is applied, the spin splitting $\Delta_{\rm AB}$ in the excited manifold differs from the spin splitting $\Delta_{12}$ in the ground manifold, as shown in Fig. 24(b) hence break the symmetry in the two transitions.

The difference in spin orientation between the first-excited manifold and the ground manifold arises from the differing spin-orbit coupling strengths in the two manifolds. As explained in Sec. 2.6, the hole experiences an intrinsic magnetic field along the $z$-axis due to spin-orbit coupling [see Eq. (37)]. This intrinsic magnetic field differs between the two manifolds because of their distinct spin-orbit coupling strengths. Consequently, when identical $x$- and/or $y$-components of the external magnetic field are applied to both manifolds, the combined magnetic fields in the two manifolds adopt different orientations, resulting in distinct spin orientations.

Fig. 25(a) shows the spin orientation angle

[TABLE]

where $|{\rm G_{\downarrow}}\rangle$ represents the lowest spin-down state in a single manifold, as a function of the $x$-component of the magnetic field strength ($B_{x}$). Fig. 25(b) illustrates the orientation difference $\theta^{\rm s}_{\rm e}-\theta^{\rm s}_{\rm g}$ between the first-excited and ground manifolds as a function of $B_{x}$. The calculations reveal that the SnV system can exhibit a maximum spin orientation difference of $0.18\pi$ at $\mu B_{x}=4\,\mathrm{meV}$.

7 Approximation Methods

The phonon, whose Hilbert space is infinitely large, introduces coupling that makes the XV center system difficult to solve exactly. Two popular approximation methods are used to simplify the treatment of electron-phonon coupling.

The Born-Oppenheimer (B-O) approximation allows one to treat the phonon motion and hole motion separately, but the approximation is valid only when the phonon energy is far lower than the electron-phonon coupling strength. This prerequisite is not fulfilled in the XV systems, hence the B-O approximation is not applicable. However, this method offers a valuable perspective for XV systems with well-developed concepts, such as the Jahn-Teller energy.

The quench-factor method treats the electron-phonon eigenstates as black boxes and approximates the electron-phonon coupling effect as merely quenching the strength of the other interactions involving the hole. The original method assumes that the strengths of all other hole-related interactions are far smaller compared to the electron-phonon interaction. However, this condition often fails, particularly in systems with strong spin-orbit coupling. In this section, we introduce an improved definition of the quench factors, enabling a more accurate description of XV systems.

7.1 Born-Oppenheimer Approximation

The Born-Oppenheimer (B-O) approximation is inspired by the assumption that the hole, which has a light effective mass and strong bonding to the atomic core, can adapt almost instantaneously to the atomic motion. Therefore, the atomic displacements are approximately constant from the perspective of the hole, i.e., the state of the hole depends only on the atomic displacements and is independent of the atomic momenta.

A rigorous derivation of the B-O approximation is provided by the adiabatic theorem [49]. It states that if a quantum system with a time-dependent Hamiltonian $H(t)$ is initially in an eigenstate $\ket{\psi_{i}(t=0)}$ of $H(0)$, it will remain in the corresponding eigenstate of $H(t)$, $\ket{\psi_{i}(t)}$, provided that the following condition holds [50]:

[TABLE]

for all other eigenstates $\left\{\ket{\psi_{j}}\right\}$ ($j\neq i$) with corresponding eigenenergies $\{E_{j}\}$. Here, $h=6.626\times 10^{-34}\ \text{J}\cdot\text{s}$ is the Planck constant.

Based on the adiabatic theorem, the hole orbital will remain in its eigenstate if the atomic displacement changes slowly enough and the gap between the hole’s eigenenergies is sufficiently large. When these conditions are fulfilled, the eigenenergy of the hole orbitals, obtained by diagonalizing the Hamiltonian $H_{\rm vo}$ (see Eq. 30) while treating $Q_{x}$ and $Q_{y}$ as constants, is given by

[TABLE]

The eigenenergies $E_{\pm}^{\rm o}$ form surfaces in the $Q_{x}$-$Q_{y}$ space, taking the shape of cones as shown in Fig. 26(a). We observe that the two eigenenergies touch at the point $(Q_{x},Q_{y})=(0,0)$. The condition (146) can only be met when the radial atomic displacement $\sqrt{Q_{x}^{2}+Q_{y}^{2}}$ is sufficiently large.

Assuming that this condition is satisfied, the harmonic potential in which the atoms move is modified by the addition of the term $E_{\pm}^{\rm o}$, resulting in the potential

[TABLE]

This potential is referred to as the adiabatic potential energy surface, as shown in Fig. 26(b). It has an upper sheet and a lower sheet. The lower sheet $E_{-}^{\rm v}$ has minima forming a ring at a distance

[TABLE]

from the origin. Additionally, there is a spike at the origin with a height

[TABLE]

This height $E_{\rm JT}$ is called the Jahn-Teller energy.

We can infer that if the phonon energy is much smaller than the spike height:

[TABLE]

the atomic displacement vector $(Q_{x},Q_{y})$ is unlikely to appear near the origin, ensuring that the adiabatic condition (146) is fulfilled. Consequently, Eq. (151) is the alternative prerequisite condition for the B-O approximation.

Lastly, Note that

[TABLE]

[see Eq. (91)]. The condition (151) is equivalent to

[TABLE]

which indicates that the electron-phonon coupling strength is much larger than the phonon energy. This is consistent with the assumption that the atoms move slowly compared to the hole orbitals.

The condition (153) implies strong electron-phonon coupling, meaning that the electron can almost instantaneously adapt to the motion of the phonons. However, this assumption does not always hold. The reason is straightforward: although a single hole generally exhibits high mobility and adaptability, in the case of the XV system, the hole is also bonded to the central X atom. This bonding increases its effective inertia, making it less responsive to the motion of the surrounding carbon atoms. As a result, the Born-Oppenheimer approximation no longer applies. Nevertheless, the concept of the Jahn-Teller energy (150) remains useful as a qualitative reference.

7.2 Quench Factors

Quench factors, or Ham factors, originally proposed by Ham [51], are used to simplify the characterization of the influence of the phonon on the interactions related to the electronic orbital states. The name of quench factors refers to two factors both less than one. It reflects the nature of the phonon effect that it effectively reduces the interaction strengths related to the electronic orbital states.

This approximation technique reduces the complexity of representing the effect of phonons, which situate in an infinite-dimensional Hilbert space, into two scalar parameters. As a result of this convenience, this method was initially utilized in the SiV system [11] and subsequently in other XV systems [13, 14]. Although widely used, few have examined its validity and precision. As shown in Sec. 2.6 and DFT simulations [13], the spin-orbit coupling strength increases along with the atomic number of the X atom. Therefore, there could be potential problems for the XV systems of larger X atomic number to inherit this approximation from the SiV system.

This section aims to establish a foundation for applying the approximate quench factor model. Firstly, we show that the quench factors based on original definitions are less precise for XV systems with strong spin-orbit coupling. Moreover, the accuracy of this approach depends on the relative electron-phonon coupling strength compared to the phonon energy. To address these issues, we introduce refined definitions of the quench factors, which are well-suited to XV systems containing elements with large atomic numbers. While these improved quench factors offer enhanced accuracy over the original formulation, their precision still depends on the ratio of the electron-phonon coupling strength to the phonon energy.

7.2.1 Original Definition

In the initial proposal of the quench factor [51], Ham considers the following Hamiltonian

[TABLE]

where $H_{\rm vib}$ and $H_{\rm vo}$ are as the ones defined in this thesis. And, $H_{\rm int}$ is the interaction Hamiltonian excluding the electron-phonon interaction. It may incorporate the electronic orbital operators $\tau_{z}^{\rm o}$, $\tau_{x}^{\rm o}$, and $\tau_{y}^{\rm o}$.

Without considering the electron spin, the Hamiltonian $H_{\rm vib}+H_{\rm vo}$ is two-fold degenerate in energy. If the interaction strength of $H_{\rm int}$ is far less than the ground splitting of $H_{\rm vib}+H_{\rm vo}$, the two degenerate ground states $\{\ket{1},\ket{2}\}$ effectively form a nearly closed subspace of the Hamiltonian $H$. Therefore, if limiting our scope to the two energetically lowest eigenstates of $H$, it is sufficient to examine the projected Hamiltonian $PHP$, where

[TABLE]

is the projection operator onto the subspace $\{\ket{1},\ket{2}\}$.

Since the states $\ket{1}$ and $\ket{2}$ are degenerate, the term $P(H_{\rm vib}+H_{\rm vo})P$ vanishes under a suitable energy gauge. The effective Hamiltonian becomes

[TABLE]

When evaluating the operators $\left\{\;P\,\tau_{i}^{\rm o}\,P\;\right\}$ ($i=x,y,z$), we find the following relationships:

[TABLE]

They can be easily proved using the symmetries as discussed in Sec. 3.2. As a result, the effective Hamiltonian can be expressed as

[TABLE]

where

[TABLE]

These are Pauli operators in the bases $\{\ket{1},\ket{2}\}$, which are entangled states involving electron orbitals and phonons. As an entity, the combination of them is also referred to as a vibron.

The effective Hamiltonian $H_{\rm Ham}$ (excluding electron spin) is two-dimensional, hence, this approximation method significantly simplifies the calculation of the ground splitting of $H$. As a concrete example, when $H_{\rm int}=H_{\rm os}$, which is the spin-orbit coupling Hamiltonian, the ground splitting is approximately given by

[TABLE]

Due to the inner product nature, $0\leq p,q\leq 1$, which gives rise to their name as quench factors. This approximation tells us that the effect of electron-phonon coupling is to quench the coupling strengths of other interactions associated with the electronic orbital states. For the concrete values of $p$ and $q$, Ham pointed out [51] that

[TABLE]

and $p$ can be expressed as functions of $E_{\rm JT}/\hbar\omega$ as shown in Fig. 27.

This approximation achieves adequate precision for the SiV system [13], where the spin-orbit coupling is much smaller than the electron-phonon coupling strength and phonon energy. However, the error of this approximation increases for the XV systems with larger spin-orbit coupling strengths. To mitigate this problem, a solution could be to separate the spin-orbit coupling term out from $H_{\rm int}$; thereby, the lowest four energy eigenstates (incorporating electron spin states) of $H_{0}=H_{\rm vib}+H_{\rm vo}+H_{\rm os}$ may form a quasi-independent subspace regarding other interaction terms associated with the electronic orbital states. In this new subspace, quench factors $p^{\prime}$ and $q^{\prime}$ can be defined similarly.

7.2.2 Refined Definition

We consider the XV system Hamiltonian in the form

[TABLE]

where $H_{0}=H_{\rm vib}+H_{\rm vo}+H_{\rm os}$, and $H^{\prime}_{\rm int}$ is the interaction Hamiltonian excluding the electron-phonon interaction and spin-orbit interaction, which can also include the electronic orbital operators $\tau_{z}^{\rm o}$, $\tau_{x}^{\rm o}$, and $\tau_{y}^{\rm o}$.

As discussed in Sec. 3.3, $H_{0}$ is two-fold degenerate. Since it commutes with the $z$-component spin operator $\sigma_{z}^{\rm s}$, the four lowest energy eigenstates can be presented as $\left\{\ket{1^{\prime}},\ket{2^{\prime}}\right\}\otimes\{\ket{\uparrow},\ket{\downarrow}\}$, where $\left\{\ket{1^{\prime}},\ket{2^{\prime}}\right\}$ are the vibronic components and $\{\ket{\uparrow},\ket{\downarrow}\}$ are the spin components, specifically, the eigenstates of $\sigma_{z}^{\rm s}$. If the energy difference between the energy of the fifth and fourth lowest energy eigenstate $\Delta_{5,4}$ is much larger than the interaction strength of $H^{\prime}_{\rm int}$, the four lowest energy eigenstates form a quasi-independent subspace.

In this subspace, we can define

[TABLE]

Then, in a suitable energy gauge, the Hamiltonian 164 can be approximated with

[TABLE]

where $\Delta_{0}$ is the ground splitting of $H_{0}$ and $\{\tau_{i}^{\rm m}\}$ ($i=x,y,z$) are the Pauli matrices in the basis $\left\{\ket{1^{\prime}},\ket{2^{\prime}}\right\}$. Note that under this refined definition, $p^{\prime}$ and $q^{\prime}$ do not obey the relationship similar to Eq. (163).

The quantity $\Delta_{0}$ can be determined in experiments. It corresponds to the ground splitting of an unstrained sample under zero magnetic field. In some literature [14], when no need to be specified, the quench factors can be absorbed into the effective interaction strengths and not shown.

To show that the quench factors $p^{\prime}$ and $q^{\prime}$ are of higher precision than $p$ and $q$ for XV species with larger spin-orbit coupling strengths, we harness the parameters of $\hbar\omega$, $E_{\rm JT}$, and $\lambda$ of various XV systems from the reported DFT calculations [13] as shown in Tab. 6). Additionally, we examine two interaction Hamiltonians, (1) $H^{\prime}_{\rm int}=\gamma_{z}\tau_{z}^{\rm o}$ for $p-p^{\prime}$ comparison and (2) $H^{\prime}_{\rm int}=\gamma_{x}\tau_{x}^{\rm o}$ for $q-q^{\prime}$ comparison. We limit our scope to

[TABLE]

Based on the quench factors $p$ and $q$, the ground splittings for the two types of $H^{\prime}_{\rm int}$ are given by

[TABLE]

respectively. The corresponding precisions of these evaluations are defined by

[TABLE]

where $\Delta_{\rm GS,num,x}$ and $\Delta_{\rm GS,num,z}$ are from the numerical diagonalization of the Hamiltonian $H_{0}+H^{\prime}_{\rm int}$, respectively. The precisions of the ground splitting evaluations based on $p^{\prime}$ and $q^{\prime}$ are defined in the same fashion and denoted by $\epsilon^{\prime}_{\rm qch,z}$ and $\epsilon^{\prime}_{\rm qch,x}$, respectively.

The comparison between the precisions of the quench factors based on the original and refined definitions for various XV system species are shown in Tab. 7. From the comparison, we see that the quench factors $p$ and $p^{\prime}$ give similar precisions. However, the precision of the quench factor $q^{\prime}$ is an order of magnitude greater than that of $q$, demonstrating an improved fit of the refined definition compared to the original one. Additionally, the precision of the quench factors $p^{\prime}$ and $q^{\prime}$ still decreases as the atomic number of X increases. This is because of the decrease of $\Delta_{5,4}$, the energy difference between the fifth and fourth lowest eigenenergies of $H_{0}$, as shown in Tab. 6.

8 Previous Studies

8.1 Results of this article

In Secs. 2.5 and 2.6, we provided order-of-magnitude estimations for the XV system parameters $\hbar\omega$, $F_{0}$, $F$, and $\lambda$, with undetermined partition factor $0<\alpha_{\rm X}<1$ and screening factor $0<\beta_{\rm X}<1$ and a relation

[TABLE]

These parameters are shown in Tab. 8.

In Sec. 4.5, we discover relations between the strain susceptibilities and the electron-displacement coupling strengths $F_{0}$ and $F$ as follows

[TABLE]

where $a_{0}=3.567\,\text{\r{A}}$ is the diamond lattice constant.

In Sec. 7.1, we show that the Jahn-Teller energy is defined as

[TABLE]

where $\mu\in\{28.0,72.6,118.7,207.2\}\,\text{Da}$ (1 Da = $1.66\times 10^{-27}$ kg) is the mass of the X nucleus.

In Sec. 3.4, we show that the electron-phonon coupling strength $J$ is given by

[TABLE]

8.2 Comparison with Hepp’s model

Hepp has pioneered modeling the SiV system [11]. His model was later extended to apply to the SnV system by Trusheim et al. [52]. This article functions as both a rectification and an expansion of Hepp’s model. Building upon first-principle reasoning, the primary accomplishment of this article is threefold.

Firstly, using group theory, Hepp analyzes the symmetries of the electronic orbital states and nuclear vibrations separately and concludes that there are four nuclear vibration modes, two sets of couples, each with $\epsilon_{g}$ symmetry, coupled to the electronic orbital doublet. Nevertheless, this article found it problematic to analyze the symmetries of the nuclear vibrations independent of the electronic orbital configurations. Instead, we classify the nuclear vibration modes according to the spatial configuration of the electronic orbital states, as presented by Eq. (27), and derive the conclusion that, in the simplest form, there are only two nuclear vibration modes dominantly coupled to the electronic orbital doublet. The four nuclear vibration modes identified by Hepp all have non-zero projections onto the two modes identified in this article, which is why they all have significant coupling to the electronic orbital doublet.

Secondly, within Hepp’s framework, the spin-orbit coupling is considered the primary interaction in terms of magnitude, while the Jahn-Teller interaction (also known as the electron-phonon interaction) is regarded as secondary. Based on his logic, the interaction Hamiltonian is expressed in the bare basis of electronic orbital and spin states, and the Jahn-Teller is presented as an interaction term. Additionally, the quench factors, as derived from the Jahn-Teller effect, only apply to terms weaker than the spin-orbit couplings, for example, the orbital gyromagnetic ratio. However, this logic can be potentially wrong, as the presumption that the spin-orbit coupling is much stronger than the Jahn-Teller interaction is not justified, and the Jahn-Teller interaction term presented in his framework is not proven to be non-zero.

What’s more, Hepp’s modeling of the Jahn-Teller effect directly contradicts the model used in Theiring’s calculation [13]. In Hepp’s model, in the presence of both the spin-orbit coupling and Jahn-Teller effect, the ground splitting of the SiV system is given by

[TABLE]

where $\Upsilon_{x}$ and $\Upsilon_{y}$ are the parameters that represent the strength of the Jahn-Teller effect. Thus, the presence of the Jahn-Teller effect increases the ground-state splitting beyond that induced by spin-orbit coupling. However, in the framework used in Theiring’s calculation, it leads to a reduction, since the Jahn-Teller effect is characterized by quench factors, which are applied to the spin-orbit coupling strengths.

Through meticulous reasoning, this article supports the model utilized in Theiring’s calculation. And, as presented in Sec. 7.2.2, we also propose a refinement on that approximate model with improved precision. In this refined model, the spin-orbit coupling and Jahn-Teller effect are considered simultaneously, and we propose refined quench factors that apply to all the other terms associated with electronic orbital states except for these two. Coincidentally, the data processing methods used in Meesala’s paper [14] align with the refined quench-factor model.

Thirdly, in Hepp’s study, he experimentally observed a branching ratio of about $\frac{1}{2}$ in the SiV emission spectrum but did not explain it. This article accomplishes this explanation, as presented in Sec. 5.4. We show that this factor originates from the Jahn-Teller effect, or the quench factors.

8.3 Partition factor and screening factor based on Thiering’s calculation

Thiering et al. [13] have calculated the absolute values of $\hbar\omega$, $E_{\rm JT}$, and $\lambda$ for the four XV systems with DFT methods. Tab. 9 shows the partition factors ($\alpha_{\rm X}$) and screening factors ($\beta_{\rm X}$) based on their calculation results.

Based on the Tab. 9, we observe larger partition factors and screening factors in the first-excited manifolds than the ground manifolds for all XV systems. This can be explained by the presence and absence of the central node in the ground and first-excited manifolds, respectively. Additionally, we note that, within the same manifold, partition factors and screening factors exhibit minimal variation among the XV species.

8.4 Lowest splitting measurement

In the Tab. 10 below, we show the measured lowest splittings in the ground and first-excited manifolds for all the XV species along with the ones calculated with the DFT method by Thiering et al. [13].

8.5 Strain Susceptibilities

Meesala et al. have claimed measuring the quenched SiV strain susceptibility parameters $d$, $f$ for both the ground and first-excited manifolds, along with two other parameters $t_{\parallel,e}-t_{\parallel,g}$ and $t_{\perp,e}-t_{\perp,g}$, the differences in parameter $t_{\parallel}$ and $t_{\perp}$ between the first-excited and ground manifolds [14]. In the data analysis section, Meesala et at haven’t taken the Jahn-Teller effect into account. However, the measured values of $d$ and $f$ can be interpreted as the quenched parameters $p^{\prime}d$, $p^{\prime}f$ based on the refined definition of the quench factors as introduced in Sec. 7.2.2 while the values of $t_{\parallel,e}-t_{\parallel,g}$ and $t_{\perp,e}-t_{\perp,g}$ do not experience any quenching effects hence being measured faithfully. The method of deriving the strain susceptibility from measuring the splitting under a controlled strain is also present in Sec. 9.

On the other hand, Guo et al. have calculated the original (un-quenched) SnV strain susceptibility parameters $d$ and $f$ with DFT methods using both PBE and SCAN functionals [15].

In the meantime, based on the relationships [Eq. (174)] between the Jahn-Teller energy and the strain susceptibilities we found in this article, we are allowed to convert Thiering’s DFT calculation results of the Jahn-Teller energies and phonon energies into the strain susceptibilities and compare them with Meesala’s and Guo’s results. We summarize these three data sets in Tab. 11 and Tab. 12.

From Tab. 11 and Tab. 12, we see that the results obtained by Meesala and Guo exhibit some agreement in both signs and the order of magnitudes, even though the XV species involved are distinct. Nonetheless, the parameters converted from Thiering’s calculation possess discrepancies with Meesala’s and Guo’s results in both signs and magnitude. Additional research endeavors are required to address this conflict.

If Thiering’s results and the connections between the Jahn-Teller energy and the strain susceptibilities outlined in this article are accurate, one possible explanation for this discrepancy lies in the intrinsic strain of the diamond sample. To explain this, when there is intrinsic strain $\gamma_{0}$ in the XV sample, the original ground splitting of the XV system is given by

[TABLE]

The artificial stress, when in certain directions, may counteract the intrinsic strain $\gamma_{0}$, resulting in diminishing the ground splitting and giving rise to opposite signs in the strain susceptibilities. Given that Meesala et al. conducted measurements on a single sample in their study, this possibility is not excluded. Furthermore, in section 9, we propose a method to evaluate the sample’s inherent strain, potentially helping to eliminate its effect.

9 Intrinsic strain measurement

Based on the refined definition of the quench factors, the Hamiltonian of an XV system under strain is given by

[TABLE]

where

[TABLE]

and $\tau_{\perp}^{\rm m}=\cos\theta\tau_{x}^{\rm m}+\sin\theta\tau_{y}^{\rm m}$ with

[TABLE]

This Hamiltonian has four dimensions and is two-fold degenerate. The splitting between the two sets of doublets is given by

[TABLE]

Additionally, when $\gamma\ll\Delta_{0}$,

[TABLE]

Based on this relation and varying the magnitude and direction of the applied stress, one may measure the quenched strain susceptibilities of the SiV system, which is what Meesala et al. obtained.

However, there is a drawback to this method. The strain of an XV system can originate from two primary sources: the intrinsic strain, which arises from adjacent lattice imperfections, and the artificial strain resulting from externally applied stress, e.g.

[TABLE]

where $\gamma_{0}$ represents the intrinsic strain and $\gamma_{\rm art}$ the artificial strain. Meesala’s calculation integrates the intrinsic strain $\gamma_{0}$ into the intrinsic splitting $\Delta_{0}$, hence not being measured.

Nevertheless, here we propose a method to measure the intrinsic strain by harnessing a magnetic field along the $y$-direction. Based on the refined quench-factor model, the Hamiltonian of an XV system in the presence of a magnetic field in the $y$-direction and strain is given by

[TABLE]

where $B_{y}$ is the strength of the magnetic field in the $y$-direction.

With a magnetic field added, the double degeneration in $H_{\rm S}$ is lifted. The four eigenenergies of $H_{\rm SM}$ are

[TABLE]

Hence, the splitting of those doublets is given by

[TABLE]

when $\gamma\ll\mu_{\rm B}B_{y}$,

[TABLE]

When no artificial strain is applied, $\gamma_{\rm art}=0$, the splitting $\Delta_{\rm SM}$ is contributed solely from the intrinsic strain $\gamma_{0}$, hence it can be measured.

10 Conclusion

We refine the theoretical framework upon previous endeavors of the XV center quantum system. Specifically, we give thorough explanations for the origin of the spin-orbit and electron-phonon interactions, and the mechanism of how the strain, electric field, and magnetic field influence the potential energy of the XV system. From those explorations, the readers may understand why the XV centers possess properties such as large ground-state splitting, robust light coupling without phonon excitation, strain tunability, resilience to charge noise, and the ability to be manipulated via magnetic fields. Original insights of this article include

• •

only two photonic modes dominantly couple to the XV orbital states,

• •

a novel symmetry, the joint-reflection symmetry, that explains the system degeneracy under a magnetic field perpendicular to the principal symmetry axis,

• •

the relationship between the electron-displacement coupling strength and the strain susceptibility,

• •

the explanation of the stark shift,

• •

the explanation of the branching ratio of light emission intensity,

• •

the examination of the quantum state addressability with the aid of the magnetic field,

• •

an improved quench-factor model with higher precision for large spin-orbit coupling strengths, and

• •

a potential proposal to measure the intrinsic strain with the aid of a magnetic field.

Besides, this article compares the existing theoretical and experimental studies, pointing out the discrepancies among them and suggesting the potential reasons for the conflicts. In summary, our studies deepen the understanding of the XV center systems and provide useful guidance to manipulate them, expanding the set of tools available for quantum information processing.

\bmhead

Acknowledgements

Parts of the research presented in Secs. 3.4, 3.5, and 5.4 were conceptualized and supervised by Dr. V. V. Dobrovitski and Dr. J. Borregaard, who also participated in the development of the methodology, analysis, investigation, and validation of the results. The author thanks Enli Chen for providing the proof of the lemma in Sec. 3.3. The author gratefully acknowledges that this work was partially supported by the collaboration between Fujitsu Limited and Delft University of Technology, co-funded by the Netherlands Enterprise Agency under project number PPS2007. The author also acknowledges funding from the NWO Gravitation Program Quantum Software Consortium (Project QSC No. 024.003.037).

Appendix A List of Used Approximations

A.1 Approximations Used in the Electron-phonon Interaction Modeling

• •

Localized Phonon Approximation: We neglect the mutual coupling between the seven nuclei within the XV system and their coupling to the remaining nuclei in the diamond bulk, resulting in seven localized independent nuclear vibrations. This approximation is based on the significant difference between the phonon frequencies of the XV system and those of the bulk, which minimizes the coupling.

• •

Imaginary Bonding Structures: The imaginary bonding structures as shown in Fig. 7 are simple models used to approximate the potential that a nucleus experiences in a certain direction. Their applicability is based on the similarity of the potential stiffness in certain directions.

• •

X-C Bond Length: The X-C bond lengths listed in Tab. 1 are based on measurements or DFT calculations for molecules or perfect crystals. Those values may differ from those in the XV color center as a defect in the diamond.

• •

Potential-Width Parameters: The potential-width parameters for the X-C bonds, $\alpha_{\rm X}$, are estimated using a guess formula (Eq. 10). This approximation assumes that the width parameter $\alpha$ for each X-C bond and C-C bond scales proportionally with its bond dissociation energy $E_{\rm 0}$ and inversely with its bond length $d_{\rm 0}$.

• •

Nearest Neighbor Interactions: We neglect the interaction of the $i$-th X-C bond orbital with atomic cores other than the $i$-th carbon atomic core. This approximation is justified by the larger distances between the $i$-th X-C bond and the other atomic cores, as well as the shielding effect provided by the electronic clouds between them.

• •

Point-Charge Model: This model approximates the middle portion of the spatially extended X–C bond orbital and the carbon atomic core as point charges with effective charges and positions. The validity of this approximation relies on the fact that the middle portion of the X–C bond orbital is distant enough from both the X and carbon nuclei. Here, the middle portion refers to the portion of the hole population that is not stationary with respect to its nearest nucleus.

• •

Small nuclear vibration amplitudes This approximation is based on the calculation presented in Sec. 2.5.2, where it shows that the typical value of the nuclear vibration amplitudes is more than 10 times smaller than the inverse of the potential-width parameter $\alpha$ and the bond length. This approximation is the basis for approximating the movement of the nuclei as harmonic oscillators. And, combined with the adaptability of the hole orbital, it is also the basis to neglect the effect caused by the nuclear displacement in the $x_{i}$ and $y_{i}$ directions because the net effect is approximately a rotation of the X–C bond.

• •

Linear Coulomb Coupling: We assume that the Coulomb energy is proportional to the atomic displacement, with the proportionality factor given by the first-order derivative of the Coulomb energy, evaluated at the equilibrium nuclear configuration, with respect to the displacements of nuclei. The validity of this approximation relies on the ”Small nuclear vibration amplitudes” approximation.

• •

Isolated X–C bond: When determining the effective position of the middle portion of the hole cloud, as approximated as a point charge, we assume the X–C bond orbital takes the shape as in an isolated X–C bond, neglecting the influence on the shape of the orbital from the other five X–C bonds in the XV system. Specifically, the presence (absence) of a node in the location of the X nucleus when the XV system is in its ground (first-excited) manifold changes the shape of the X–C bond orbital, hence the effective position of the middle portion of the hole cloud. Consequently, our model gives orders-of-magnitude estimations of the electron-phonon coupling strength for both the ground and first-excited manifolds without specifying the difference between the two cases.

A.2 Approximations Used in the Spin-Orbit Interaction Modeling

• •

Isolated X atom: We assume an isolated X atom when calculating the spin-orbit coupling strength of the hole in the XV system. This approximation is valid because the majority of the hole population is confined within the region enclosed by the six carbon nuclei of the XV system. Adopting the classical picture, the hole orbits around only the X nucleus and its inner-shell electrons. The electric fields generated by the charge outside the hole’s trajectory have no divergence within the area enclosed by the hole’s trajectory, hence contribute zero net spin-orbit coupling strength.

• •

Neglected contribution from the carbon $sp^{3}$ orbitals: The hole is distributed across the six X–C bonds, where each bond involves an $sp^{3}d^{2}$ orbital on the X atom and an $sp^{3}$ orbital on the carbon atom. However, we neglect the contribution of the spin-orbit coupling strength from the population located in the $sp^{3}$ orbitals for two reasons: 1) the spin-orbit coupling strength is reciprocal cubic decay in distance from the X nucleus and (2) the population in the $sp^{3}$ orbitals suffers a more substantial screening effect on the electric field of the X nucleus than that in the $sp^{3}d^{2}$ orbitals. Regarding the ratio between the population in the $sp^{3}d^{2}$ and $sp^{3}$ orbitals, it is smaller in the ground manifold than in the first-excited manifold because of the presence and absence of the node in the electronic orbital eigenstates, respectively. However, we assume this ratio to be one for both cases as a rough estimation.

Appendix B Deduction of the Matrix Elements of the displacement operator111The method in this appendix is based on Proc. R. Soc. Lond. A 1958 244, 1-16[36].

This appendix aims to project the nuclear displacement operators $\hat{Q}_{x}$ and $\hat{Q}_{y}$ onto the basis set $\{|nm\rangle\}$, as introduced in Sec. 3.4.

To improve the readability, we replace the nuclear displacement operators $\hat{Q}_{x}$ and $\hat{Q}_{y}$ with simpler notations $\hat{x}$ and $\hat{y}$ in the following. Please don’t mix them with the hole orbital displacements in the main text.

To achieve the goal, we need to calculate the matrix elements $\langle n^{\prime}m^{\prime}|x|nm\rangle$ and $\langle n^{\prime}m^{\prime}|y|nm\rangle$. Instead of performing cumbersome integrals with the analytical expressions of $|nm\rangle$, we adopt an algebraic approach, deducing the result from three fundamental facts:

[TABLE]

In Eq. (190), the potential energy of the harmonic oscillator is known to account for half of the total energy. At first glance, the connection between these equations and the desired matrix elements might seem indirect. However, by carefully applying logical reasoning, we bridge this gap and derive the desired results. Let the logic do the magic.

First, we aim to identify which matrix elements are non-zero. To achieve this, we begin by examining a seemingly unrelated calculation:

By evaluating the commutator $\left[x,H_{\textrm{v}}\right]$, we derive the following relation:

[TABLE]

If we replace $x$ in the equation above with

[TABLE]

the relation remains valid. By sandwiching this relation with $\langle n^{\prime}m^{\prime}|$ and $|nm\rangle$ and defining

[TABLE]

we arrive at:

[TABLE]

Since this relation must always hold, it implies that for $Q_{\pm\>nm}^{\ \ n^{\prime}m^{\prime}}$ to be non-zero, the following condition must be satisfied:

[TABLE]

This result indicates that we only need to compute two types of matrix elements: $Q_{\pm\>nm}^{\ \ n+1\;m^{\prime}}$ and $Q_{\pm\>nm}^{\ \ n-1\;m^{\prime}}$.

Next, by calculating the commutators $[x,L_{z}]$ and $[y,L_{z}]$, we derive the following relation:

[TABLE]

Sandwiching this relation with $\langle n^{\prime}m^{\prime}|$ and $|nm\rangle$ and applying similar reasoning as above, we obtain:

[TABLE]

In summary, only the following four types of matrix elements can potentially be non-zero:

[TABLE]

Furthermore, noting that, based on the definition of the sandwiching operation,

[TABLE]

only two of the elements, $Q_{+\>nm}^{\ \ n\pm 1\;m+1}$, need to be calculated. In fact, we can write:

[TABLE]

because, based on the analytical expression of the eigenwavefunctions of $|nm\rangle$ [35], the matrix elements are all real. This property will be used again in subsequent derivations.

To proceed, we utilize the relation $x^{2}+y^{2}=Q_{-}Q_{+}$. Sandwiching this relation with $\langle nm|$ and $|nm\rangle$, we find:

[TABLE]

Recalling Eq. (190), we derive the first key relation for determining $Q_{+\>nm}^{\ \ n\pm 1\;m+1}$:

[TABLE]

Next, we derive a second relation by examining the following:

[TABLE]

which can be proven using the commutators employed earlier. The first few essential steps are:

[TABLE]

Sandwiching the equation with $\langle nm|$ and $|nm\rangle$, we obtain:

[TABLE]

Combining Eqs. (203) and (205), we solve for the matrix elements:

[TABLE]

For clarity, $Q_{+}$ and $Q_{-}$ can be expressed as the sum of submatrices with the bases explicitly denoted:

[TABLE]

[TABLE]

Note that for cases where $m=\pm n$, $|nm\rangle$ may couple to undefined states, such as $|n-1\ n+1\rangle$ or $|n-1\ -n-1\rangle$. However, since the corresponding matrix elements are zero, this inclusion does not cause confusion.

Returning to the original question of expressing $x$ and $y$ in matrix form, the operators $x$ and $y$ are given by $x=\frac{1}{2}\left(Q_{+}+Q_{-}\right)$ and $y=\frac{1}{2\mathrm{i}}\left(Q_{+}-Q_{-}\right)$. Thus, we can write:

[TABLE]

[TABLE]

For the second-order Jahn-Teller effect, we need to compute $x^{2}-y^{2}$ and $2xy$. The results are as follows:

[TABLE]

[TABLE]

For the Jahn-Teller operators $Q_{\pm\pm}\coloneqq\frac{1}{\sqrt{2}}\left[\left(x^{2}-y^{2}\right)\mp\mathrm{i}2xy\right]$, the expressions are:

[TABLE]

[TABLE]

With these results, all matrix elements required for this thesis have been derived.