Package reference

class spinsim.Device(value)

Bases: Enum

The target device that the integrator is being compiled for.

CPU = 'cpu'

Use the numba.jit() LLVM compiler to compile the integrator to run on all CPU cores, in parallel.

Note

To use this device option, the user defined field function must be numba.jit() compilable. See Supported Python features for compilable python features, and Supported Numpy features for compilable numpy features.

CPU_SINGLE = 'cpu_single'

Use the numba.jit() LLVM compiler to compile the integrator to run on a single CPU core.

Note

To use this device option, the user defined field function must be numba.jit() compilable. See Supported Python features for compilable python features, and Supported Numpy features for compilable numpy features.

CUDA = 'cuda'

Use the numba.cuda.jit() LLVM compiler to compile the integrator to run on an Nvidia cuda compatible GPU, in parallel.

Note

To use this device option, the user defined field function must be numba.cuda.jit() compilable. See Supported CUDA Python features for compilable python features.

PYTHON = 'python'

Use pure python interpreted code for the integrator, ie, don’t compile the integrator.

ROC = 'roc'

Use the numba.roc.jit() LLVM compiler to compile the integrator to run on an AMD ROCm compatible GPU, in parallel.

Warning

Work in progress, not currently functional!

class spinsim.ExponentiationMethod(value)

Bases: Enum

The implementation to use for matrix exponentiation within the integrator.

Parameters:

  • value (str) – A text label that can be used for archiving.

  • index (int) – A reference number, used when compiling the integrator, where higher level objects like enums cannot be interpreted.

ANALYTIC = 'analytic'

Analytic expression of the matrix exponential. For spin-half SpinQuantumNumber.HALF systems only. See Utilities.matrix_exponential_analytic() for more information.

LIE_TROTTER = 'lie_trotter'

Approximation using the Lie Trotter theorem, using the Pauli matrices and a single quadratic operator. See Utilities.matrix_exponential_lie_trotter() for more information.

LIE_TROTTER_8 = 'lie_trotter_8'

Approximation using the Lie Trotter theorem, using all basis elements of su(3). For spin-one SpinQuantumNumber.HALF systems only. See Utilities.matrix_exponential_lie_trotter_8() for more information.

class spinsim.IntegrationMethod(value)

Bases: Enum

Options for describing which method is used during the integration.

Parameters:

value (str) – A text label that can be used for archiving.

EULER = 'euler'

Euler integration method.

HEUN = 'heun'

Integration method from AtomicPy. Makes two Euler integration steps, one sampling the field from the start of the time step, one sampling the field from the end of the time step. The equivalent of the trapezoidal method.

MAGNUS_CF4 = 'magnus_cf4'

Commutator free, fourth order Magnus based integrator.

class spinsim.Results(time: ndarray, time_evolution: ndarray, state: ndarray, spin_calculator: callable)

Bases: object

The results of an evaluation of the integrator.

time

The times that state was evaluated at.

Type:

numpy.ndarray of numpy.float64 (time_index)

time_evolution

The evaluated time evolution operator between each time step. See Numba architecture for some information.

Type:

numpy.ndarray of numpy.float128 (time_index, y_index, x_index)

state

The evaluated quantum state of the spin system over time, written in terms of the eigenstates of the spin projection operator in the z direction.

Type:

numpy.ndarray of numpy.complex128 (time_index, magnetic_quantum_number)

spin

The expected spin projection (Bloch vector) over time. This is calculated just in time using the JITed callable spin_calculator.

Type:

numpy.ndarray of numpy.float64 (time_index, spatial_direction)

spin_calculator

Calculates the expected spin projection (Bloch vector) over time for a given time series of a quantum state. Used to calculate spin the first time it is referenced by the user.

Parameters:

  • state (numpy.ndarray of numpy.complex128 (time_index, magnetic_quantum_number)) - The quantum state of the spin system over time, written in terms of the eigenstates of the spin projection operator in the z direction.

Returns:

Type:

callable

Parameters:

  • time (numpy.ndarray of numpy.float64 (time_index)) – The times that state was evaluated at.

  • time_evolution (numpy.ndarray of numpy.float128 (time_index, y_index, x_index)) – The evaluated time evolution operator between each time step. See Numba architecture for some information.

  • state (numpy.ndarray of numpy.complex128 (time_index, magnetic_quantum_number)) – The evaluated quantum state of the spin system over time, written in terms of the eigenstates of the spin projection operator in the z direction.

  • spin_calculator (callable) –

    Calculates the expected spin projection (Bloch vector) over time for a given time series of a quantum state. Used to calculate spin the first time it is referenced by the user.

    Parameters:

    • state (numpy.ndarray of numpy.complex128 (time_index, magnetic_quantum_number)) - The quantum state of the spin system over time, written in terms of the eigenstates of the spin projection operator in the z direction.

    Returns:

