# NV Center in Diamond

In this tutorial we will go over the main steps of running CCE calculations for the NV center in diamond with the **PyCCE** module. Those include:

Generating the spin bath using the

`pycce.BathCell`

instance.Setting up properties of the

`pycce.Simulator`

instance.Running the calculations with the

`Simulator.compute`

function.

We will compute the Hahn-echo coherence function (with decoupling \(\pi\)-pulse applied) using the following available methods:

Conventional CCE.

Generalized CCE (gCCE).

gCCE with Monte-Carlo bath sampling.

Finally, we will run a simulation on how different bath polarization will impact Hahn-echo signal.

```
```

```
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys
import pycce as pc
import ase
from mpl_toolkits import mplot3d
seed = 8805
np.random.seed(seed)
np.set_printoptions(suppress=True, precision=5)
```

## Generate nuclear spin bath

Building a supercell of nuclear spins from the `ase.Atoms`

object.

### Build BathCell

To generate cell it from `ase.atoms`

object, use classmethod `BathCell.from_ase`

.

```
```

```
from ase.build import bulk
# Generate unitcell from ase
diamond = bulk('C', 'diamond', cubic=True)
diamond = pc.read_ase(diamond)
```

The following attributes are created with this initiallization:

`.cell`

is ndarray containing information of lattice vectors. Each**column**is a lattice vector in cartesian coordinates.`.atoms`

is a dictionary with keys corresponding to the atom name, and each item is a list of the coordinates in cell coordinates.

```
```

```
print('Cell\n', diamond.cell)
print('\nAtoms\n', diamond.atoms)
```

```
Cell
[[3.57 0. 0. ]
[0. 3.57 0. ]
[0. 0. 3.57]]
Atoms
defaultdict(<class 'list'>, {'C': [array([0., 0., 0.]), array([0.25, 0.25, 0.25]), array([0. , 0.5, 0.5]), array([0.25, 0.75, 0.75]), array([0.5, 0. , 0.5]), array([0.75, 0.25, 0.75]), array([0.5, 0.5, 0. ]), array([0.75, 0.75, 0.25])]})
```

### Populate BathCell with isotopes

The **PyCCE** package uses EasySpin database of the concentrations of all common stable isotopes with non-zero spin, however the user can proide custom concentrations.

Use function `BathCell.add_isotopes`

to add one (or several) isotopes of the element. Each isotope is initiallized with `tuple`

containing name of the isotope and its concentration.

Name of the isotope includes the number and element symbol, provided in the `atoms`

object. As an output, the `BathCell.add_isotopes`

method returns view on dictionary `BathCell.isotopes`

which can be modified directly. Structure of the dictionary-like object:

```
{element_1: {isotope_1: concentration, isotope_2: concentration},
element_2: {isotope_3: concentration ...}}
```

```
```

```
# Add types of isotopes
diamond.add_isotopes(('13C', 0.011))
```

```
```

```
defaultdict(dict, {'C': {'13C': 0.011}})
```

Isotopes may also be directly added to `BathCell.isotopes`

. For example, below we are adding an isotope without the nuclear spin:

```
```

```
diamond.isotopes['C']['14C'] = 0.001
```

### Set z-direction of the bath (optional)

In the `Simulator`

object everything is set in \(S_z\) basis. When the quantization axis of the defect does not allign with the (0, 0, 1) direction of the crystal axis, the user needs to define the axis.

If one wants to specify the complete rotation of cartesian axes, one can provide a rotation matrix to rotate the cartesian reference frame with respect to the cell coordinates by calling the `BathCell.rotate`

method.

```
```

```
# set z direction of the defect
diamond.zdir = [1, 1, 1]
```

### Generate spin bath

To generate the spin bath, use the `BathCell.gen_supercell`

method. First argument is the linear size of the supercell (minimum distance between any two faces of the supercell is equal to or larger than this parameter). Additional keyword arguments are `remove`

and `add`

.

`remove`

takes a `tuple`

or `list of tuples`

