`core.py`: Core features¶

The core features of this library

simba.core.is_transfer_matrix_physically_realisable(expr: sympy.matrices.dense.MutableDenseMatrix, d_matrix: Optional[sympy.matrices.dense.MutableDenseMatrix] = None) → bool[source]¶

Check if transfer matrix $\mathbf{G}(s)$ given by expr and direct-feed matrix d_matrix is possible to physically realise by the conditions given in [transfer-function], i.e. that it obeys,

$\mathbf{G}^\sim(s) J \mathbf{G}(s) = J,\ \ D J D^\dagger = J,$

where $\mathbf{G}^\sim(s) \equiv \mathbf{G}^\dag(-s^*)$ .

Parameters

expr – a sympy Matrix representing the transfer Matrix $\mathbf{G}(s)$
d_matrix – a sympy Matrix representing the direct-feed Matrix, defaulting to the identity matrix

simba.core.transfer_func_coeffs_to_state_space(numer: List, denom: List) → simba.core.StateSpace[source]¶: See StateSpace.from_transfer_function_coeffs.

simba.core.transfer_function_to_state_space(expr: sympy.core.expr.Expr) → simba.core.StateSpace[source]¶: See StateSpace.from_transfer_function.

simba.core.tf2ss(expr: sympy.core.expr.Expr) → simba.core.StateSpace[source]¶: See transfer_function_to_state_space

simba.core.transfer_function_to_realisable_state_space(expr: sympy.core.expr.Expr) → simba.core.StateSpace[source]¶: Convert given transfer function to physically realisable state space if possible.

simba.core.tf2rss(expr: sympy.core.expr.Expr) → simba.core.StateSpace[source]¶: See transfer_function_to_realisable_state_space

simba.core.tf2network(expr: sympy.core.expr.Expr) → simba.core.SplitNetwork[source]¶: Calculate split network directly from transfer function

class simba.core.StateSpace(a, b, c, d, *, paired_operator_form=False)[source]¶

Represents a dynamical quantum state-space which describes the time-domain evolution of a system.

$\dot{x} &= a x + b u, \\ y &= c x + d u.$

where the state vectors are bounded linear operators usually in doubled-up form, (use reorder_to_paired_form to convert to paired operator form)

$x \in \mathbb{L}^{n\times 1},\ u \in \mathbb{L}^{m\times 1},\ y \in \mathbb{L}^{l\times 1},$

and, the “system matrices” are

$a \in \mathbb{C}^{n\times n},\ b \in \mathbb{C}^{n\times m},\ c \in \mathbb{C}^{l\times n},\ d \in \mathbb{C}^{l\times m}.$

If the dimensions do not match up, raises DimensionError.

For SISO quantum systems $m = l = 2$ (usually the operator and its conjugate operator, e.g. $u = (\hat{u}, \hat{u}^\dagger)^T$ ).

Attributes:

a: Internal dynamics
b: Input coupling
c: Internal coupling to output
d: Direct-feed (input to output)
paired_operator_form: Whether or not system is in paired operator form.

Note that the matrices are stored using sympy.ImmutableMatrix, and so are immutable. New StateSpaces should be created from modified matrices instead.

classmethod from_transfer_function_coeffs(numer: List, denom: List) → simba.core.StateSpace[source]¶

Return the SISO controllable canonical form state space for the given list of numerators and denominators of a pole-zero form transfer function, given in order of ascending powers, assuming complex ladder operators $(a, a^\dagger)$ are used.

The coefficients are defined via the transfer function between the input $u(s)$ and the output $y(s)$ , where $s$ is the complex Laplace frequency, 2

$\frac{y(s)}{u(s)} = \frac{b_0 s^n + b_1 s^{n-1} + \dots + b_{n-1} s + b_n} {s^n + a_1 s^{n-1} + \dots + a_{n-1} s + a_n}.$

Note that we assume the coefficients are normalized with respect to the highest order term in the denominator.

Raises CoefficientError if lengths of coefficient lists are wrong.

Reference: https://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf

TODO: implement for MIMO systems

Parameters

numer – The numerator coefficients: $[b_n, \dots, b_0]$
denom – The denominator coefficients: $[a_n, \dots, a_1]$

Returns

StateSpace for the given system

classmethod from_transfer_function(expr: sympy.core.expr.Expr) → simba.core.StateSpace[source]¶: Call from_transfer_function_coeffs passing the expression to transfer_function_to_coeffs.

to_transfer_function() → sympy.core.expr.Expr[source]¶: Calculate transfer function matrix for the system using the convention given by 2.