class spinsim.Simulator(get_field: callable, spin_quantum_number: SpinQuantumNumber, device: Optional[Device] = None, exponentiation_method: Optional[ExponentiationMethod] = None, use_rotating_frame: bool = True, integration_method: IntegrationMethod = IntegrationMethod.MAGNUS_CF4, number_of_squares: int = 24, threads_per_block: int = 64, max_registers: Optional[int] = None, number_of_threads: Optional[int] = None)

Bases: object

spin_quantum_number

The option to select whether the simulator will integrate a spin-half SpinQuantumNumber.HALF, or spin-one SpinQuantumNumber.ONE quantum system.

Type:

SpinQuantumNumber

threads_per_block

The size of each thread block (workgroup), in terms of the number of threads (workitems) they each contain, when running on the GPU target devices Device.CUDA (Device.ROC). Defaults to 64. Modifying might be able to increase execution time for different GPU models.

Type:

int

device

The option to select which device will be targeted for integration. That is, whether the integrator is compiled for a CPU or GPU. Defaults to Device.CUDA if the system it is being run on is Nvidia Cuda compatible, and defaults to Device.CPU otherwise. See Device for all options and more details.

Type:

Device

number_of_threads

The number of CPU threads to use when running on a CPU device.

Type:

int

get_time_evolution

The internal function for evaluating the time evolution operator in parallel. Compiled for chosen device on object constrution.

Parameters:

  • sweep_parameters (numpy.ndarray of numpy.float64) - The input to the get_field() function supplied by the user. Modifies the field function so the integrator can be used for many experiments, without the need for slow recompilation. For example, if the sweep_parameters is used to define the bias field strength in get_field(), then one can run many simulations, sweeping through bias values, by calling this method multiple times, each time varying sweep_parameters.

  • time_coarse (numpy.ndarray of numpy.float64 (time_index)) - The times that state was evaluated at.

  • time_end_points (numpy.ndarray of numpy.float64 (start/end)) - The time offset that the experiment is to start at, and the time that the experiment is to finish at. Measured in s.

  • time_step_integration (float) - The integration time step. Measured in s.

  • time_step_output (float) - The sample resolution of the output timeseries for the state. Must be a whole number multiple of time_step_integration. Measured in s.

  • time_evolution_output (numpy.ndarray of numpy.float128 (time_index, y_index, x_index)) - The evaluated time evolution operator between each time step. See Numba architecture for some information.

Type:

callable

spin_calculator

Calculates the expected spin projection (Bloch vector) over time for a given time series of a quantum state. This callable is passed to the Results object returned from Simulator.evaluate(), and is executed there just in time if the spin property is needed. Compiled for chosen device on object constrution.

Parameters:

  • state (numpy.ndarray of numpy.complex128 (time_index, magnetic_quantum_number)) - The quantum state of the spin system over time, written in terms of the eigenstates of the spin projection operator in the z direction.

Returns:

Type:

callable

Parameters:

  • get_field (callable) –

    A python function that describes the field that the spin system is being put under. It must have three arguments:

    • time_sample (float) - the time to sample the field at, in units of s.

    • simulation_index (int) - a parameter that can be swept over when multiple simulations need to be run. For example, it is used to sweep over dressing frequencies during the simulations that spinsim was designed for.

    • field_sample (numpy.ndarray of numpy.float64 (spatial_index)) the returned value of the field. This is a four dimensional vector, with the first three entries being x, y, z spatial directions (to model a magnetic field, for example), and the fourth entry being the amplitude of the quadratic shift (only appearing, and required, in spin-one systems).

    Note

    This function must be compilable for the device that the integrator is being compiled for. See Device for more information and links.

  • spin_quantum_number (SpinQuantumNumber) – The option to select whether the simulator will integrate a spin-half SpinQuantumNumber.HALF, or spin-one SpinQuantumNumber.ONE quantum system.

  • device (Device) – The option to select which device will be targeted for integration. That is, whether the integrator is compiled for a CPU or GPU. Defaults to Device.CUDA if the system it is being run on is Nvidia Cuda compatible, and defaults to Device.CPU otherwise. See Device for all options and more details.

  • exponentiation_method (ExponentiationMethod) – Which method to use for matrix exponentiation in the integration algorithm. Defaults to ExponentiationMethod.LIE_TROTTER when spin_quantum_number is set to SpinQuantumNumber.ONE, and defaults to ExponentiationMethod.ANALYTIC when spin_quantum_number is set to SpinQuantumNumber.HALF. See ExponentiationMethod for more details.

  • use_rotating_frame (bool) –

    Whether or not to use the rotating frame optimisation. Defaults to True. If set to True, the integrator moves into a frame rotating in the z axis by an amount defined by the field in the z direction. This removes the (possibly large) z component of the field, which increases the accuracy of the output since the integrator will on average take smaller steps.

    Note

    The use of a rotating frame is commonly associated with the use of a rotating wave approximation, a technique used to get approximate analytic solutions of spin system dynamics. This is not done when this option is set to True - no such approximations are made, and the output state in given out of the rotating frame. One can, of course, use spinsim to integrate states in the rotating frame, using the rating wave approximation: just define get_field() with field functions that use the rotating wave approximation in the rotating frame.

  • integration_method (IntegrationMethod) – Which integration method to use in the integration. Defaults to IntegrationMethod.MAGNUS_CF4. See IntegrationMethod for more details.

  • number_of_squares (int) – The number of squares made by the matrix exponentiator, if ExponentiationMethod.LIE_TROTTER is chosen.

  • threads_per_block (int) – The size of each thread block (workgroup), in terms of the number of threads (workitems) they each contain, when running on the GPU target devices Device.CUDA (Device.ROC). Defaults to 64. Modifying might be able to increase execution time for different GPU models.

  • max_registers (int) –

    The maximum number of registers allocated per thread when using Device.CUDA as the target device, and can be modified to increase the execution speed for a specific GPU model.

    Raising this value allocates more registers (fast memory) to each thread, out of a maximum number for the whole GPU, for each specific GPU model. This means that if more registers are allocated than are available for the GPU model, the GPU must run fewer threads concurrently than it has Cuda cores, meaning some cores are inactive, and the GPU is said to have less occupancy. Lowering the value increases GPU occupancy, meaning more threads run concurrently, at the expense of fewer resgiters being avaliable to each thread, meaning slower memory must be used. Thus, there will be an optimal value of max_registers for each model of GPU running spinsim, balancing more threads vs faster running threads, and changing this value could increase performance for your GPU. See Achieved Occupancy for Nvidia’s official explanation.

  • number_of_threads (int) – The number of CPU threads to use when running on a CPU device.