as an argument. First element of each tuple is the name of the **atom** at that location, second element - coordinates in unit cell coordinates. If such atoms are found in the supercell, they are removed from it.

`add`

takes a `tuple`

or `list of tuples`

as an argument. First element of each tuple is the name of the **isotope** at that location, second element - coordinates in unit cell coordinates. Each of the specified isotopes will be added in the final supercell at specified locations.

```
```

```
# Add the defect. remove and add atoms at the positions (in cell coordinates)
atoms = diamond.gen_supercell(200, remove=[('C', [0., 0, 0]),
('C', [0.5, 0.5, 0.5])],
add=('14N', [0.5, 0.5, 0.5]),
seed=seed)
```

```
/home/onizhuk/midway/codes_development/pyCCE/pycce/bath/array.py:222: UserWarning: Spin type for 14C was not provided and was not found in common isotopes.
obj[n] = array[n]
```

Note, that because the `14C`

isotope doesn’t have a spin, **PyCCE** does not find it in common isotopes, and raises a warning. We have to provide `SpinType`

for it separately, or define the properties as follows:

```
```

```
atoms['14C'].gyro = 0
atoms['14C'].spin = 0
```

## BathArray Structure

The bath spins are stored in the `BathArray`

object - a subclass of `np.ndarray`

with fixed datastructure:

`N`

field`dtype('<U16')`

contains the names of bath spins.`xyz`

field`dtype('<f8', (3,))`

contains the positions of bath spins (in A).`A`

field`dtype('<f8', (3, 3))`

contains the hyperfine coupling of bath spins (in kHz).`Q`

field`dtype('<f8', (3, 3))`

contains the quadrupole tensor of bath spins (in kHz) (Relevant for spin >= 1).

All of the fields are accesible as attributes of `BathArray`

. Additionally, the subarrays of the specific spins are accessible with their name as indicated above.

Upon generation of the array from the cell, the `Q`

and `A`

fields are empty. The Hyperfine couplings will be automatically computed by the `Simulator`

object, however the quadrupole couplings must be set by the user.

The additional attributes allow one to access `SpinType`

properties:

`name`

returns the spin name or array of spin names;`spin`

returns the value of the spin or array of ones;`gyro`

returns gyromagnetic ratios of the spins;`q`

returns quadrupole constants of the spins;`h`

returns a dictionary with user-defined additions to the Hamiltonian.`detuning`

returns detunings of the spins (See definition below).

For example, below we print out the attributes of the first two spins in the `BathArray`

.

```
```

```
print('Names\n', atoms[:2].N)
print('\nCoordinates\n', atoms[:2].xyz)
print('\nHyperfine tensors\n', atoms[:2].A)
print('\nQuadrupole tensors\n',atoms[:2].Q)
```

```
Names
['13C' '13C']
Coordinates
[[-13.97678 -1.48178 -92.75132]
[ 27.89939 42.17939 -45.86038]]
Hyperfine tensors
[[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]]
Quadrupole tensors
[[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]
[[0. 0. 0.]
[0. 0. 0.]
[0. 0. 0.]]]
```

The properties of spin types (gyromagnetic ratio, quadrupole moment, etc) are stored in the `BathArray.types`

attribute, which is an instance of `SpinDict`

containing `SpinType`

classes. For most known isotopes `SpinType`

can be found in the `pycce.common_isotopes`

dictionary, and is set by default (including electron spin-1/2, which is denoted by setting `N = e`

). The user can add additional `SpinType`

objects, by calling `BathArray.add_type`

method or setting elements of
`SpinDict`

directly. For details of the first approach see documentation of `SpinDict.add_type`

method.

The direct setting of types is rather simple. The user can set elements of `SpinDict`

with `tuple`

, containing:

**(spin, gyromagnetic ratio, quadrupole moment *(optional)*, detuning *(optional)*, )**

**OR**

**(isotope, spin, gyromagnetic ratio, quadrupole moment *(optional)*, detuning *(optional)*, )**

where:

**isotope**(*str*) is the name of the given spin (same one as in`N`

field of`BathArray`

) to define new`SpinType`