extended_to_quantum() → simba.core.StateSpace[source]¶: Extend SISO state-space to quantum MIMO state space in doubled-up ordering (see 1). Returns extended StateSpace. Does not modify original.

reorder_to_paired_form() → simba.core.StateSpace[source]¶

Return a new StateSpace with the system matrices reordered so that the state vectors, inputs, and outputs are converted from doubled-up form,

$(a_1, a_2, \dots, a_n; a_1^\dagger, a_2^\dagger, \dots, a_n^\dagger)^T,$

to paired operator form,

$(a_1, a_1^\dagger; a_2, a_2^\dagger; \dots; a_n, a_n^\dagger)^T,$

Does nothing if self.paired_operator_form is True

Returns: StateSpace in paired up form.

find_transformation_to_physically_realisable() → sympy.matrices.dense.MutableDenseMatrix[source]¶

Return the $T$ matrix that transforms the state space into a physically realisable one.

Raise StateSpaceError if system is not possible to physically realise.

Raise ResultError if there was some other unexpected error during finding T.

to_physically_realisable() → simba.core.StateSpace[source]¶

Return copy of state space transformed to a physically realisable state-space, or just return self if already physically realisable.

Transforms to paired operator form first if needed.

Raise DimensionError if system does not have even number of degrees of freedom.

property is_physically_realisable: bool¶: Test physical realisability conditions using the paired operator form of the state-space.

$AJ + JA^\dagger + BJB^\dagger &= 0, \\ JC^\dagger + BJD^\dagger &= 0.$

raise_error_if_not_possible_to_realise() → None[source]¶

Raises StateSpaceError if system cannot be physically realised according to the conditions given in is_transfer_matrix_physically_realisable.

Note this does not imply that the system matrices $(A, B, C, D)$ are physically realisable, just that they can be transformed to a physically realisable form.

to_skr() → simba.core.SKR[source]¶: Convert state space to SLH form as discussed in [synthesis], specifically returning the matrices $(S, K, R)$ . Assume physically realisable but won’t error if it’s not.

to_slh(symbol='a') → simba.core.SLH[source]¶: Create StateSpace.to_skr returning a SLH object using given symbol name symbol, defaulting to ‘a’.

property num_degrees_of_freedom: int¶: Returns num degrees of freedom if classical, or 2 * num degrees of freedom if quantum.

property num_inputs: int¶: Returns num inputs if classical, or 2 * num inputs if quantum.

property num_outputs: int¶: Returns num outputs if classical, or 2 * num outputs if quantum.

pprint() → None[source]¶: Pretty print the system matrices for debug or interactive programming purposes.

__iter__() → iter[source]¶: Returns iterator holding tuple with the four system matrices. Use for unpacking.

__eq__(other) → bool[source]¶: Equality for state spaces means that all the ABCD matrices are equal and both are or aren’t quantum.

__str__() → str[source]¶: Prettify the equation.

simba.core.j_matrix(num_dof: int) → sympy.matrices.dense.MutableDenseMatrix[source]¶

Return quantum $J$ matrix for a paired operator form StateSpace with given num_dof,

$J = \text{diag}(1, -1;\dots;1, -1) \in \mathbb{R}^{2n\times 2n}.$

Raises ValueError if num_dof is not even.

simba.core.transfer_function_to_coeffs(expr: sympy.core.expr.Expr, flip_s=True) → simba.core.Coefficients[source]¶

Extract transfer function coefficients from the given expression which is a function of frequency $s$ and a ratio of two polynomials.

Returns a namedtuple containing two lists of length n and n-1 for the numerator and denominator respectively in order of ascending powers, where n is the order of expr in $s$ .

The numerator list will be padded if the denominator is higher order than the numerator.

A NotImplementedError will be raised if the denominator is lower order than the numerator.

The coefficients are normalised such that the coefficient of the highest order term in the denominator is one.

Parameters

expr – ratio of two polynomials in $s$
flip_s – if True, set s to -s, that is use the Laplace convention defined in 2

Returns

Coefficients instance

class simba.core.Coefficients(numer: List, denom: List)[source]¶

Represents the transfer function coefficients as returned by transfer_function_to_coeffs

numer: List¶: Alias for field number 0

denom: List¶: Alias for field number 1

simba.core.concat(a: simba.core.SLH, b: simba.core.SLH) → simba.core.SLH[source]¶

Concatenate two SLH systems using the concatenation product. [synthesis]