compile_time_evolver(get_field: callable, spin_quantum_number: SpinQuantumNumber, device: Device, use_rotating_frame: bool = True, integration_method: IntegrationMethod = IntegrationMethod.MAGNUS_CF4, exponentiation_method: Optional[ExponentiationMethod] = None, number_of_squares: int = 24, threads_per_block: int = 64, max_registers: Optional[int] = None)

Compiles the integrator and spin calculation functions of the simulator.

Parameters:

  • get_field (callable) –

    A python function that describes the field that the spin system is being put under. It must have three arguments:

    • time_sample (float) - the time to sample the field at, in units of s.

    • sweep_parameters (numpy.ndarray of numpy.float64) - an array of parameters that can be swept over when multiple simulations need to be run. For example, it is used to sweep over dressing frequencies during the magnetometry experiments that spinsim was designed for.

    • field_sample (numpy.ndarray of numpy.float64 (spatial_index)) the returned value of the field. This is a four dimensional vector, with the first three entries being x, y, z spatial directions (to model a magnetic field, for example), and the fourth entry being the amplitude of the quadratic shift (only appearing, and required, in spin-one systems).

    Note

    This function must be compilable for the device that the integrator is being compiled for. See Device for more information and links.

  • spin_quantum_number (SpinQuantumNumber) – The option to select whether the simulator will integrate a spin-half SpinQuantumNumber.HALF, or spin-one SpinQuantumNumber.ONE quantum system.

  • device (Device) – The option to select which device will be targeted for integration. That is, whether the integrator is compiled for a CPU or GPU. Defaults to Device.CUDA if the system it is being run on is Nvidia Cuda compatible, and defaults to Device.CPU otherwise. See Device for all options and more details.

  • exponentiation_method (ExponentiationMethod) – Which method to use for matrix exponentiation in the integration algorithm. Defaults to ExponentiationMethod.LIE_TROTTER when spin_quantum_number is set to SpinQuantumNumber.ONE, and defaults to ExponentiationMethod.ANALYTIC when spin_quantum_number is set to SpinQuantumNumber.HALF. See ExponentiationMethod for more details.

  • use_rotating_frame (bool) –

    Whether or not to use the rotating frame optimisation. Defaults to True. If set to True, the integrator moves into a frame rotating in the z axis by an amount defined by the field in the z direction. This removes the (possibly large) z component of the field, which increases the accuracy of the output since the integrator will on average take smaller steps.

    Note

    The use of a rotating frame is commonly associated with the use of a rotating wave approximation, a technique used to get approximate analytic solutions of spin system dynamics. This is not done when this option is set to True - no such approximations are made, and the output state in given out of the rotating frame. One can, of course, use spinsim to integrate states in the rotating frame, using the rating wave approximation: just define get_field() with field functions that use the rotating wave approximation in the rotating frame.

  • integration_method (IntegrationMethod) – Which integration method to use in the integration. Defaults to IntegrationMethod.MAGNUS_CF4. See IntegrationMethod for more details.

  • number_of_squares (int) – The number of squares made by the matrix exponentiator, if ExponentiationMethod.LIE_TROTTER is chosen.

  • threads_per_block (int) – The size of each thread block (workgroup), in terms of the number of threads (workitems) they each contain, when running on the GPU target devices Device.CUDA (Device.ROC). Defaults to 64. Modifying might be able to increase execution time for different GPU models.

  • max_registers (int) –

    The maximum number of registers allocated per thread when using Device.CUDA as the target device, and can be modified to increase the execution speed for a specific GPU model. Defaults to 63 (optimal for GTX1070, the device used for testing. Note that one extra register per thread is always added to the number specified for control, so really this number is 64).

    Raising this value allocates more registers (fast memory) to each thread, out of a maximum number for the whole GPU, for each specific GPU model. This means that if more registers are allocated than are available for the GPU model, the GPU must run fewer threads concurrently than it has Cuda cores, meaning some cores are inactive, and the GPU is said to have less occupancy. Lowering the value increases GPU occupancy, meaning more threads run concurrently, at the expense of fewer registers being avaliable to each thread, meaning slower memory must be used. Thus, there will be an optimal value of max_registers for each model of GPU running spinsim, balancing more threads vs faster running threads, and changing this value could increase performance for your GPU. See Achieved Occupancy for Nvidia’s official explanation.