object. The key of`SpinDict`

**has**to be the correct name of the spin (“isotope” field in the tuple).**spin**(*float*) is the total spin of the given bath spin.**gyromagnetic ratio**(*float*) is the gyromagnetic ratio of the given bath spin.**quadrupole moment**(*float*) is the quadrupole moment of the given bath spin. Relevant only when electric field gradient are used to generate quadrupole couplings for spins, stored in the`BathArray`

, with`BathArray.from_efg`

method.**detuning**(float) is an additional energy splitting for model spins, included as an extra \(+\omega \hat S_z\) term in the Hamiltonian, where \(\omega\) is the detuning.

Units of gyromagnetic ratio are rad / ms / G, quadrupole moments are given in barn, detunings are given in kHz.

```
```

```
# Several ways to set SpinDict elements
atoms.types['14C'] = 0, 0, 0
atoms.types['Y'] = ('Y', 0, 0, 0)
atoms.types['A'] = pc.SpinType('A', 0, 0, 0)
print(atoms.types)
```

```
SpinDict(13C: (0.5, 6.7283), 14N: (1.0, 1.9338, 0.0204), 14C: (0.0, 0.0000), ...)
```

## Simulator class

The parameters of the CCE simulator engine.

Main parameters to consider:

`spin`

— Either instance of the`CenterArray`

or float - total spin of the central spin (assuming one central spin).`bath`

— spin bath in any specified format. Can be either:Instance of

`BathArray`

class;`ndarray`

with`dtype([('N', np.unicode_, 16), ('xyz', np.float64, (3,))])`

containing names of bath spins (same ones as stored in self.ntype) and positions of the spins in angstroms;The name of the

`.xyz`

text file containing 4 columns: name of the bath spin and xyz coordinates in A.`r_bath`

— cutoff radius around the central spin for the bath.`order`

— maximum size of the cluster.`r_dipole`

— cutoff radius for the pairwise distance to consider two nuclear spins to be connected.`magnetic_field`

— applied magnetic field. Can also be provided during the simulation run.`pulses`

— number of pulses in Carr-Purcell-Meiboom-Gill (CPMG) sequence or the pulse sequence itself.

For the full description see the documentation of the `Simulator`

object.

First we setup a “mock” instance of `Simulator`

to visualize the smaller part of the bath around the central spin.

```
```

```
# Setting the runner engine
mock = pc.Simulator(spin=1, position=[0,0,0],
bath=atoms, r_bath=20,
r_dipole=6, order=3)
```

During the initiallization, depending on the provided keyword arguments several methods may be called:

`Simulator.read_bath`

is called if keyword`bath`

is provided. It may take several additional arguments:`r_bath`

- cutoff distance from the qubit for the bath.`skiprows`

- if`bath`

is provided as`.xyz`

file, this argument tells how many rows to skip when reading the file.`external_bath`

-`BathArray`

instance, which contains bath spins with pre defined hyperfines to be used.`hyperfine`

- defines the way to compute hyperfine couplings. If it is not given and`bath`

doesn’t contain any predefined hyperfines (`bath['A'].any() == False`

) the point dipole approximation is used. Otherwise it can be an instance of`pc.Cube`

object, or callable with signature`func(coord, gyro, central_gyro)`

, where`coord`

is an array of the bath spin coordinate, gyro is the gyromagnetic ratio of bath spin, central_gyro is the gyromagnetic ratio of the central bath spin.`types`

- instance of`SpinDict`

or input to create one.`error_range`

- maximum allowed distance between positions in`bath`

and`external_bath`

for two spins to be considered the same.`ext_r_bath`

- cutoff distance from the qubit for the`external_bath`

. Useful if`external_bath`

has very assymetric shape and user wants to keep the precision level of the hyperfine at different distances consistent.`imap`

- instance of the`pc.InteractionMap`

class, which contain tensor of bath spin interactions. If not provided, interactions between bath spins are assumed to be the same as one of point dipoles.

Generates `BathArray`