Let $G_1 = (S_1, K_1 x_{1,0}, \frac{1}{2} x_{1,0}^\dag R_2 x_{1,0}$ ) and $G_2 = (S_2, K_2 x_{2,0}, \frac{1}{2} x_{2,0}^\dag R_2 x_{2,0})$ , where $x_{k,0} \equiv x_k(t=0)$ . Then the concatenation product is defined as,

$G_1 \boxplus G_2 = \left(S_{1\boxplus2}, (K_1x_{1,0},K_2x_{2,0})^T, \frac{1}{2} x_{1,0}^\dag R_2 x_{1,0}+\frac{1}{2} x_{2,0}^\dag R_2 x_{2,0}\right),$

where,

$S_{1\boxplus2} = \begin{bmatrix}S_1 & 0 \\ 0 & S_2\end{bmatrix}.$

For simplicity we assume the system’s share no degrees of freedom, i.e. $x_{1,0} \neq x_{2,0}$ .

Parameters

a – SLH object representing generalised open oscillator $G_1$
b – SLH object representing generalised open oscillator $G_2$

Returns

SLH object representing concatenation of both

simba.core.series(g_to, g_from)[source]¶

Not yet implemented

Series product representing the feeding of the output of g_from into the input of g_to. The arguments are in this order to match the notation below. The generalised open oscillators $G_1$ and $G_2$ are as defined in concat.

The series product is then defined as,

$G_2 \triangleleft G_1 = \bigg(&S_2S_1, K_2x_{2,0} + S_2 K_1x_{1,0}, \\ &\frac{1}{2} x_{1,0}^\dag R_1 x_{1,0} + \frac{1}{2} x_{2,0}^\dag R_2 x_{2,0} + \frac{1}{2i} x_{2,0}^\dag (K_2^\dag S_2 K_1 - K_2^T S_2^* K_1^*) x_{1,0}\bigg).$

Parameters

g_from – SLH object representing generalised open oscillator $G_1$
g_to – SLH object representing generalised open oscillator $G_2$

Returns

SLH object representing concatenation of both

class simba.core.SKR(s: sympy.matrices.dense.MutableDenseMatrix, k: sympy.matrices.dense.MutableDenseMatrix, r: sympy.matrices.dense.MutableDenseMatrix)[source]¶

Represents SKR matrices as returned by StateSpace.to_skr

s: sympy.matrices.dense.MutableDenseMatrix¶: Alias for field number 0

k: sympy.matrices.dense.MutableDenseMatrix¶: Alias for field number 1

r: sympy.matrices.dense.MutableDenseMatrix¶: Alias for field number 2

class simba.core.SLH(s, k, r, x0=None)[source]¶

Represents a generalised open oscillator in the SLH formalism. [synthesis]

Attributes: (all stored as sympy ImmutableMatrix objects)

s: scattering matrix
k: linear coupling matrix
r: internal system Hamiltonian matrix
x0: vector of system state Symbols

split() → simba.core.SplitNetwork[source]¶: Returns split_system(self).

property interaction_hamiltonian¶: Returns the interaction Hamiltonian for the system via interaction_hamiltonian_from_linear_coupling_operator(self.k * self.x0)

simba.core.interaction_hamiltonian_from_linear_coupling_operator(l_operator: sympy.matrices.dense.MutableDenseMatrix) → sympy.core.expr.Expr[source]¶

Calculate the idealised interaction hamiltonian for the given linear coupling operator.

$H_\text{int} = i[L^\dagger\ -L^\dagger] u,$

where $u = (u_1, u_1^\dagger; \dots; u_m, u_m^\dagger)^T$ and $L \in \mathbb{L}^{m\times1}$

Raises DimensionError if not l_operator not a column vector.

simba.core.linear_coupling_operator_from_k_matrix(k_matrix: sympy.matrices.dense.MutableDenseMatrix, symbol: str = 'a') → sympy.matrices.dense.MutableDenseMatrix[source]¶

Calculate symbolic linear coupling operator from K matrix assuming paired operator form.

$L = K x_0,$

where $x_0 = (a_1, a_1^\dagger; \dots; a_n, a_n^\dagger)^T$ .

symbol is the symbol name to use for the state variables.

Raises DimensionError if not an even number of dimensions.

simba.core.hamiltonian_from_r_matrix(r_matrix: sympy.matrices.dense.MutableDenseMatrix, symbol: str = 'a') → sympy.core.expr.Expr[source]¶

Calculate symbolic internal Hamiltonian from R matrix assuming paired operator form.

$H = x_0^\dagger R x_0,$

where $x_0 = (a_1, a_1^\dagger; \dots; a_n, a_n^\dagger)$ and $R \in \mathbb{R}^{2n\times2n}$ .