evaluate(time_start: float64, time_end: float64, time_step_integration: float64, time_step_output: float64, state_init: ndarray, sweep_parameters: ndarray = [0]) → Results

Integrates the time dependent Schroedinger equation and returns the quantum state of the spin system over time.

Parameters:

  • sweep_parameters (numpy.ndarray of numpy.float64) – The input to the get_field() function supplied by the user. Modifies the field function so the integrator can be used for many experiments, without the need for slow recompilation. For example, if the sweep_parameters is used to define the bias field strength in get_field(), then one can run many simulations, sweeping through bias values, by calling this method multiple times, each time varying sweep_parameters.

  • time_start (float) – The time offset that the experiment is to start at. Measured in s.

  • time_end (float) – The time that the experiment is to finish at. Measured in s. The duration of the experiment is time_end - time_start.

  • time_step_integration (float) – The integration time step. Measured in s.

  • time_step_output (float) – The sample resolution of the output timeseries for the state. Must be a whole number multiple of time_step_integration. Measured in s.

  • state_init (numpy.ndarray of numpy.complex128 (magnetic_quantum_number)) – The initial quantum state of the spin system, written in terms of the eigenstates of the spin projection operator in the z direction.

Returns:

results – An object containing the results of the simulation.

Return type:

Results

static get_state(state_init: ndarray, state: ndarray, time_evolution: ndarray)

Use the stepwise time evolution operators in succession to find the quantum state timeseries of the 3 level atom.

Parameters:

  • state_init (numpy.ndarray of numpy.complex128) – The state (spin wavefunction) of the system at the start of the simulation.

  • state (numpy.ndarray of numpy.complex128 (time_index, state_index)) – The state (wavefunction) of the spin system in the lab frame, for each time sampled.

  • time_evolution (numpy.ndarray of numpy.complex128 (time_index, bra_state_index, ket_state_index)) – The evaluated time evolution operator between each time step. See Numba architecture for some information.

class spinsim.SpinQuantumNumber(value)

Bases: Enum

Options for the spin quantum number of a system.

Parameters:
HALF = 0.5

For two level systems.

ONE = 1

For three level systems.

class spinsim.Utilities(spin_quantum_number: SpinQuantumNumber, device: Device, threads_per_block: int, number_of_squares: int)

Bases: object

A on object that contains definitions of all of the device functions (functions compiled for use on the target device) used in the integrator. These device functions are compiled for the chosen target device on construction of the object.

conj(z)

Conjugate of a complex number.

\[\begin{split}\begin{align*} (a + ib)^* &= a - ib\\ a, b &\in \mathbb{R} \end{align*}\end{split}\]

Parameters:

  • z (numpy.complex128) - The complex number to take the conjugate of.

Returns

  • cz (numpy.complex128) - The conjugate of z.

Type:

callable

expm1i(b)

For real input \(b\), returns \(\exp(ib) - 1\), while avoiding floating point cancellation errors.

Parameters:

  • b (numpy.float64) - The imaginary component to exponentiate.

Returns

  • em1i (numpy.complex128) - The evalauted output.

Type:

callable

cos_exp_m1(a, b)

For real input \(a\), \(b\), returns \(\cos(a)\exp(ib) - 1\), while avoiding floating point cancellation errors.

Parameters:

  • a (numpy.float64) - The real component to take the cosine of.

  • b (numpy.float64) - The imaginary component to exponentiate.

Returns

  • cem1 (numpy.complex128) - The evalauted output.

Type:

callable

cos_m1(a, b)

For real input \(a\), returns \(\cos(a) - 1\), while avoiding floating point cancellation errors.

Parameters:

  • a (numpy.float64) - The real component to take the cosine of.

Returns

  • cm1 (numpy.complex128) - The evalauted output.

Type:

callable

set_to(operator, result)

Copy the contents of one matrix into another.

\[(A)_{i, j} = (B)_{i, j}\]