object with hyperfine tensors to be used in the calculation.

`Simulator.generate_clusters`

is called if`order`

and`r_dipole`

are provided. It produces`dict`

object, which contains the indexes of the bath spins in the clusters.

We implemented the following procedure to determine the clusters:

Each bath spin \(i\) forms a cluster of one. Bath spins \(i\) and \(j\) form cluster of two if there is an edge between them (distance \(d_{ij} \le\) `r_dipole`

). Bath spins \(i\), \(j\), and \(k\) form a cluster of three if enough edges connect them (e.g., there are two edges \(ij\) and \(jk\)) and so on. In general, we assume that spins \(\{i..n\}\) form clusters if they form a connected graph. Only clusters up to the size indicated by the `order`

parameter (equal to CCE order) are included.

We use `matplotlib`

to visualize the spatial distribution of the spin bath. The grey lines show connected pairs of nuclear spins, red dashed lines show clusters of three. You can try to increase `r_dipole`

, `r_bath`

parameters, or increase `order`

and visuallize.

```
```

```
# add 3D axis
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(projection='3d')
# We want to visualize the smaller bath
data = mock.bath
# First plot the positions of the bath
colors = np.abs(data.A[:,2,2]/data.A[:,2,2].max())
ax.scatter3D(data.x, data.y, data.z, c=colors, cmap='viridis');
# Plot all pairs of nuclear spins, which are contained
# in the calc.clusters dictionary under they key 2
for c in mock.clusters[2]:
ax.plot3D(data.x[c], data.y[c], data.z[c], color='grey')
# Plot all triplets of nuclear spins, which are contained
# in the calc.clusters dictionary under they key 3
for c in mock.clusters[3]:
ax.plot3D(data.x[c], data.y[c], data.z[c], color='red', ls='--', lw=0.5)
ax.set(xlabel='x (A)', ylabel='y (A)', zlabel='z (A)');
```

Now we setup `Simulator`

object for the actual simulation.

```
```

```
# Parameters of CCE calculations engine
# Order of CCE aproximation
order = 2
# Bath cutoff radius
r_bath = 40 # in A
# Cluster cutoff radius
r_dipole = 8 # in A
```

We will use the `CenterArray`

object to store the properties of the central spin, however for simple usecases one can provide the corresponding keywords to the `Simulator`

object directly (see examples below).

```
```

```
# position of central spin
position = [0, 0, 0]
# Qubit levels (in Sz basis)
alpha = [0, 0, 1]; beta = [0, 1, 0]
# ZFS Parametters of NV center in diamond
D = 2.88 * 1e6 # in kHz
E = 0 # in kHz
nv = pc.CenterArray(spin=1, position=position, D=D, E=E, alpha=alpha, beta=beta)
```

The code already knows the properties of the most common nuclear spins and of elecron spin (accessible under the name `'e'`

), however the user can provide their own by calling `BathArray.add_type`

method. The way to initiallize `SpinType`

objects is the same as in `SpinDict`

above.

```
```

```
# The code already knows most exsiting isotopes.
# Bath spin types
# name spin gyro quadrupole (for s>1/2)
spin_types = [('14N', 1, 1.9338, 20.44),
('13C', 1 / 2, 6.72828),
('29Si', 1 / 2, -5.3188),]
atoms.add_type(*spin_types)
```

### Setting the `Simulator`

object

All of the kwargs can be provided at the moment of creation. If all of the kwargs are provided, several methods of the `Simulator`

class are called:

`Simulator.read_bath`

;`Simulator.generate_clusters`

.

The details are available in the `Simulator`

methods description.

```
```

```
# Setting the runner engine
calc = pc.Simulator(spin=nv, bath=atoms,
r_bath=r_bath, r_dipole=r_dipole, order=order)
```

Taking advantage of subclassing `np.ndarray`

we can change *in situ* the quadrupole tensor of the Nitrogen nuclear spin.

```
```