symbol is the symbol name to use for the state variables.

Raises DimensionError if not an even number of dimensions or if not square.

simba.core.make_complex_ladder_state(num_dofs: int, symbol: str = 'a') → sympy.matrices.dense.MutableDenseMatrix[source]¶

Return matrix of complex ladder operators with 2 * num_dofs elements.

Use symbol keyword arg to set alternative symbol for variables instead of using a

For example, for num_dofs == 2, result is $(a_1, a_1^\dagger; a_2, a_2^\dagger)^T$ .

If num_dofs == 1 then just returns $(a, a^\dagger)^T$

class simba.core.States(states: sympy.matrices.dense.MutableDenseMatrix)[source]¶

Represents a collection of state variables that can be queried to extract different variables.

get_symbol(name: str) → Optional[sympy.core.symbol.Symbol][source]¶: Get symbol from self.states with given name. To get a conjugated variable, e.g. $a_1^\dagger$ , use conjugate(a_1). To get a primed variable write a'_1 rather than a_1'.

get_symbols(names: List[str]) → List[sympy.core.symbol.Symbol][source]¶: Same as get_symbol but takes and returns a list.

index_of(symbol: sympy.core.symbol.Symbol) → Optional[int][source]¶: Return index of symbol in self.states, or None if not found.

class simba.core.SplitNetwork(gs: List[simba.core.SLH], h_d: sympy.matrices.dense.MutableDenseMatrix)[source]¶

Represents a tuple of ([G_1, …, G_n], H^d), as returned by split_system.

Each $G_i$ is assumed to be series fed into $G_{i+1}$ , however if $L = 0$ for $G_i$ it can be considered as not being series connected as the external field will only couple to the auxiliary degree of freedom which is adiabatically eliminated.

__iter__()[source]¶

__str__()[source]¶: Return str(self).

property states: sympy.matrices.dense.MutableDenseMatrix¶: Return the state symbols for the main and auxiliary mode for the entire network.

property input_output_symbols: sympy.matrices.dense.MutableDenseMatrix¶: Return the symbols for the input and output fields.

class InteractionHamiltonian(h: sympy.matrices.dense.MutableDenseMatrix)[source]¶

Represents the interaction Hamiltonian as a Matrix in order

$(a_1, a_1^\dagger, \dots, a_n, a_n^\dagger; a_1', a_1'^\dagger,\dots, a_n', a_n'^\dagger),$

where $a_i$ is the annihilation operator for the main cavity mode of the i-th system and $a_i'$ is for the corresponding auxiliary mode.

property states: sympy.matrices.dense.MutableDenseMatrix¶: Return the state symbols for the main and auxiliary modes for the network.

property dynamical_matrix: sympy.matrices.dense.MutableDenseMatrix¶: Compute dynamical matrix as shown in notes/eqns-of-motion-from-hamiltonian-matrix.pdf.

property equations_of_motion: sympy.matrices.dense.MutableDenseMatrix¶

Compute the Heisenberg equations of motion of the interacting terms (not including the Langevin equations) , using the Bosonic commutation relations

$[a_i, a_j] = 0, \quad [a_i, a_j^\dagger] = \delta_{i, j},$

and

$[a_i', a_j'] = 0, \quad [a_i', a_j'^\dagger] = \delta_{i, j}$

Returns eqns in order $(a_1, a_1^\dagger, \dots, a_n, a_n^\dagger; a_1', a_1'^\dagger,\dots, a_n', a_n'^\dagger)$ .

:returns column vector of equations of motion

property interaction_hamiltonian: simba.core.SplitNetwork.InteractionHamiltonian¶: Compute the interaction Hamiltonian between internal degrees of freedom in the network, not including the external continuum fields.

property input_output_eqns: sympy.matrices.dense.MutableDenseMatrix¶: Calculate the input-output equations for the auxiliary fields as column vector.

property state_vector: sympy.matrices.dense.MutableDenseMatrix¶

Returns column vector of states as follows:

$(a_1, a_1^\dagger, \dots, a_n, a_n^\dagger; a_1', a_1'^\dagger,\dots, a_n', a_n'^\dagger; a_\text{in 1}, a_\text{in 1}^\dagger, a_\text{in m},a_\text{in m}^\dagger; a_\text{out 1}, a_\text{out 1}^\dagger, a_\text{out m}, a_\text{out m}^\dagger)^T$

for system with n internal degrees of freedom and m inputs and outputs.