Parameters:

  • operator (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix to copy from.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix to copy to.

Type:

callable

set_to_one(operator)

Make a matrix the multiplicative identity, ie, \(1\).

\[\begin{split}\begin{align*} (A)_{i, j} &= \delta_{i, j}\\ &= \begin{cases} 1,&i = j\\ 0,&i\neq j \end{cases} \end{align*}\end{split}\]

Parameters:

  • operator (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix to set to \(1\).

Type:

callable

set_to_zero(operator)

Make a matrix the zero matrix.

\[\begin{align*} (A)_{i, j} &= 0 \end{align*}\]

Parameters:

  • operator (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix to set to \(0\).

Type:

callable

matrix_multiply(left, right, result)

Multiply matrices left and right together, to be returned in result.

\[\begin{align*} (LR)_{i,k} = \sum_j (L)_{i,j} (R)_{j,k} \end{align*}\]

Parameters:

  • left (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix to left multiply by.

  • right (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix to right multiply by.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - A matrix to be filled with the result of the product.

Type:

callable

matrix_square_m1(operator, result)

For matrix \(A = 1 + a\) \(S = A^2 = 1 + s\). Here the input is the residuals \(a\), and the output is \(s\). This is a way to evaluate \(s\) without floating point cancellation error. Specifically,

\[\begin{split}\begin{align*} s &= S - 1\\ &= A^2 - 1\\ &= (2\cdot 1 + a)a \end{align*}\end{split}\]

Parameters:

  • operator (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The residual of the matrix to square.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - A matrix to be filled with the residual of the result of the product.

Type:

callable

matrix_multiply_m1(left, right, result)

For matrices \(L = 1 + l\) and \(R = 1 + r\), evaluates \(O = LR = 1 + o\). Here the inputs are the residuals \(l\) and \(r\), and the output is \(o\). This is a way to evaluate \(o\) without floating point cancellation error. Specifically,

\[\begin{split}\begin{align*} o &= O - 1\\ &= LR - 1\\ &= l + r + lr \end{align*}\end{split}\]

Parameters:

  • left (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The residual of the matrix to left multiply by.

  • right (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The residual of the matrix to right multiply by.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - A matrix to be filled with the residual of the result of the product.

Type:

callable

matrix_exponential_analytic(field_sample, result)

Calculates a \(\mathfrak{su}(2)\) matrix exponential based on its analytic form.

Warning

Only available for use with spin-half systems. Will not work with spin-one systems.

Assumes the exponent is an imaginary linear combination of \(\mathfrak{su}(2)\), being,

\[\begin{align*} A &= -i(\omega_x J_x + \omega_y J_y + \omega_z J_z), \end{align*}\]

with

\[\begin{split}\begin{align*} J_x &= \frac{1}{2}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix},& J_y &= \frac{1}{2}\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix},& J_z &= \frac{1}{2}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{align*}\end{split}\]

Then the exponential can be calculated as

\[\begin{split}\begin{align*} \exp(A) &= \exp(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z)\\ &= \begin{pmatrix} \cos(\frac{\omega_r}{2}) - i\frac{\omega_z}{\omega_r}\sin(\frac{\omega_r}{2}) & -\frac{\omega_y + i\omega_x}{\omega_r}\sin(\frac{\omega_r}{2})\\ \frac{\omega_y - i\omega_x}{\omega_r}\sin(\frac{\omega_r}{2}) & \cos(\frac{\omega_r}{2}) + i\frac{\omega_z}{\omega_r}\sin(\frac{\omega_r}{2}) \end{pmatrix} \end{align*}\end{split}\]

with \(\omega_r = \sqrt{\omega_x^2 + \omega_y^2 + \omega_z^2}\).

Parameters:

  • field_sample (numpy.ndarray of numpy.float64, (y_index, x_index)) - The values of \(\omega_x\), \(\omega_y\) and \(\omega_z\) respectively, as described above.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix which the result of the exponentiation is to be written to.

Type:

callable

matrix_exponential_lie_trotter(field_sample, result)

Calculates a matrix exponential based on the Lie Product Formula,

\[\exp(A + B) = \lim_{c \to \infty} \left(\exp\left(\frac{1}{c}A\right) \exp\left(\frac{1}{c}B\right)\right)^c.\]

For spin-half systems:

Assumes the exponent is an imaginary linear combination of a subspace of \(\mathfrak{su}(2)\), being,

\[\begin{align*} A &= -i(\omega_x J_x + \omega_y J_y + \omega_z J_z), \end{align*}\]

with

\[\begin{split}\begin{align*} J_x &= \frac{1}{2}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix},& J_y &= \frac{1}{2}\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix},& J_z &= \frac{1}{2}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{align*}\end{split}\]

Then the exponential can be approximated as, for large \(\tau\),

\[\begin{split}\begin{align*} \exp(A) =& \exp\left(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z\right)\\ =& \exp\left(2^{-\tau}\left(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z\right)\right)^{2^\tau}\\ \approx& \biggl(\exp\left(-i\frac12 2^{-\tau} \omega_z J_z\right)\exp\left(-i\left(2^{-\tau} \omega_\phi J_\phi\right)\right)\exp\left(-i\frac12 2^{-\tau} \omega_z J_z\right)\biggr)^{2^\tau}\\ =& \begin{pmatrix} \cos\left(\frac{\Phi}{2}\right)e^{-iz} & -i\sin\left(\frac{\Phi}{2}\right) e^{i\phi}\\ -i\sin\left(\frac{\Phi}{2}\right) e^{-i\phi} & \cos\left(\frac{\Phi}{2}\right)e^{iz} \end{pmatrix}^{2^\tau}\\ =& T^{2^\tau}. \end{align*}\end{split}\]

Here \(z = 2^{-\tau}\frac{\omega_z}{2}\), \(\Phi = 2^{-\tau}\sqrt{\omega_x^2 + \omega_y^2}\), and \(\phi = \mathrm{atan}2(\omega_y, \omega_x)\).

For spin-one systems

Assumes the exponent is an imaginary linear combination of a subspace of \(\mathfrak{su}(3)\), being,

\[\begin{align*} A &= -i(\omega_x J_x + \omega_y J_y + \omega_z J_z + \omega_q Q), \end{align*}\]

with

\[\begin{split}\begin{align*} J_x &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix},& J_y &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0 \end{pmatrix},\\ J_z &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix},& Q &= \frac{1}{3}\begin{pmatrix} 1 & 0 & 0 \\ 0 & -2 & 0 \\ 0 & 0 & 1 \end{pmatrix} \end{align*}\end{split}\]

Then the exponential can be approximated as, for large \(\tau\),

\[\begin{split}\begin{align*} \exp(A) =& \exp\left(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z - i\omega_q Q\right)\\ =& \exp\left(2^{-\tau}\left(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z - i\omega_q Q\right)\right)^{2^\tau}\\ \approx& \biggl(\exp\left(-i\frac12\left(2^{-\tau} \omega_z J_z + 2^{-\tau}\omega_q Q\right)\right)\nonumber\\ &\cdot\exp\left(-i\left(2^{-\tau} \omega_\phi J_\phi\right)\right)\nonumber\\ &\cdot\exp\left(-i\frac12\left(2^{-\tau} \omega_z J_z + 2^{-\tau} \omega_q Q\right)\right)\biggr)^{2^\tau}\\ =& \begin{pmatrix} \left(\cos\left(\frac{\Phi}{2}\right) e^{-iz}e^{-iq}\right)^2 & \frac{-i}{\sqrt{2}} \sin(\Phi)e^{iq}e^{-iz}e^{-i\phi} & -\left(\sin\left(\frac{\Phi}{2}\right)e^{iq}e^{-i\phi}\right)^2\\ \frac{-i}{\sqrt{2}} \sin(\Phi)e^{iq}e^{-iz}e^{i\phi} & \cos(\Phi)e^{i4q} & \frac{-i}{\sqrt{2}} \sin(\Phi)e^{iq}e^{iz}e^{-i\phi}\\ -\left(\sin\left(\frac{\Phi}{2}\right)e^{-iq}e^{i\phi}\right)^2 & \frac{-i}{\sqrt{2}} \sin(\Phi)e^{iq}e^{iz}e^{i\phi} & \left(\cos\left(\frac{\Phi}{2}\right) e^{iz}e^{-iq}\right)^2 \end{pmatrix}^{2^\tau}.\\ \end{align*}\end{split}\]

Here \(z = 2^{-\tau}\frac{\omega_z}{2}\), \(q = 2^{-\tau}\frac{\omega_q}{6}\), \(\Phi = 2^{-\tau}\sqrt{\omega_x^2 + \omega_y^2}\), and \(\phi = \mathrm{atan}2(\omega_y, \omega_x)\). Once \(T\) is calculated, it is then recursively squared \(\tau\) times to obtain \(\exp(A)\).

Parameters:

  • field_sample (numpy.ndarray of numpy.float64, (y_index, x_index)) - The values of x, y and z (and q for spin-one) respectively, as described above.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix which the result of the exponentiation is to be written to.

  • number_of_squares (int) - The number of squares to make to the approximate matrix (\(\tau\) above).

Type:

callable

matrix_exponential_lie_trotter_8(field_sample, result)

Calculates a matrix exponential based on the Lie Product Formula,

\[\exp(A + B) = \lim_{c \to \infty} \left(\exp\left(\frac{1}{c}A\right) \exp\left(\frac{1}{c}B\right)\right)^c.\]

Warning

Only available for use with spin-one systems. Will not work with spin-half systems.

Assumes the exponent is an imaginary linear combination elements of \(\mathfrak{su}(3)\), being,

\[\begin{align*} A &= -i(\omega_x J_x + \omega_y J_y + \omega_z J_z + \omega_q Q + \omega_{u1} U_1 + \omega_{u2} U_2 + \omega_{v1} V_1 + \omega_{v2} V_2), \end{align*}\]

with

\[\begin{split}\begin{align*} J_x &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix},& J_y &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & -i & 0 \\ i & 0 & -i \\ 0 & i & 0 \end{pmatrix},\\ J_z &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix},& Q &= \frac{1}{3}\begin{pmatrix} 1 & 0 & 0 \\ 0 & -2 & 0 \\ 0 & 0 & 1 \end{pmatrix},\\ U_1 &= \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{pmatrix},& U_2 &= \begin{pmatrix} 0 & 0 & -i \\ 0 & 0 & 0 \\ i & 0 & 0 \end{pmatrix},\\ V_1 &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & -1 \\ 0 & -1 & 0 \end{pmatrix},& V_2 &= \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & -i & 0 \\ i & 0 & i \\ 0 & -i & 0 \end{pmatrix}.\\ \end{align*}\end{split}\]