```
nspin = calc.bath
# Set model quadrupole tensor at N atom
quad = np.asarray([[-2.5, 0, 0],
[0, -2.5, 0],
[0, 0, 5.0]]) * 1e3 * 2 * np.pi
nspin['Q'][nspin['N'] == '14N'] = quad
```

Note, that we need to apply the boolean mask **second** because of how structured arrays work.

### Compute coherence function with conventional CCE

The general interface to compute any property with PyCCE is implemented through the `Simulator.compute`

method. It takes two keyword arguments to determine which quantity to compute and how:

`method`

can take ‘cce’ or ‘gcce’ values, and determines which method to use - conventional or generalized CCE.`quantity`

can take ‘coherence’ or ‘noise’ values, and determines which quantity to compute - coherence function or autocorrelation function of the noise.

Each of the methods can be performed with Monte Carlo bath state sampling (if `nbstates`

keyword is non zero) and with interlaced averaging (If `interlaced`

keyword is set to `True`

).

In the first example we use the conventional CCE method without Monte Carlo bath state sampling. In the conventional CCE method the Hamiltonian is projected on the qubit levels, and the coherence is computed from the overlap of the bath evolution, entangled with two different qubit states.

The conventional CCE requires one argument:

`timespace`

— time points at which the coherence function is computed.

Additionally, one can provide the following arguments now, instead of when initiallizing `Simulator`

object:

`pulses`

— number of pulses in CPMG sequence (0 - FID, 1 - HE etc., default 0) or explicit sequence of pulses as`Sequence`

class instance.`mangetic_field`

— magnetic field along z-axis or vector of the magnetic field. Default (0, 0, 0).

```
```

```
# Time points
time_space = np.linspace(0, 2, 201) # in ms
# Number of pulses in CPMG seq (0 = FID, 1 = HE)
n = 1
# Mag. Field (Bx By Bz)
b = np.array([0, 0, 500]) # in G
l_conv = calc.compute(time_space, pulses=n, magnetic_field=b,
method='cce', quantity='coherence', as_delay=False)
```

```
```

```
%%timeit
calc.compute(time_space, pulses=n, magnetic_field=b,
method='cce', quantity='coherence', as_delay=False)
```

```
705 ms ± 9.74 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

### Generalized CCE (gCCE)

In contrast to the conventional CCE method, in generalized CCE arpproach each cluster includes the central spin explicitly.

`Simulator`

can take `pulses`

argument as an actual pulse sequence with an iterable of `Pulse`

objects.

For example:

```
p1 = pc.Pulse('x', np.pi)
p2 = pc.Pulse('y', np.pi)
seq = [p1, p2, p1, p2]
```

`seq`

will define XY-4 pulse sequence.

An integer number to define the number of pulses is also accepted as in the case of conventional CCE. If the integer is provided, the code assumes the CPMG sequence.

```
```

```
# Hahn-echo pulse sequence
pulse_sequence = [pc.Pulse('x', np.pi)]
# Calculate coherence function
l_generatilze = calc.compute(time_space, magnetic_field=b,
pulses=pulse_sequence,
method='gcce', quantity='coherence')
```

```
```

```
%%timeit
calc.compute(time_space, magnetic_field=b,
pulses=pulse_sequence,
method='gcce',
quantity='coherence')
```

```
1.92 s ± 24.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

### gCCE with random sampling of bath states

Using this approach, one may carry out generalized CCE calculations for the set of random bath states. This functionality can be turned on by by setting the keyword argument `nbstates`

to a number of bath states to sample over. Recommended number of bath states is above 100, but the convergence should be checked for each system. Note, that this computation is roughly `nbstates`

times longer than an equivalent generalized CCE calculation, as it computes everything `nbstates`

times.

For details see `help(calc.compute)`

.

```
```

```
# Number of random bath states to sample over
n_bath_states = 20
# Calculate coherence function
l_gcce = calc.compute(time_space, magnetic_field=b,
pulses=pulse_sequence,
nbstates=n_bath_states,
method='gcce', quantity='coherence', seed=seed)
```

```
```