class DynamicalMatrix(matrix: sympy.matrices.dense.MutableDenseMatrix, states: Union[sympy.matrices.dense.MutableDenseMatrix, simba.core.States])[source]¶

Represents the dynamical matrix as returned by SplitNetwork.dynamical_matrix.

Used to calculate the transfer functions between any two quantities in the system.

class TransferMatrix(matrix: sympy.matrices.dense.MutableDenseMatrix, states: Union[sympy.matrices.dense.MutableDenseMatrix, simba.core.States])[source]¶

Represents the results of DynamicalMatrix.transfer_matrix. Used to extract the transfer functions.

closed_loop(excitation: Union[sympy.core.symbol.Symbol, str], variable: Union[sympy.core.symbol.Symbol, str]) → sympy.core.expr.Expr[source]¶

Get the closed loop transfer function from the excitation to the variable.

The variables can be given as strings (which are passed to States.get_symbol) or as sympy Symbols.

E.g. closed_loop("ain", "aout") is the closed-loop transfer function from ain to aout.

closed_loop_gain(excitation: Union[sympy.core.symbol.Symbol, str])[source]¶

Calculate the closed loop gain from the excitation to the corresponding variable.

I.e. if excitation is $\hat{a}_\text{exc}$ then returns the transfer function from $\hat{a}_\text{exc}$ to $\hat{a}$ .

open_loop(excitation: Union[sympy.core.symbol.Symbol, str], variable: Union[sympy.core.symbol.Symbol, str]) → sympy.core.expr.Expr[source]¶: Calculate the open loop transfer function from the excitation to the variable (the closed loop transfer function divided by the closed loop gain).

property transfer_matrix: simba.core.SplitNetwork.DynamicalMatrix.TransferMatrix¶: Calculate the transfer matrix for the system by adding an excitation to each mode.

property eqns: sympy.matrices.dense.MutableDenseMatrix¶: Calculate matrix of RHS of frequency-domain equations $v = Mv$ .

property dynamical_matrix: simba.core.SplitNetwork.DynamicalMatrix¶

Calculate the full frequency-domain dynamical matrix $M$ for the system, include input and output equations, such that $v = M v$ where $v$ is the full frequency-domain state vector in the order returned by state_vector.

TODO: use a sparse matrix

property transfer_matrix: simba.core.SplitNetwork.DynamicalMatrix.TransferMatrix¶: Shortcut for self.dynamical_matrix.transfer_matrix

property tfm: simba.core.SplitNetwork.DynamicalMatrix.TransferMatrix¶: Shortcut for self.dynamical_matrix.transfer_matrix

property aux_coupling_constants: List[sympy.core.symbol.Symbol]¶

Return a list of the Sympy symbols use for the coupling constants for the auxiliary modes.

Examples:

>>> from sympy import symbols
>>> s = symbols('s')
>>> gamma_f, omega_s = symbols('gamma_f omega_s', real=True, positive=True)
>>> network = tf2network((s**2 + s * gamma_f + omega_s**2) / (s**2 - s * gamma_f + omega_s**2))
>>> gamma_1, gamma_2 = network.aux_coupling_constants

simba.core.split_system(open_osc: simba.core.SLH) → simba.core.SplitNetwork[source]¶

Split n degree of freedom open oscillator into n one degree of freedom open oscillators and a direct interaction Hamiltonian matrix.

Parameters: open_osc – the n degree of freedom open oscillator to split
Returns: a SplitNetwork which represents a tuple of ([G_1, …, G_n], H^d)

Footnotes

1

Quantum in this case meaning that the following transformation is applied:

$\dot{x} = a x + b u,\quad y = c x + d u$

to,

$\begin{bmatrix}\dot{x} \\ \dot{x}^\#\end{bmatrix} &= \begin{bmatrix}a & 0 \\ 0 & a^\#\end{bmatrix} \begin{bmatrix}x \\ x^\#\end{bmatrix} + \begin{bmatrix}b & 0 \\ 0 & b^\#\end{bmatrix} \begin{bmatrix}u \\ u^\#\end{bmatrix}, \\ \begin{bmatrix}y \\ y^\#\end{bmatrix} &= \begin{bmatrix}c & 0 \\ 0 & c^\#\end{bmatrix} \begin{bmatrix}x \\ x^\#\end{bmatrix} + \begin{bmatrix}d & 0 \\ 0 & d^\#\end{bmatrix} \begin{bmatrix}u \\ u^\#\end{bmatrix},$

where for a matrix $m$ , the notation $m^\#$ means “take the adjoint of each element”. Effectively each vector is in doubled-up form, as discussed in [squeezing-components].

2(1,2,3)