Then the exponential can be approximated as, for large \(\tau\),

\[\begin{split}\begin{align*} \exp(A) =& \exp\biggl(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z - i\omega_q Q\\ &- i\omega_{u1} U_1 - i\omega_{u2} U_2 - i\omega_{v1} V_1 - i\omega_{v2} V_2\biggr)\\ & \exp\biggl(2^{-\tau}\biggl(-i\omega_x J_x - i\omega_y J_y - i\omega_z J_z - i\omega_q Q\\ &- i\omega_{u1} U_1 - i\omega_{u2} U_2 - i\omega_{v1} V_1 - i\omega_{v2} V_2\biggr)\biggr)^{2^\tau}\\ \approx& \biggl(\exp\left(-i2^{-\tau} \omega_\phi J_\phi\right)\exp\left(-i2^{-\tau} \omega_{u\phi} U_{u\phi}\right)\\ &\cdot\exp\left(-i2^{-\tau} \omega_{v\phi} V_{v\phi}\right)\exp\left(-i2^{-\tau} \omega_z J_z -i2^{-\tau} \omega_q Q \right)\biggr)^{2^\tau}\\ =& \biggl(\begin{pmatrix} \cos^2\left(\frac{\Phi}{2}\right) & \frac{-i}{\sqrt{2}} \sin(\Phi)e^{-i\phi} & -\left(\sin\left(\frac{\Phi}{2}\right)e^{-i\phi}\right)^2\\ \frac{-i}{\sqrt{2}} \sin(\Phi)e^{i\phi} & \cos\left(\Phi\right) & \frac{-i}{\sqrt{2}} \sin(\Phi)e^{-i\phi}\\ -\left(\sin\left(\frac{\Phi}{2}\right)e^{i\phi}\right)^2 & \frac{-i}{\sqrt{2}} \sin(\Phi)e^{i\phi} & \cos^2\left(\frac{\Phi}{2}\right) \end{pmatrix}\\ &\cdot \begin{pmatrix} \cos\left(\Phi_u\right) & 0 & -i \sin\left(\Phi_u\right)e^{-i\phi_u}\\ 0 & 1 & 0\\ -i \sin\left(\Phi_u\right)e^{i\phi_u} & 0 & \cos\left(\Phi_u\right) \end{pmatrix}\\ &\cdot \begin{pmatrix} \cos^2\left(\frac{\Phi_v}{2}\right) & \frac{-i}{\sqrt{2}} \sin(\Phi_v)e^{-i\phi_v} & \left(\sin\left(\frac{\Phi_v}{2}\right)e^{-i\phi_v}\right)^2\\ \frac{-i}{\sqrt{2}} \sin(\Phi_v)e^{i\phi_v} & \cos\left(\Phi_v\right) & \frac{i}{\sqrt{2}} \sin(\Phi_v)e^{-i\phi_v}\\ \left(\sin\left(\frac{\Phi_v}{2}\right)e^{i\phi_v}\right)^2 & \frac{i}{\sqrt{2}} \sin(\Phi_v)e^{i\phi_v} & \cos^2\left(\frac{\Phi_v}{2}\right) \end{pmatrix}\\ &\cdot \begin{pmatrix} e^{-iz - iq} & 0 & 0\\ 0 & e^{i2q} & 0\\ 0 & 0 & e^{iz - iq} \end{pmatrix}\biggr)^{2^\tau}\\ =& T^{2^\tau}. \end{align*}\end{split}\]

