"""
This module contains P1 class for simulating P1 centers in diamond,
being a subclass of QSys.
"""
from typing import Literal
import numpy as np
from qutip import Qobj, jmat, qeye, tensor
from .constants import gamma_e, gamma_N14, gamma_N15
from .qsys import QSys
__all__ = ["P1"]
####################################################################################################
[docs]
class P1(QSys):
"""
P1 class contains attributes and methods to simulate the neutral substitutional nitrogen center
in diamond.
Due to Jahn-Teller effect this defect has a C_{3v} point symmetry, which means there are four
isolated orientation families of P1 centers in diamond. In simulations the lab reference frame
is chosen in such a way that it is a principal axis system (PAS) for one of the four families
(rotation index 0). Hamiltonians of other orientation families are obtained by the transition
from PAS to the corresponding orientation using rotation matrices.
Attributes
----------
B0 : float
Magnetic field intensity.
rot_index : int
Rotation index, integer between 0 and 3.
R : list
List of rotation matrices to go from the PAS (principal axis system) frame to the lab frame.
N : int
Nitrogen isotope, or 0 for no nuclear spin.
units_B0 : str
Units of the magnetic field (T, mT or G)
theta : float
Angle of the static magnetic field vector with respect to the z axis
phi_r : float
Azimutal angle of the static magnetic field vector within the xy plane
theta_1 : float
Angle of the MW magnetic field vector with respect to the z axis
phi_r_1 : float
Azimutal angle of the MW magnetic field vector within the xy plane
units_angles : str
Units of the angles (deg or rad)
B0_vector : np.ndarray
Unit vector of the static magnetic field in the lab frame
B1_vector : np.ndarray
Unit vector of the MW magnetic field in the lab frame
h1 : Qobj
Standard microwave Hamiltonian for the P1 center corresponding to the electronic spin
transitions.
eigenstates : np.ndarray
Array containing the eigenstates of the Hamiltonian.
dim_add_spin : int
Dimension of the added spin, if any.
Methods
-------
_rot_pas_to_lab
Defines the rotation matrices to go from the PAS (principal axis system) frame to the lab
frame.
electron_zeeman
Calculates the electron Zeeman Hamiltonian term.
hyperfine
Calculates the hyperfine Hamiltonian term.
quadrupole
Calculates the quadrupole Hamiltonian term.
nuclear_zeeman
Calculates the nuclear Zeeman Hamiltonian term.
_set_h1
Sets the standard microwave Hamiltonian for the P1 center corresponding to the electronic
spin transitions.
add_spin
Adds an extra spin to the system and updates the microwave Hamiltonian accordingly.
Notes
-----
The P1 class inherits the methods and attributes from the QSys class.
"""
def __init__(
self,
B0: float,
*,
rot_index: Literal[0, 1, 2, 3],
N: Literal[15, 14, 0] | None = None,
rho0: Qobj | np.ndarray | None = None,
c_ops: Qobj | list[Qobj] | None = None,
observable: Qobj | list[Qobj] | None = None,
units_B0: Literal["T", "mT", "G"] = "mT",
theta: float = 0.0,
phi_r: float = 0.0,
theta_1: float = 90.0,
phi_r_1: float = 0.0,
units_angles: Literal["rad", "deg"] = "deg",
) -> None:
"""
Constructor for the P1 class.
Takes the nitrogen isotope, the rotation index, the magnetic field intensity as inputs and
calculates the Hamiltonian with all relevant attributes.
Parameters
----------
B0 : float
Magnetic field intensity.
rot_index : int
Rotation index, integer between 0 and 3.
N : 15 | 14 | 0 | None
Nitrogen isotope, or 0 for no nuclear spin
rho0 : Qobj | array
Initial state of the system. Can be a Qobj, an array or an index number indicating the
system eigenstates
c_ops : Qobj | list(Qobj)
List of collapse operators
observable : Qobj | list(Qobj)
Observable to be measured
units_B0 : str
Units of the magnetic field (T, mT or G)
theta : float
Angle of the static magnetic field vector with respect to the z axis
phi_r : float
Azimutal angle of the static magnetic field vector within the xy plane
theta_1 : float
Angle of the MW magnetic field vector with respect to the z axis
phi_r_1 : float
Azimutal angle of the MW magnetic field vector within the xy plane
units_angles : str
Units of the angles (deg or rad)
"""
self.B0, self.units_B0 = self._check_B0(B0, units_B0)
self.theta, self.phi_r, self.units_angles = self._check_angles(theta, phi_r, units_angles)
self.theta_1, self.phi_r_1, _ = self._check_angles(theta_1, phi_r_1, units_angles)
self.B0_vector = np.array(
[
[
np.sin(self.theta) * np.cos(self.phi_r),
np.sin(self.theta) * np.sin(self.phi_r),
np.cos(self.theta),
]
]
)
self.B1_vector = np.array(
[
[
np.sin(self.theta_1) * np.cos(self.phi_r_1),
np.sin(self.theta_1) * np.sin(self.phi_r_1),
np.cos(self.theta_1),
]
]
)
if rot_index not in range(4):
raise ValueError(
f"Invalid value for rotation index r. Expected an integer between 0 and 3, got {rot_index}."
)
self.rot_index = rot_index
self._rot_pas_to_lab()
self.N = N
if N == 14:
H0 = (
self.electron_zeeman()
+ self.hyperfine()
+ self.quadrupole()
+ self.nuclear_zeeman()
)
elif N == 15:
H0 = self.electron_zeeman() + self.hyperfine() + self.nuclear_zeeman()
elif N == 0 or N is None:
H0 = self.electron_zeeman()
else:
raise ValueError(
f"Invalid value for nitrogen isotope N. Expected either 14 or 15, got {N}."
)
self.eigenstates = np.array([psi * psi.dag() for psi in H0.eigenstates()[1]])
# check if the obervable is an integer or a list of integers
# if yes then takes the eigenstates' observables corresponding to the integer
if isinstance(observable, int) and observable in range(len(self.eigenstates)):
observable = self.eigenstates[observable]
# for one single observable, the initial state is assigned to it
if rho0 is None:
rho0 = observable
elif isinstance(observable, (list, np.ndarray)) and all(
isinstance(obs, int) and obs in range(len(self.eigenstates)) for obs in observable
):
observable = [self.eigenstates[obs] for obs in observable]
super().__init__(H0, rho0=rho0, c_ops=c_ops, observable=observable, units_H0="MHz")
self._set_h1()
def _rot_pas_to_lab(self) -> None:
"""
Defines the rotation matrices to go from the PAS (principal axis system)
frame to the lab frame.
Start by defining the basis vectors of the PAS frame in the lab frame.
Then generates the R list containing the rotation matrices.
"""
basis = []
basis.append(
[
1 / np.sqrt(6) * np.array([[1, 1, -2]]).T,
1 / np.sqrt(2) * np.array([[1, -1, 0]]).T,
1 / np.sqrt(3) * np.array([[-1, -1, -1]]).T,
]
)
basis.append(
[
1 / np.sqrt(6) * np.array([[2, 1, 1]]).T,
1 / np.sqrt(2) * np.array([[0, 1, -1]]).T,
1 / np.sqrt(3) * np.array([[-1, 1, 1]]).T,
]
)
basis.append(
[
1 / np.sqrt(6) * np.array([[1, 2, 1]]).T,
1 / np.sqrt(2) * np.array([[-1, 0, 1]]).T,
1 / np.sqrt(3) * np.array([[1, -1, 1]]).T,
]
)
basis.append(
[
1 / np.sqrt(6) * np.array([[1, 1, 2]]).T,
1 / np.sqrt(2) * np.array([[1, -1, 0]]).T,
1 / np.sqrt(3) * np.array([[1, 1, -1]]).T,
]
)
R_1 = np.c_[basis[0][0], basis[0][1], basis[0][2]]
R_2 = np.c_[basis[1][0], basis[1][1], basis[1][2]]
R_3 = np.c_[basis[2][0], basis[2][1], basis[2][2]]
R_4 = np.c_[basis[3][0], basis[3][1], basis[3][2]]
R_11 = R_1.T @ R_1
R_12 = R_1.T @ R_2
R_13 = R_1.T @ R_3
R_14 = R_1.T @ R_4
self.R = [R_11, R_12, R_13, R_14]
[docs]
def electron_zeeman(self) -> Qobj:
"""
Electron Zeeman Hamiltonian term,
rotated according to the rotation index r.
Returns
-------
Zeeman Hamiltonian : Qobj
"""
H_ez = (
gamma_e
* self.B0
* (
(self.R[self.rot_index] @ self.B0_vector.T)[0][0] * jmat(1 / 2, "x")
+ (self.R[self.rot_index] @ self.B0_vector.T)[1][0] * jmat(1 / 2, "y")
+ (self.R[self.rot_index] @ self.B0_vector.T)[2][0] * jmat(1 / 2, "z")
)
)
return self._tensor_product_N(H_ez, self.N)
[docs]
def hyperfine(self) -> Qobj:
"""
Get the hyperfine term
Returns
-------
Hyperfine Hamiltonian : Qobj
"""
if self.N == 14:
return 114.03 * tensor(jmat(1 / 2, "z"), jmat(1, "z")) + 81.32 * (
tensor(jmat(1 / 2, "x"), jmat(1, "x")) + tensor(jmat(1 / 2, "y"), jmat(1, "y"))
)
elif self.N == 15:
return -159.73 * tensor(jmat(1 / 2, "z"), jmat(1 / 2, "z")) + -113.84 * (
tensor(jmat(1 / 2, "x"), jmat(1 / 2, "x"))
+ tensor(jmat(1 / 2, "y"), jmat(1 / 2, "y"))
)
else:
raise ValueError(
f"Invalid value for nitrogen isotope N. Expected either 14 or 15, got {self.N}."
)
[docs]
def quadrupole(self) -> Qobj | int:
"""
Get the quadrupole term
Returns
-------
Quadrupole Hamiltonian : Qobj
"""
if self.N == 14:
return -3.97 * tensor(qeye(2), jmat(1, "z") ** 2)
elif self.N in (0, 15):
return 0
else:
raise ValueError(
f"Invalid value for nitrogen isotope N. Expected either 14 or 15, got {self.N}."
)
[docs]
def nuclear_zeeman(self) -> Qobj | None:
"""
Nuclear Zeeman Hamiltonian term
Returns
-------
Zeeman Hamiltonian : Qobj
"""
if self.N == 14:
return (
-gamma_N14
* self.B0
* tensor(
qeye(2),
(self.R[self.rot_index] @ self.B0_vector.T)[0][0] * jmat(1, "x")
+ (self.R[self.rot_index] @ self.B0_vector.T)[1][0] * jmat(1, "y")
+ (self.R[self.rot_index] @ self.B0_vector.T)[2][0] * jmat(1, "z"),
)
)
elif self.N == 15:
return (
-gamma_N15
* self.B0
* tensor(
qeye(2),
(self.R[self.rot_index] @ self.B0_vector.T)[0][0] * jmat(1 / 2, "x")
+ (self.R[self.rot_index] @ self.B0_vector.T)[1][0] * jmat(1 / 2, "y")
+ (self.R[self.rot_index] @ self.B0_vector.T)[2][0] * jmat(1 / 2, "z"),
)
)
else:
return None
def _set_h1(self) -> None:
"""
Sets the standard control Hamiltonian for the P1 center corresponding to the electronic spin
transitions.
"""
h1_e = 2 * (
(self.R[self.rot_index] @ self.B1_vector.T)[0][0] * jmat(1 / 2, "x")
+ (self.R[self.rot_index] @ self.B1_vector.T)[1][0] * jmat(1 / 2, "y")
+ (self.R[self.rot_index] @ self.B1_vector.T)[2][0] * jmat(1 / 2, "z")
)
if self.N == 15:
self.h1 = tensor(h1_e, qeye(2))
elif self.N == 14:
self.h1 = tensor(h1_e, qeye(3))
elif self.N == 0 or self.N is None:
self.h1 = h1_e
else:
raise ValueError(f"Invalid value for Nitrogen. Expected either 14 or 15, got {self.N}.")
[docs]
def add_spin(self, H_spin: Qobj) -> None:
"""
Overwrites the parent class method by calling it and updating the h1 attribute
Parameters
----------
H_spin : Qobj
Hamiltonian of the extra spin
"""
super().add_spin(H_spin)
self.h1 = tensor(self.h1, qeye(self.dim_add_spin))