```
%%timeit
n_bath_states = 5
calc.compute(time_space, magnetic_field=b,
pulses=pulse_sequence,
nbstates=n_bath_states,
method='gcce', quantity='coherence', seed=seed)
```

```
11.9 s ± 379 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

Take a look at the results of three different methods, and check that they produce similar coherence decay. Note that the results obtained using gCCE with bath states sampling deviates from other ones (generalized and conventional CCE), as the chosen number of states (20) is not enough to converge.

```
```

```
plt.plot(time_space, l_conv.real,
label='conventional CCE')
plt.plot(time_space, l_generatilze.real,
label='generalized CCE', ls='--')
plt.plot(time_space, l_gcce.real,
label='gCCE with bath sampling', ls='', marker='.')
plt.legend()
plt.xlabel('Time (ms)')
plt.ylabel('Coherence');
```

## Convergence parameters

Having confirmed that all methods produce the same results, we check the convergence of the conventional CCE with respect to `order`

, `r_bath`

, `r_dipole`

parameters of the `Simulator`

object.

First, define all of the parameters.

```
```

```
parameters = dict(
order=2, # CCE order
r_bath=40, # Size of the bath in A
r_dipole=8, # Cutoff of pairwise clusters in A
position=[0, 0, 0], # Position of central Spin
alpha=[0, 0, 1],
beta=[0, 1, 0],
pulses = 1, # N pulses in CPMG sequence
magnetic_field=[0,0,500]
) # Qubit levels)
time_space = np.linspace(0, 2, 201) # Time points in ms
```

We can define a little helper function to streamline the process. Note that resetting the parameters automatically recomputes the properties of the bath.

```
```

```
def runner(variable, values):
invalue = parameters[variable]
calc = pc.Simulator(spin=1, bath=atoms, **parameters)
ls = []
for v in values:
setattr(calc, variable, v)
l = calc.compute(time_space, method='cce',
quantity='coherence')
ls.append(l.real)
parameters[variable] = invalue
ls = pd.DataFrame(ls, columns=time_space, index=values).T
return ls
```

Now we can compute the coherence function at different values of the parameters:

```
```

```
orders = runner('order', [1, 2, 3, 4])
rbs = runner('r_bath', [20, 30, 40, 50, 60])
rds = runner('r_dipole', [4, 6, 8, 10])
```

We can visualize the convergence of the coherence function with respect to different parameters:

```
```

```
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
orders.plot(ax=axes[0], title='order')
rbs.plot(ax=axes[1], title='r_bath')
rds.plot(ax=axes[2], title='r_dipole')
for ax in axes:
ax.set(xlabel='Time (ms)', ylabel='Coherence')
fig.tight_layout()
```

## Bath polarization

To study different bath polarization we will modify `BathArray.state`

attribute, which contains spin states for each bath spin. For simplicity we will assume the gaussian profile of the polarization.

```
```

```
def polarize(bath, gamma=5):
# Polarizations of each bath spin
if gamma > 0:
polos = np.exp(-(bath.dist()/gamma)**2) * 0.5
else:
polos = np.zeros(bath.size)
for a, pol in zip(bath, polos):
# Skip 14N
if a.N != '13C':
continue
# Generate density matrix
dm = np.zeros((2, 2), dtype=np.complex128)
dm[0,0] = 0.5 + pol
dm[1,1] = 0.5 - pol
a.state = dm
return
# Use already optimized parameters
calc = pc.Simulator(spin=1, bath=atoms, **parameters)
# Standard deviations of the polarization gaussian profile
gammas = [0, 1, 2, 5, 10, 20, 30, 40, 60, 80]
ls = []
ts = np.linspace(0, 5, 501)
for gamma in gammas:
polarize(calc.bath, gamma=gamma)
l = calc.compute(ts)
ls.append(l.real)
df = pd.DataFrame(ls, columns=ts, index=gammas).T
```

With increased polarization in the bath the Hahn-echo signal decays significantly slower.

```
```

```
fig, ax = plt.subplots()
df.plot(cmap='magma', ax=ax)
ax.set(xlabel='Time (ms)', ylabel='Coherence');
```