Here \(z = 2^{-\tau}\frac{\omega_z}{2}\), \(q = 2^{-\tau}\frac{\omega_q}{6}\), \(\Phi = 2^{-\tau}\sqrt{\omega_x^2 + \omega_y^2}\), \(\phi = \mathrm{atan}2(\omega_y, \omega_x)\), \(\Phi_u = 2^{-\tau}\sqrt{\omega_{u1}^2 + \omega_{u2}^2}\), \(\phi_u = \mathrm{atan}2(\omega_{u1}, \omega_{u2})\), \(\Phi_v = 2^{-\tau}\sqrt{\omega_{v1}^2 + \omega_{v2}^2}\), and \(\phi_v = \mathrm{atan}2(\omega_{v1}, \omega_{v2})\). Once \(T\) is calculated, it is then recursively squared \(\tau\) times to obtain \(\exp(A)\).

Parameters:

  • field_sample (numpy.ndarray of numpy.float64, (y_index, x_index)) - The values of x, y, z, q, u1, u2, v1 and v2 respectively, as described above.

  • result (numpy.ndarray of numpy.complex128, (y_index, x_index)) - The matrix which the result of the exponentiation is to be written to.

  • number_of_squares (int) - The number of squares to make to the approximate matrix (\(\tau\) above).

Type:

callable

Parameters:

  • spin_quantum_number (SpinQuantumNumber) – The option to select whether the simulator will integrate a spin-half SpinQuantumNumber.HALF, or spin-one SpinQuantumNumber.ONE quantum system.

  • device (Device) – The option to select which device will be targeted for integration. That is, whether the integrator is compiled for a CPU or GPU. Defaults to Device.CUDA if the system it is being run on is Nvidia Cuda compatible, and defaults to Device.CPU otherwise. See Device for all options and more details.

  • threads_per_block (int) – The size of each thread block (workgroup), in terms of the number of threads (workitems) they each contain, when running on the GPU target devices Device.CUDA (Device.ROC). Defaults to 64. Modifying might be able to increase execution time for different GPU models.

spinsim.generate_interpolation_sampler(time: ndarray, amplitude: ndarray, order=3, zero_boundary=True, device: Optional[Device] = None)

Uses a user provided time series to create a field sampling function that is compatible with the Simulator object. Samples in the time series do not have to be equally spaced, and as such the user much specify a numpy.ndarray for both time and amplitude.

Note

These interpolators can be set up in a variety of ways, such pulse trains where the output modulates a sinusoid. See the examples to learn more.

The sampler can be defined to use one of three kinds of interpolation, defined by the parameter order. Below we use \(t_j\) and \(b_j\) for the time and amplitude samples respectively, \(t\) for the sampling time, and \(\omega(t)\) for the field value at the sampling time.

Use order = 0 to use the Zero Order Hold (ZOH) protocol:

\[\omega_a(t) += b_j \quad \text{for } t_j < t < t_{j + 1}\]

This method can be used for simulating hard pulse trains.

Use order = 1 to use the linear interpolation (LERP) protocol:

\[\begin{split}\tau =& \frac{t - t_j}{t_{j + 1} - t_j},\\ \omega_a(t) +=& (1 - \tau) b_j + \tau b_{j + 1} \quad \text{for } t_j < t < t_{j + 1}\end{split}\]

This method can be used for simulating linear ramps.

Use order = 3 to use the cubic Hermite spline using the finite difference gradient:

\[\begin{split}\tau =& \frac{t - t_j}{t_{j + 1} - t_j},\\ m_j =& \frac12 \left( \frac{b_{j + 1} - b_j}{t_{j + 1} - t_j} + \frac{b_j - b_{j - 1}}{t_j - t_{j - 1}} \right) \\ \omega_a(t) +=& (2\tau^3 - 3\tau^2 + 1) b_j + (-2\tau^3 + 3\tau^2) b_{j + 1} + (\tau^3 - 2\tau^2 + \tau) (t_{j + 1} - t_j) m_j \\ &+ (\tau^3 - \tau^2) (t_{j + 1} - t_j) m_{j + 1} \quad \text{for } t_j < t < t_{j + 1}\end{split}\]

This method is best for simulating based off of recorded traces or smooth functions that cannot be constructed with basic functions in math or numpy.

For any value of order, if \(t < t_0\) or \(t > t_\text{end}\) then there are two options. If zero_boundary is set to True (default), then \(\omega_a(t) = 0\). Otherwise, if \(t < t_0\), then \(\omega_a(t) = b_0\), and if \(t > t_\text{end}\), then \(\omega_a(t) = b_\text{end}\).

Parameters:

  • time (numpy.ndarray of numpy.float64 (time_index)) – The times that correspond to the sampled amplitudes in amplitude. Note that times do not have to be equally spaced.

  • amplitude (numpy.ndarray of numpy.float64 (time_index)) – The sampled amplitudes to be interpolated between.

  • order (int) – The interpolation method to be used by the interpolator. Defaults to 3 (for cubic interpolation). For more information, see above.

  • zero_boundary (bool) – Whether or not the boundaries (ie times before the first and after the final sample) should be zero padded. Defaults to True (they will be zero padded). For more information, see above.

  • device (Device) – The device that the interpolator should be compiled to. If None (the default), it will compile to CUDA if a CUDA compatible device is found, and multithread CPU otherwise. To be used with a Simulator instance, you must select the same device as the one selected by that. If the interpolator used in a context outside of spinsim (plotting the interpolator output, for example), then you must select Device.CPU or Device.CPU_SINGLE.

Returns:

sampler – An interpolating function based off the time series that can be used as or within a get_field() function for a Simulator instance.

Parameters:

  • time_sample (float) - The time at which to sample the timeseries.

  • user_parameters (numpy.ndarray of numpy.float64, (parameter_index)) - Unused here, but a required format for the Simulator instance to read.

  • field_sample (numpy.ndarray of numpy.float64, (y_index, x_index)) - The array of field values that this function adds to.

Returns:

  • field_sample (float) - The interpolated value of the amplitude.

Return type:

callable