API Reference

This page documents the public modules and objects exposed by the package.

Initial States

class gaussian_systems.initial_state.GaussianCVState(mean_vector, covariance_matrix)

Bases: object

Container for an n-mode Gaussian continuous-variable state.

This class represents a Gaussian CV state by its mean vector and covariance matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The state is validated on construction and internal copies of the supplied arrays are stored. The mean_vector and covariance_matrix properties also return copies, so external code cannot mutate the internal state accidentally.

Parameters:

  • mean_vector (numpy.ndarray) – The Gaussian mean vector of shape (2n,).

  • covariance_matrix (numpy.ndarray) – The Gaussian covariance matrix of shape (2n, 2n).

Raises:

  • TypeError – If the inputs are not valid NumPy arrays.

  • ValueError – If the mean vector and covariance matrix are not valid or are dimensionally inconsistent.

Notes

This class stores only first and second moments. It therefore represents Gaussian states completely, but does not encode non-Gaussian information.

copy_state()

Create a deep copy of the Gaussian state.

This method returns a new GaussianCVState instance with copies of the mean vector and covariance matrix. The returned state is completely independent of the original.

Returns:

A deep copy of the current state.

Return type:

GaussianCVState

Notes

Modifications to the returned state will not affect the original state.

property covariance_matrix: ndarray[tuple[Any, ...], dtype[float64]]

Return a copy of the covariance matrix.

Returns:

A copy of the covariance matrix.

Return type:

numpy.ndarray

property mean_vector: ndarray[tuple[Any, ...], dtype[float64]]

Return a copy of the mean vector.

Returns:

A copy of the mean vector in x-then-p ordering.

Return type:

numpy.ndarray

plot_state(ax=None, n_std=2)

Plot single-mode phase-space marginals of the Gaussian state.

This method visualizes each mode of the Gaussian state by plotting its marginal covariance ellipse in the local phase space

(x_i, p_i),

where the full state uses the global x-then-p ordering

(x_1, …, x_n, p_1, …, p_n).

For each mode, the plotted ellipse represents an n_std-standard- deviation contour derived from the corresponding 2×2 marginal covariance matrix, and the mode mean is shown as a point.

Parameters:

  • ax (matplotlib.axes.Axes, optional) – The axes on which to draw the plot. If None, a new figure and axes are created.

  • n_std (Real, optional) – The number of standard deviations used to scale each covariance ellipse. Must be strictly positive. Default is 2.

Returns:

The axes containing the plot.

Return type:

matplotlib.axes.Axes

Raises:

  • TypeError – If n_std is not a real scalar.

  • ValueError – If n_std is not finite or not strictly positive.

Notes

Each mode is plotted independently using its local marginal mean and covariance. Inter-mode correlations are not visualized by this method.

The ellipse orientation is determined by the eigenvectors of the local covariance matrix, and the axis lengths are proportional to the square roots of its eigenvalues.

single_mode_displacement(displacement, mode_id)

Apply a single-mode displacement in place.

This method applies a phase-space displacement to one mode of the state, updating the internal mean vector and leaving the covariance matrix unchanged. The state uses the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

Parameters:

  • displacement (Complex) – The complex displacement amplitude. Must be finite.

  • mode_id (Integral) – The target mode index (1-based).

Returns:

The current state instance after in-place update.

Return type:

GaussianCVState

Raises:

  • TypeError – If displacement or mode_id has invalid type.

  • ValueError – If displacement is not finite or if mode_id is out of bounds.

Notes

This method mutates the current object and returns self to support chained operations.

single_mode_squeeze(squeeze_tuple, mode_id)

Apply a single-mode squeezing unitary in place.

This method applies a single-mode Gaussian squeezing transformation to one mode of the state, updating both the mean vector and covariance matrix in the x-then-p phase-space ordering.

Parameters:

  • squeeze_tuple (tuple of (Real, Real)) – The squeezing parameters (squeeze_magnitude, squeeze_angle).

  • mode_id (Integral) – The target mode index (1-based).

Returns:

The current state instance after in-place update.

Return type:

GaussianCVState

Raises:

  • TypeError – If the squeezing parameters or mode_id have invalid type.

  • ValueError – If the squeezing parameters are invalid or if mode_id is out of bounds.

Notes

This method mutates the current object and returns self to support chained operations.

single_mode_thermal_reset(nbar, mode_id)

Reset a single mode to an uncorrelated thermal state.

This method replaces the selected mode with a thermal state of occupation nbar by removing all correlations involving that mode and setting its local quadrature variances to

nbar + 1/2.

In the global x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n),

a single mode is represented by the two noncontiguous indices corresponding to x_i and p_i. This method therefore zeros the rows and columns associated with those indices and then sets the diagonal entries for that mode to nbar + 1/2.

Parameters:

  • nbar (Real) – The thermal occupation number. Must be non-negative and finite.

  • mode_id (Integral) – The target mode index (1-based).

Returns:

The current state instance after in-place update.

Return type:

GaussianCVState

Raises:

  • TypeError – If nbar or mode_id has invalid type.

  • ValueError – If nbar is negative or not finite, or if mode_id is out of bounds.

Notes

This is a destructive, non-unitary reset. It removes all correlations involving the selected mode and leaves the mean vector unchanged.

The state is modified in place and the method returns self to support chained operations.

state_to_vector()

Convert the Gaussian state to a flattened vector representation.

This method serializes the current state by concatenating the mean vector and the flattened covariance matrix into a one-dimensional array:

[mean_vector, covariance_matrix.flatten()]

where the covariance matrix is flattened in row-major (C) order.

Returns:

A one-dimensional array of length 2n + (2n)^2 representing the state.

Return type:

numpy.ndarray

Notes

The state is assumed to be ordered in the x-then-p convention

(x_1, …, x_n, p_1, …, p_n).

This representation is useful for numerical solvers such as solve_ivp or other vectorized workflows. Reconstruction can be performed using extract_mean_covariance.

classmethod thermal(n, nbars=None)

Construct an n-mode product thermal state.

Parameters:

  • n (Integral) – The number of modes. Must be strictly positive.

  • nbars (numpy.ndarray or None, optional) – A one-dimensional array of thermal occupation numbers. If fewer than n values are provided, the remaining modes are assigned zero thermal occupation. If None, all modes are vacuum.

Returns:

The n-mode thermal product state.

Return type:

GaussianCVState

Raises:

  • TypeError – If n or nbars has an invalid type.

  • ValueError – If n is invalid or if nbars contains invalid values.

two_mode_mix(coupling_tuple, mode_ids)

Apply a two-mode passive mixing unitary in place.

This method applies a two-mode passive Gaussian unitary to the selected pair of modes, updating both the mean vector and covariance matrix. The selected subsystem is ordered as

(x_i, x_j, p_i, p_j)

in the global x-then-p convention

(x_1, …, x_n, p_1, …, p_n).

Parameters:

  • coupling_tuple (tuple of (Real, Real)) – The mixing parameters (mixing_magnitude, mixing_angle).

  • mode_ids (tuple of (Integral, Integral)) – The pair of target mode indices (1-based). The two indices must be distinct.

Returns:

The current state instance after in-place update.

Return type:

GaussianCVState

Raises:

  • TypeError – If the coupling parameters or mode indices have invalid type.

  • ValueError – If the coupling parameters are invalid, if mode_ids does not have length 2, if the indices are not distinct, or if any mode index is out of bounds.

Notes

This method mutates the current object and returns self to support chained operations.

two_mode_squeeze(squeeze_tuple, mode_ids)

Apply a two-mode squeezing unitary in place.

This method applies a two-mode active Gaussian unitary to the selected pair of modes, updating both the mean vector and covariance matrix. The selected subsystem is ordered as

(x_i, x_j, p_i, p_j)

in the global x-then-p convention

(x_1, …, x_n, p_1, …, p_n).

Parameters:

  • squeeze_tuple (tuple of (Real, Real)) – The squeezing parameters (squeeze_magnitude, squeeze_angle).

  • mode_ids (tuple of (Integral, Integral)) – The pair of target mode indices (1-based). The two indices must be distinct.

Returns:

The current state instance after in-place update.

Return type:

GaussianCVState

Raises:

  • TypeError – If the squeezing parameters or mode indices have invalid type.

  • ValueError – If the squeezing parameters are invalid, if mode_ids does not have length 2, if the indices are not distinct, or if any mode index is out of bounds.

Notes

This method mutates the current object and returns self to support chained operations.

classmethod vacuum(n)

Construct an n-mode vacuum state.

Parameters:

n (Integral) – The number of modes. Must be strictly positive.

Returns:

The n-mode vacuum state.

Return type:

GaussianCVState

Raises:

  • TypeError – If n is not an integer.

  • ValueError – If n is not strictly positive.

gaussian_systems.initial_state.apply_1_mode_displacement(mean_vector_covariance_matrix_tuple, displacement, mode_id)

Apply a single-mode displacement to a Gaussian state.

This function applies a phase-space displacement to one mode of a Gaussian state represented by a mean vector and covariance matrix. In the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n),

a complex displacement alpha acts on the selected mode as

Δx = sqrt(2) * Re(alpha), Δp = sqrt(2) * Im(alpha).

The mean vector is shifted accordingly, while the covariance matrix is unchanged.

Parameters:

  • mean_vector_covariance_matrix_tuple (tuple of (numpy.ndarray, numpy.ndarray)) – The Gaussian state given as (mean_vector, covariance_matrix).

  • displacement (Complex) – The complex displacement amplitude. Must be finite.

  • mode_id (Integral) – The target mode index (1-based).

Returns:

The displaced Gaussian state as (new_mean_vector, covariance_matrix).

Return type:

tuple of (numpy.ndarray, numpy.ndarray)

Raises:

  • TypeError – If the state is invalid, if displacement is not complex-valued, or if mode_id is not an integer index.

  • ValueError – If displacement is not finite or if mode_id is outside the valid range.

Notes

Only the mean vector is modified. The covariance matrix is returned unchanged.

gaussian_systems.initial_state.apply_1_mode_squeeze_unitary(mean_vector_covariance_matrix_tuple, squeeze_magnitude_angle_tuple, mode_id)

Apply a single-mode squeezing unitary to a Gaussian state.

This function applies a single-mode squeezing transformation to one mode of an n-mode Gaussian state represented by a mean vector and covariance matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

If the squeezing parameters are (r, φ), the single-mode symplectic transformation is

S(r, φ) = R(φ/2) @ diag(exp(-r), exp(r)) @ R(φ/2)^T,

where R is the single-mode rotation matrix. This 2×2 transformation is embedded into the full 2n-dimensional phase space and applied as

mean -> S_full @ mean cov -> S_full @ cov @ S_full.T

Parameters:

  • mean_vector_covariance_matrix_tuple (tuple of (numpy.ndarray, numpy.ndarray)) – The Gaussian state given as (mean_vector, covariance_matrix).

  • squeeze_magnitude_angle_tuple (tuple of (Real, Real)) – The squeezing parameters (squeeze_magnitude, squeeze_angle). The squeezing magnitude must be non-negative and the squeezing angle must be finite.

  • mode_id (Integral) – The target mode index (1-based).

Returns:

The transformed Gaussian state as (new_mean_vector, new_covariance_matrix).

Return type:

tuple of (numpy.ndarray, numpy.ndarray)

Raises:

  • TypeError – If the state, squeezing parameters, or mode index have invalid types.

  • ValueError – If the state is invalid, if the squeezing magnitude is negative, if the squeezing angle is not finite, or if mode_id is outside the valid range.

Notes

The transformed covariance matrix is symmetrized numerically before being returned. Only the selected mode is squeezed; all other modes are left unchanged.

gaussian_systems.initial_state.apply_2_mode_mix_unitary(mean_vector_covariance_matrix_tuple, coupling_magnitude_angle_tuple, mode_ids)

Apply a two-mode passive mixing unitary to a Gaussian state.

This function applies a two-mode passive Gaussian unitary to an n-mode Gaussian state represented by a mean vector and covariance matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

If the mixing parameters are (r, φ), the corresponding two-mode symplectic transformation is the 4×4 real phase-space representation returned by two_mode_mixing_matrix(r, φ). This transformation is embedded into the full 2n-dimensional phase space and applied as

mean -> S_full @ mean cov -> S_full @ cov @ S_full.T

Parameters:

  • mean_vector_covariance_matrix_tuple (tuple of (numpy.ndarray, numpy.ndarray)) – The Gaussian state given as (mean_vector, covariance_matrix).

  • coupling_magnitude_angle_tuple (tuple of (Real, Real)) – The mixing parameters (mixing_magnitude, mixing_angle). The mixing magnitude must be non-negative and the mixing angle must be finite.

  • mode_ids (tuple of (Integral, Integral)) – The pair of mode indices (1-based) defining the two-mode subsystem. The two indices must be distinct.

Returns:

The transformed Gaussian state as (new_mean_vector, new_covariance_matrix).

Return type:

tuple of (numpy.ndarray, numpy.ndarray)

Raises:

  • TypeError – If the state, mixing parameters, or mode indices have invalid types.

  • ValueError – If the state is invalid, if the mixing magnitude is negative, if the mixing angle is not finite, if the tuple does not have length 2, if the indices are not distinct, or if any mode index is out of bounds.

Notes

The transformed covariance matrix is symmetrized numerically before being returned. The selected subsystem is ordered as

(x_i, x_j, p_i, p_j)

for mode_ids = (i, j).

gaussian_systems.initial_state.apply_2_mode_squeeze_unitary(mean_vector_covariance_matrix_tuple, squeeze_magnitude_angle_tuple, mode_ids)

Apply a two-mode squeezing unitary to a Gaussian state.

This function applies a two-mode active Gaussian unitary to an n-mode Gaussian state represented by a mean vector and covariance matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

If the squeezing parameters are (r, φ), the corresponding two-mode symplectic transformation is the 4×4 real phase-space representation returned by two_mode_squeezing_matrix(r, φ). This transformation is embedded into the full 2n-dimensional phase space and applied as

mean -> S_full @ mean cov -> S_full @ cov @ S_full.T

Parameters:

  • mean_vector_covariance_matrix_tuple (tuple of (numpy.ndarray, numpy.ndarray)) – The Gaussian state given as (mean_vector, covariance_matrix).

  • squeeze_magnitude_angle_tuple (tuple of (Real, Real)) – The squeezing parameters (squeeze_magnitude, squeeze_angle). The squeezing magnitude must be non-negative and the squeezing angle must be finite.

  • mode_ids (tuple of (Integral, Integral)) – The pair of mode indices (1-based) defining the two-mode subsystem. The two indices must be distinct.

Returns:

The transformed Gaussian state as (new_mean_vector, new_covariance_matrix).

Return type:

tuple of (numpy.ndarray, numpy.ndarray)

Raises:

  • TypeError – If the state, squeezing parameters, or mode indices have invalid types.

  • ValueError – If the state is invalid, if the squeezing magnitude is negative, if the squeezing angle is not finite, if mode_ids does not have length 2, if the two indices are not distinct, or if any mode index is outside the valid range.

Notes

The selected subsystem is ordered as

(x_i, x_j, p_i, p_j)

for mode_ids = (i, j). The transformed covariance matrix is symmetrized numerically before being returned.

gaussian_systems.initial_state.single_mode_squeeze_matrix(squeeze_magnitude, squeeze_angle)

Construct a single-mode squeezing matrix.

This function returns the 2×2 symplectic matrix representing a single-mode squeezing operation in the (x, p) quadrature basis.

The transformation is given by

S(r, φ) = R(φ/2) @ diag(exp(-r), exp(r)) @ R(φ/2)^T,

where r is the squeezing magnitude and φ is the squeezing angle. Here R(θ) is the phase-space rotation matrix.

Parameters:

  • squeeze_magnitude (Real) – The squeezing magnitude r. Must be non-negative.

  • squeeze_angle (Real) – The squeezing angle φ in radians. Must be finite.

Returns:

The single-mode squeezing matrix.

Return type:

numpy.ndarray of shape (2, 2)

Raises:

  • TypeError – If inputs are not real scalars.

  • ValueError – If squeeze_magnitude is negative or if inputs are not finite.

Notes

The matrix acts on phase-space vectors ordered as (x, p). For r = 0, the identity transformation is returned.

gaussian_systems.initial_state.thermal_vacuum_covariance(n, nbars=None)

Construct the covariance matrix of an n-mode thermal-vacuum product state.

This function returns the covariance matrix for an n-mode product state with thermal occupations nbars in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

For mode occupations nbar_i, the covariance matrix is diagonal with entries

nbar_i + 1/2

appearing once in the x block and once in the p block. Thus the returned matrix has the form

diag(v_1, …, v_n, v_1, …, v_n),

where v_i = nbar_i + 1/2.

Parameters:

  • n (Integral) – The number of modes. Must be strictly positive.

  • nbars (numpy.ndarray or None, optional) – A one-dimensional array of thermal occupation numbers for the modes. Entries must be real, finite, and non-negative. If None, all modes are taken to be in vacuum, corresponding to zero thermal occupation.

Returns:

The covariance matrix of the thermal-vacuum product state.

Return type:

numpy.ndarray of shape (2n, 2n)

Raises:

  • TypeError – If n is not an integer or if nbars is not a valid NumPy array when provided.

  • ValueError – If n is not strictly positive, if nbars contains invalid values, or if more than n thermal occupations are supplied.

Warns:

UserWarning – If fewer than n thermal occupations are provided. Missing modes are assigned zero thermal occupation.

Notes

This function constructs only the covariance matrix. The corresponding mean vector is zero for all modes.

gaussian_systems.initial_state.thermal_vacuum_mean(n)

Construct the mean vector of an n-mode thermal-vacuum state.

This function returns the mean vector for an n-mode Gaussian state with zero displacement. In the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n),

the mean vector is identically zero.

Parameters:

n (Integral) – The number of modes. Must be strictly positive.

Returns:

The zero mean vector.

Return type:

numpy.ndarray of shape (2n,)

Raises:

  • TypeError – If n is not an integer.

  • ValueError – If n is not strictly positive.

Notes

This corresponds to the mean vector of a thermal or vacuum state.

gaussian_systems.initial_state.two_mode_mixing_matrix(coupling_magnitude, coupling_angle)

Construct a two-mode passive mixing (beam-splitter) matrix.

This function returns the real symplectic representation of a two-mode passive Gaussian unitary. The transformation is defined by a complex 2×2 unitary matrix

U = [[cos(r), exp(i φ) sin(r)],

[-exp(-i φ) sin(r), cos(r)]],

where r is the mixing magnitude and φ is a phase.

This unitary is mapped to a real 4×4 symplectic matrix acting on the phase-space vector ordered as

(x_1, x_2, p_1, p_2),

via

S = [[Re(U), -Im(U)],

[Im(U), Re(U)]].

Parameters:

  • coupling_magnitude (Real) – The mixing magnitude r. Must be non-negative.

  • coupling_angle (Real) – The phase φ in radians. Must be finite.

Returns:

The two-mode mixing symplectic matrix.

Return type:

numpy.ndarray of shape (4, 4)

Raises:

  • TypeError – If inputs are not real scalars.

  • ValueError – If coupling_magnitude is negative or if inputs are not finite.

Notes

This transformation is passive (number-conserving) and corresponds to a beam-splitter–type interaction between two modes.

gaussian_systems.initial_state.two_mode_squeezing_matrix(squeeze_magnitude, squeeze_angle)

Construct a two-mode squeezing (active) symplectic matrix.

This function returns the 4×4 real symplectic matrix representing a two-mode squeezing transformation in the x-then-p phase-space ordering

(x_1, x_2, p_1, p_2).

The transformation is parameterized by a squeezing magnitude r and squeezing angle φ, with

μ = cosh(r), ν = sinh(r) * exp(i φ).

Writing ν = ν_r + i ν_i, the resulting matrix is

[[ μ, ν_r, 0, ν_i],

[ν_r, μ, ν_i, 0 ], [ 0, ν_i, μ, -ν_r], [ν_i, 0, -ν_r, μ ]].

Parameters:

  • squeeze_magnitude (Real) – The squeezing magnitude r. Must be non-negative.

  • squeeze_angle (Real) – The squeezing angle φ in radians. Must be finite.

Returns:

The two-mode squeezing symplectic matrix.

Return type:

numpy.ndarray of shape (4, 4)

Raises:

  • TypeError – If inputs are not real scalars.

  • ValueError – If squeeze_magnitude is negative or if inputs are not finite.

Notes

This transformation is active (non-number-conserving) and generates correlations and entanglement between the two modes.

Systems

class gaussian_systems.systems.GaussianCVSystem(hamiltonian_matrix, lindblad_matrix)

Bases: object

Container for an n-mode Gaussian continuous-variable dynamical system.

This class represents a Gaussian CV system by a quadratic Hamiltonian matrix and a Lindblad Gram matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The Hamiltonian matrix defines the coherent quadratic dynamics, while the Lindblad Gram matrix defines the dissipative Gaussian channel structure.

Parameters:

  • hamiltonian_matrix (numpy.ndarray) – The quadratic Hamiltonian matrix of shape (2n, 2n).

  • lindblad_matrix (numpy.ndarray) – The Lindblad Gram matrix of shape (2n, 2n).

Raises:

  • TypeError – If the inputs are not valid NumPy arrays.

  • ValueError – If the Hamiltonian and Lindblad matrices are invalid or dimensionally incompatible.

Notes

This class stores the system generator, not the state. Time evolution of Gaussian states is defined relative to both matrices.

beamsplitter_coupling(subsystem, coupling)

Add a beamsplitter-type coupling to the Hamiltonian.

This adds the quadratic interaction

H += g (x_i x_j + p_i p_j),

which corresponds to a number-conserving bilinear coupling between modes.

Parameters:

  • subsystem (tuple of (Integral, Integral)) – The pair of mode indices (1-based).

  • coupling (Real) – The coupling strength.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

copy_system()

Create a deep copy of the Gaussian CV system.

This method returns a new GaussianCVSystem instance with copies of the Hamiltonian and Lindblad matrices. The returned system is fully independent of the original.

Returns:

A new system instance with identical parameters.

Return type:

GaussianCVSystem

Notes

The Hamiltonian and Lindblad matrices are copied, so modifying the returned system does not affect the original.

evolve_state(state, t_eval)

Evolve a Gaussian state over a specified time grid.

This method evolves the supplied GaussianCVState under the current GaussianCVSystem at each time in t_eval using the closed-form Gaussian channel generated from the system drift and diffusion matrices.

The evolution is performed in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n),

with the mean and covariance propagated independently through the precompiled Gaussian channel representation.

Parameters:

  • state (GaussianCVState) – The initial Gaussian state to evolve.

  • t_eval (numpy.ndarray) – A one-dimensional array of evaluation times. Entries must be real, finite, non-negative, non-empty, and strictly increasing.

Returns:

A solution object containing the evaluation times, evolved mean vectors, and evolved covariance matrices.

Return type:

GaussianSolution

Raises:

  • TypeError – If state or t_eval has invalid type.

  • ValueError – If the state and system are incompatible, if t_eval is invalid, or if any evolved covariance matrix fails the quantum physicality test.

Notes

The input state is not modified. The covariance evolution is implemented using an augmented vectorized representation of the form

(vec(covariance), 1),

so that affine covariance evolution can be handled by a linear matrix exponential.

classmethod free_evolution(n, frequency_array=None)

Construct a free-evolution Gaussian system with no dissipation.

This classmethod returns an n-mode Gaussian system whose Hamiltonian consists only of diagonal self-energy terms and whose Lindblad Gram matrix is identically zero.

Parameters:

  • n (Integral) – The number of modes. Must be strictly positive.

  • frequency_array (numpy.ndarray or None, optional) – A one-dimensional array of mode frequencies. If fewer than n frequencies are supplied, missing entries are set to zero according to the rotated-frame convention. If None, all frequencies are taken to be zero.

Returns:

The free-evolution system.

Return type:

GaussianCVSystem

Raises:

  • TypeError – If n or frequency_array has invalid type.

  • ValueError – If the frequency specification is invalid.

gaussian_channel()

Construct the linear generators for Gaussian moment evolution.

This method returns the matrices governing the evolution of the mean vector and the vectorized covariance matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The mean evolves as

d/dt mean = A @ mean,

where A is the drift matrix.

The covariance evolves according to

d/dt vec(cov) = K @ vec(cov) + vec(D),

where

K = A ⊗ I + I ⊗ A,

and D is the diffusion matrix. This affine evolution is embedded into a linear system by augmenting the state vector as

(vec(cov), 1),

yielding

d/dt [vec(cov); 1] = A_covariance @ [vec(cov); 1].

Parameters:

None

Return type:

tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]]

Returns:

  • tuple of (numpy.ndarray, numpy.ndarray)

  • The pair (A, A_covariance), where

    • A has shape (2n, 2n) and governs mean evolution

    • A_covariance has shape (4n^2 + 1, 4n^2 + 1) and governs the

  • augmented covariance evolution

Notes

The augmented covariance generator encodes both drift and diffusion in a linear form suitable for matrix-exponential propagation.

generate_drift_and_diffusion()

Generate the Gaussian drift and diffusion matrices.

This method constructs the matrices governing first- and second-moment evolution for the current Gaussian CV system in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The returned matrices (A, D) define the Gaussian moment equations

d/dt mean = A @ mean, d/dt cov = A @ cov + cov @ A.T + D.

They are computed from the system Hamiltonian matrix H, Lindblad Gram matrix M, and canonical symplectic form Ω as

A = Ω @ (H + Im(M)), D = Ω @ Re(M) @ Ω.T.

Parameters:

None

Returns:

The drift matrix A and diffusion matrix D, each of shape (2n, 2n).

Return type:

tuple of (numpy.ndarray, numpy.ndarray)

Notes

Any residual imaginary parts arising from numerical roundoff are removed with np.real before returning the result.

property hamiltonian_matrix: ndarray[tuple[Any, ...], dtype[float64]]

Return a copy of the Hamiltonian matrix.

Returns:

A copy of the quadratic Hamiltonian matrix.

Return type:

numpy.ndarray

property lindblad_matrix: ndarray[tuple[Any, ...], dtype[complex128]]

Return a copy of the Lindblad Gram matrix.

Returns:

A copy of the Lindblad Gram matrix.

Return type:

numpy.ndarray

momentum_coupling(subsystem, coupling)

Add a momentum-momentum coupling term to the Hamiltonian.

This adds a quadratic interaction of the form

H += g p_i p_j

Parameters:

  • subsystem (tuple of (Integral, Integral)) – The pair of mode indices (1-based).

  • coupling (Real) – The coupling strength.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

multi_annihilation_dissipator(subsystem, decay)

Add a collective annihilation dissipator to the system.

This method adds a single-environment dissipative channel with collective jump operator proportional to

L ∝ sum_i sqrt(decay) a_i,

where the sum runs over the specified subsystem modes.

Parameters:

  • subsystem (tuple of Integral) – The mode indices (1-based) included in the collective dissipator.

  • decay (Real) – The decay rate associated with each participating mode.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

Raises:

  • TypeError – If the inputs have invalid type.

  • ValueError – If the subsystem indices are invalid or if the decay rate is invalid.

Notes

This is a collective dissipator: all selected modes couple to the same environmental channel.

multi_position_dissipator(subsystem, decay)

Add a collective position dissipator to the system.

This method adds a single-environment dissipative channel with collective jump operator proportional to

L ∝ sum_i sqrt(decay) x_i,

where the sum runs over the specified subsystem modes.

Parameters:

  • subsystem (tuple of Integral) – The mode indices (1-based) included in the collective dissipator.

  • decay (Real) – The decay rate associated with each participating mode.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

Raises:

  • TypeError – If the inputs have invalid type.

  • ValueError – If the subsystem indices are invalid or if the decay rate is invalid.

Notes

This is a collective dissipator: all selected modes couple to the same environmental channel.

multi_thermal_dissipator(subsystem, decay, thermal_occupation)

Add a collective thermal dissipator to the system.

This method adds a collective thermal bath acting on the specified modes, with emission and absorption contributions corresponding to thermal occupation thermal_occupation. The resulting dissipator is built from two independent collective channels:

L_emission ∝ sum_i sqrt(decay * (thermal_occupation + 1)) a_i, L_absorption ∝ sum_i sqrt(decay * thermal_occupation) a_i^†.

Parameters:

  • subsystem (tuple of Integral) – The mode indices (1-based) included in the collective dissipator.

  • decay (Real) – The base decay rate.

  • thermal_occupation (Real) – The bath thermal occupation number. Must be non-negative.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

Raises:

  • TypeError – If the inputs have invalid type.

  • ValueError – If the subsystem indices are invalid or if decay or thermal_occupation is invalid.

Notes

This is a collective thermal dissipator: all selected modes couple to the same thermal environment. Emission and absorption are added as separate Lindblad channels.

position_coupling(subsystem, coupling)

Add a position-position coupling term to the Hamiltonian.

This adds a quadratic interaction of the form

H += g x_i x_j

in the x-then-p phase-space ordering.

Parameters:

  • subsystem (tuple of (Integral, Integral)) – The pair of mode indices (1-based).

  • coupling (Real) – The coupling strength.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

position_difference_coupling(subsystem, coupling)

Add a position-difference quadratic term to the Hamiltonian.

This adds the interaction

H += g (x_i - x_j)^2

= g (x_i^2 + x_j^2 - 2 x_i x_j),

implemented through a combination of quadratic terms.

Parameters:

  • subsystem (tuple of (Integral, Integral)) – The pair of mode indices (1-based).

  • coupling (Real) – The coupling strength.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

position_i_momentum_j_coupling(subsystem, coupling)

Add a cross-quadrature coupling term to the Hamiltonian.

This adds a quadratic interaction of the form

H += g x_i p_j

Parameters:

  • subsystem (tuple of (Integral, Integral)) – The ordered pair of mode indices (1-based).

  • coupling (Real) – The coupling strength.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

squeezer_coupling(subsystem, coupling)

Add a two-mode squeezing coupling to the Hamiltonian.

This adds the quadratic interaction

H += g (x_i x_j - p_i p_j),

corresponding to a non-number-conserving two-mode squeezing interaction.

Parameters:

  • subsystem (tuple of (Integral, Integral)) – The pair of mode indices (1-based).

  • coupling (Real) – The coupling strength.

Returns:

The updated system (self).

Return type:

GaussianCVSystem

gaussian_systems.systems.single_pole_ou_embedding(state, system, subsystem, coupling_types, memory_rate, env_freq, decay_rate, thermal_occupation=0)

Construct the single-pole Ornstein–Uhlenbeck pseudomode embedding.

This function embeds an n-mode Gaussian state and system into an (n+1)-mode representation by appending a pseudomode as mode n+1. The phase-space convention is x-then-p ordering,

(x_1, …, x_n, x_{n+1}, p_1, …, p_n, p_{n+1}).

The embedding implements a single-pole Ornstein–Uhlenbeck (OU) environment via a Markovian pseudomode model. The appended mode represents the environmental memory degree of freedom.

Construction steps

  • Embed the mean vector, covariance matrix, Hamiltonian, and Lindblad Gram matrix into the enlarged (n+1)-mode phase space.

  • Initialize the pseudomode in a thermal state with occupation thermal_occupation.

  • Add a thermal dissipator acting on the pseudomode with decay rate

    γ_pseudo = 2 * memory_rate.

  • Add system–pseudomode interaction terms derived from the rotating-wave approximation (RWA) interaction

    H_int = L c^† + L^† c,

    where c is the pseudomode annihilation operator and L is the system coupling operator.

Coupling conventions

Each entry of coupling_types specifies the form of the system operator L for the corresponding mode in subsystem:

  • "position":

    L = x_i

    H_int = x_i (c^† + c) = sqrt(2) x_i x_{n+1}

    Implemented as a position-position coupling with an additional factor of sqrt(2) applied to the coupling strength.

  • "momentum":

    L = p_i

    H_int = p_i (c^† + c) = sqrt(2) p_i x_{n+1}

    Implemented as a cross-quadrature coupling between p_i and the pseudomode position x_{n+1}, with an additional factor of sqrt(2).

  • "annihilation":

    L = a_i

    H_int = a_i c^† + a_i^† c

    Implemented as a beamsplitter-type coupling

    x_i x_{n+1} + p_i p_{n+1}.

The effective coupling strength is

g = sqrt(decay_rate * memory_rate / 2),

and the sqrt(2) factor required for Hermitian couplings (position and momentum) is applied explicitly in those branches.

The pseudomode self-energy is added with frequency env_freq.

type state:

GaussianCVState

param state:

The initial n-mode Gaussian state.

type state:

GaussianCVState

type system:

GaussianCVSystem

param system:

The corresponding n-mode Gaussian system.

type system:

GaussianCVSystem

type subsystem:

tuple[Integral, ...]

param subsystem:

The system mode indices (1-based) coupled to the pseudomode.

type subsystem:

tuple of Integral

type coupling_types:

tuple[str, ...]

param coupling_types:

The coupling type for each subsystem mode. Must contain only "position", "momentum", or "annihilation".

type coupling_types:

tuple of str

type memory_rate:

Real

param memory_rate:

The OU memory rate. Must be strictly positive.

type memory_rate:

Real

type env_freq:

Real

param env_freq:

The pseudomode frequency.

type env_freq:

Real

type decay_rate:

Real

param decay_rate:

The effective system–environment coupling rate.

type decay_rate:

Real

type thermal_occupation:

Real

param thermal_occupation:

The thermal occupation number of the environment. Must be non-negative.

type thermal_occupation:

Real

returns:

The embedded Gaussian state and system in the enlarged (n+1)-mode space.

rtype:

tuple of (GaussianCVState, GaussianCVSystem)

raises TypeError:

If the inputs have invalid types.

raises ValueError:

If the state and system are incompatible, if the embedding parameters are invalid, or if any coupling type is not recognized.

Notes

The appended pseudomode is always mode n+1. The embedding produces a Markovian system whose reduced dynamics reproduce a single-pole OU environment for the specified subsystem couplings.

Solution

class gaussian_systems.solution.GaussianSolution(t_eval, mean_vectors, covariance_matrices)

Bases: object

Container for the time-resolved output of Gaussian CV evolution.

This class stores the time grid together with the corresponding mean vectors and covariance matrices produced by GaussianCVSystem.evolve_state. It also provides convenience methods for evaluating Gaussian CV metrics along the trajectory.

The phase-space convention throughout is x-then-p ordering,

(x_1, …, x_n, p_1, …, p_n).

Parameters:

  • t_eval (array-like) – The time grid used for evolution.

  • mean_vectors (sequence of numpy.ndarray) – The evolved mean vectors at each time in t_eval.

  • covariance_matrices (sequence of numpy.ndarray) – The evolved covariance matrices at each time in t_eval.

Notes

This class is primarily intended as the output container of GaussianCVSystem.evolve_state.

entanglement_time_trace(subsystem=(1, 2))

Compute logarithmic negativity along the trajectory.

This method evaluates the logarithmic negativity of the specified two-mode subsystem for each covariance matrix in the trajectory.

Parameters:

subsystem (tuple of (int, int), optional) – The two mode indices (1-based) defining the subsystem. Default is (1, 2).

Returns:

The logarithmic negativity evaluated at each time point.

Return type:

list of Real

Notes

This method uses compute_logarithmic_negativity on each covariance matrix in the trajectory.

entropy_time_trace(subsystem=None)

Compute Rényi-2 entropy along the trajectory.

This method evaluates the Rényi-2 entropy of either the full state or a selected subsystem for each covariance matrix in the trajectory.

Parameters:

subsystem (tuple of int or None, optional) – The subsystem mode indices (1-based). If None, the entropy of the full state is returned.

Returns:

The Rényi-2 entropy evaluated at each time point.

Return type:

list of Real

fidelity_time_trace_fixed(state2, subsystem=None)

Compute fidelity to a fixed Gaussian state along the trajectory.

This method evaluates the Gaussian fidelity between each evolved state in the trajectory and a fixed comparison state state2. If a subsystem is specified, the same subsystem reduction is applied to both the trajectory state and the fixed state before computing fidelity.

Parameters:

  • state2 (GaussianCVState) – The fixed comparison Gaussian state.

  • subsystem (tuple of int or None, optional) – The subsystem mode indices (1-based). If None, fidelity is computed using the full states.

Returns:

The Gaussian fidelity evaluated at each time point.

Return type:

list of Real

Notes

This method assumes that the fixed state and each evolved state are dimensionally compatible after applying the same subsystem selection.

purity_time_trace(subsystem=None)

Compute purity along the trajectory.

This method evaluates the Gaussian purity of either the full state or a selected subsystem for each covariance matrix in the trajectory.

Parameters:

subsystem (tuple of int or None, optional) – The subsystem mode indices (1-based). If None, the purity of the full state is returned.

Returns:

The purity evaluated at each time point.

Return type:

list of Real

Metrics

gaussian_systems.metrics.compute_gaussian_fidelity(mean_covariance_tuple1, mean_covariance_tuple2)

Compute the standard Uhlmann fidelity between two Gaussian states.

This function evaluates the Gaussian-state fidelity between two states specified by their mean vectors and covariance matrices in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The function automatically selects the appropriate implementation based on the number of modes:

  • n = 1: one-mode closed-form formula

  • n = 2: two-mode closed-form formula

  • n >= 3: general n-mode auxiliary-matrix formula

The returned value is the standard Uhlmann fidelity, not the root fidelity convention often used in parts of the Gaussian-state literature.

Parameters:

  • mean_covariance_tuple1 (tuple of (numpy.ndarray, numpy.ndarray)) – The first Gaussian state, given as (mean_vector, covariance_matrix).

  • mean_covariance_tuple2 (tuple of (numpy.ndarray, numpy.ndarray)) – The second Gaussian state, given as (mean_vector, covariance_matrix).

Returns:

The standard Uhlmann fidelity between the two Gaussian states.

Return type:

numbers.Real

Raises:

  • TypeError – If either input state has invalid type.

  • ValueError – If either state is invalid or if the two states are dimensionally incompatible.

Notes

This function dispatches to specialized low-mode formulas when available and otherwise uses the general n-mode Gaussian fidelity expression.

gaussian_systems.metrics.compute_logarithmic_negativity(covariance_matrix, subsystem=(1, 2))

Compute the logarithmic negativity of a two-mode Gaussian subsystem.

This function evaluates the logarithmic negativity of the two-mode subsystem specified by subsystem from a physical Gaussian covariance matrix in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The selected two-mode covariance matrix is ordered as

(x_i, x_j, p_i, p_j),

where subsystem = (i, j). Logarithmic negativity is computed from the smallest symplectic eigenvalue of the partially transposed covariance matrix:

E_N = max(0, -log(2 * nu_min_tilde)),

where nu_min_tilde is the minimum symplectic eigenvalue after partial transposition.

Parameters:

  • covariance_matrix (numpy.ndarray) – A physical Gaussian covariance matrix of shape (2n, 2n).

  • subsystem (tuple of (int, int), optional) – The two mode indices (1-based) defining the subsystem. Default is (1, 2).

Returns:

The logarithmic negativity of the selected two-mode subsystem.

Return type:

numbers.Real

Raises:

  • TypeError – If subsystem does not have the correct type.

  • ValueError – If covariance_matrix is not physical, if subsystem does not have length 2, or if the subsystem indices are invalid.

Notes

Partial transposition is implemented in the reduced two-mode x-then-p ordering via the fixed matrix

diag(1, 1, -1, 1),

corresponding to a sign flip of one momentum quadrature. The formula assumes the convention that vacuum quadrature variance equals 1/2.

gaussian_systems.metrics.n_mode_gaussian_fidelity(mean_covariance_tuple1, mean_covariance_tuple2)

Compute the standard Uhlmann fidelity between two n-mode Gaussian states.

This function evaluates the fidelity between two Gaussian states specified by their mean vectors and covariance matrices in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

The implementation uses the general auxiliary-matrix formula for the root Gaussian fidelity prefactor,

F_0(V1, V2) = F_tot / det(V1 + V2)^(1/4),

where

F_tot = prod_k [w_k^aux + sqrt((w_k^aux)^2 - 1)]^(1/2),

and the auxiliary matrix is defined by

W = -2 V i Ω, W_aux = -(W1 + W2)^(-1) (I + W2 W1).

The returned value is the standard Uhlmann fidelity, i.e. the square of the root fidelity convention often used in the Gaussian-state literature.

Parameters:

  • mean_covariance_tuple1 (tuple of (numpy.ndarray, numpy.ndarray)) – The first Gaussian state, given as (mean_vector, covariance_matrix).

  • mean_covariance_tuple2 (tuple of (numpy.ndarray, numpy.ndarray)) – The second Gaussian state, given as (mean_vector, covariance_matrix).

Returns:

The standard Uhlmann fidelity between the two Gaussian states.

Return type:

numbers.Real

Raises:

  • TypeError – If either input state has invalid type.

  • ValueError – If either state is invalid, if the two states are dimensionally incompatible, or if the auxiliary spectrum falls outside the expected physical domain.

Notes

The calculation is assembled in the log-domain as

log F = 2 log(F_tot) - (1/2) log(det(V1 + V2))

  • (1/2) du^T (V1 + V2)^(-1) du,

and exponentiated only at the end for improved numerical stability.

gaussian_systems.metrics.one_mode_gaussian_fidelity(mean_covariance_tuple1, mean_covariance_tuple2)

Compute the standard Uhlmann fidelity between two one-mode Gaussian states.

This function evaluates the fidelity between two one-mode Gaussian states specified by their mean vectors and covariance matrices in the (x, p) phase-space ordering.

The implementation uses the one-mode closed-form Gaussian fidelity formula in terms of

Delta = det(V1 + V2), Lambda = 2^(2n) det(V1 + (i/2) Ω) det(V2 + (i/2) Ω),

with n = 1. The returned value is the standard Uhlmann fidelity, i.e. the square of the root fidelity often used in the Gaussian-state literature.

The denominator is evaluated in rationalized form,

1 / (sqrt(Delta + Lambda) - sqrt(Lambda)) = (sqrt(Delta + Lambda) + sqrt(Lambda)) / Delta,

which is the form used internally.

Parameters:

  • mean_covariance_tuple1 (tuple of (numpy.ndarray, numpy.ndarray)) – The first one-mode Gaussian state, given as (mean_vector, covariance_matrix).

  • mean_covariance_tuple2 (tuple of (numpy.ndarray, numpy.ndarray)) – The second one-mode Gaussian state, given as (mean_vector, covariance_matrix).

Returns:

The standard Uhlmann fidelity between the two one-mode Gaussian states.

Return type:

numbers.Real

Raises:

  • TypeError – If either input state has invalid type.

  • ValueError – If either state is invalid, if the two states are dimensionally incompatible, or if either covariance matrix is not 2x2.

Notes

This function is restricted to one-mode states. The calculation is assembled partly in the log-domain for improved numerical stability, then exponentiated at the end.

gaussian_systems.metrics.renyi_two_entropy(covariance_matrix, subsystem=None)

Compute the Rényi-2 entropy of a Gaussian state or subsystem.

This function evaluates the Rényi-2 entropy

S_2 = -log(Tr(rho^2)),

where Tr(rho^2) is the purity of the state. The purity is computed from the covariance matrix under the convention that vacuum quadrature variance equals 1/2.

The covariance matrix is assumed to follow the x-then-p ordering

(x_1, …, x_n, p_1, …, p_n).

If a subsystem is specified, the entropy is computed for the corresponding reduced Gaussian state.

Parameters:

  • covariance_matrix (numpy.ndarray) – A physical Gaussian covariance matrix of shape (2n, 2n).

  • subsystem (tuple of int or None, optional) – The mode indices (1-based) defining the subsystem. If None, the entropy of the full state is returned.

Returns:

The Rényi-2 entropy of the selected state or subsystem.

Return type:

numbers.Real

Raises:

ValueError – If covariance_matrix is not physical or if the subsystem indices are invalid.

Notes

The Rényi-2 entropy is non-negative and equals zero for pure Gaussian states.

gaussian_systems.metrics.state_purity(covariance_matrix, subsystem=None)

Compute the Gaussian purity of a state or subsystem.

This function evaluates the purity

purity = Tr(rho^2),

for a Gaussian state specified by its covariance matrix. In the convention where vacuum quadrature variance equals 1/2, the purity of a subsystem with covariance matrix V is

purity = 1 / sqrt(det(2 V)).

The full covariance matrix is assumed to follow the x-then-p ordering

(x_1, …, x_n, p_1, …, p_n).

If a subsystem is specified, the corresponding reduced covariance matrix is extracted in the ordering

(x_{i_1}, …, x_{i_k}, p_{i_1}, …, p_{i_k}).

Parameters:

  • covariance_matrix (numpy.ndarray) – A physical Gaussian covariance matrix of shape (2n, 2n).

  • subsystem (tuple of int or None, optional) – The mode indices (1-based) defining the subsystem whose purity is to be computed. If None, the purity of the full state is returned.

Returns:

The purity of the selected state or subsystem.

Return type:

numbers.Real

Raises:

ValueError – If covariance_matrix is not physical or if the subsystem indices are invalid.

Notes

For a pure Gaussian state, the purity equals 1. Mixed states satisfy 0 < purity < 1.

gaussian_systems.metrics.two_mode_gaussian_fidelity(mean_covariance_tuple1, mean_covariance_tuple2)

Compute the standard Uhlmann fidelity between two two-mode Gaussian states.

This function evaluates the fidelity between two two-mode Gaussian states specified by their mean vectors and covariance matrices in the x-then-p phase-space ordering

(x_1, x_2, p_1, p_2).

The implementation uses the closed-form two-mode Gaussian fidelity prefactor

F_0^2(V1, V2) = 1 / (sqrt(Gamma) + sqrt(Lambda) - sqrt((sqrt(Gamma) + sqrt(Lambda))^2 - Delta)),

where

Delta = det(V1 + V2), Lambda = 2^(2n) det(V1 + (i/2) Ω) det(V2 + (i/2) Ω), Gamma = 2^(2n) det(Ω V1 Ω V2 - (1/4) I),

with n = 2.

The denominator is evaluated in rationalized form,

F_0^2(V1, V2) = (sqrt(Gamma) + sqrt(Lambda) + sqrt((sqrt(Gamma) + sqrt(Lambda))^2 - Delta)) / Delta,

and combined with the Gaussian displacement contribution. The returned value is the standard Uhlmann fidelity, not the root fidelity.

Parameters:

  • mean_covariance_tuple1 (tuple of (numpy.ndarray, numpy.ndarray)) – The first two-mode Gaussian state, given as (mean_vector, covariance_matrix).

  • mean_covariance_tuple2 (tuple of (numpy.ndarray, numpy.ndarray)) – The second two-mode Gaussian state, given as (mean_vector, covariance_matrix).

Returns:

The standard Uhlmann fidelity between the two two-mode Gaussian states.

Return type:

numbers.Real

Raises:

  • TypeError – If either input state has invalid type.

  • ValueError – If either state is invalid, if the two states are dimensionally incompatible, or if either covariance matrix is not 4x4.

Notes

This function is restricted to two-mode states. The calculation is assembled partly in the log-domain for improved numerical stability, then exponentiated at the end.

Conventions

gaussian_systems.conventions.compress_mean_covariance(mean_vector, covariance_matrix)

Flatten a Gaussian state into a single vector.

This function concatenates a mean vector and covariance matrix into a one-dimensional array. The output is given by

[mean_vector, covariance_matrix.flatten()],

where the covariance matrix is flattened in row-major (C) order.

Parameters:

  • mean_vector (numpy.ndarray)

  • (2n (The covariance matrix of shape)

  • ).

  • covariance_matrix (numpy.ndarray)

  • (2n

  • 2n).

Returns:

A one-dimensional array of length 2n + (2n)^2 containing the concatenated data.

Return type:

numpy.ndarray

Raises:

  • TypeError – If inputs are not NumPy arrays.

  • ValueError – If the mean vector and covariance matrix are not valid or dimensionally consistent.

Notes

The covariance matrix is flattened in row-major (C) order. To reconstruct the original matrix, the number of modes n (or dimension 2n) must be known.

gaussian_systems.conventions.covariance_subsystem(covariance_matrix, indices=None)

Extract a subsystem covariance matrix in x-then-p ordering.

This function returns the covariance matrix corresponding to a subset of modes from an n-mode Gaussian state. The input covariance matrix is assumed to be ordered as (x_1, …, x_n, p_1, …, p_n).

For a subsystem with modes (i_1, …, i_k), the returned matrix corresponds to the quadratures (x_{i_1}, …, x_{i_k}, p_{i_1}, …, p_{i_k}).

Parameters:

  • covariance_matrix (numpy.ndarray) – The full covariance matrix of shape (2n, 2n).

  • indices (tuple of Integral or None) – The mode indices (1-based) defining the subsystem.

Returns:

The subsystem covariance matrix of shape (2k, 2k). If indices is None returns a copy of the covariance matrix.

Return type:

numpy.ndarray

Raises:

  • TypeError – If covariance_matrix is not a NumPy array.

  • ValueError – If covariance_matrix is not valid or if indices are out of bounds.

Notes

The returned matrix follows the (x-subsystem, p-subsystem) ordering.

gaussian_systems.conventions.extract_mean_covariance(mean_covariance_vector)

Reconstruct a Gaussian state from a flattened representation.

This function extracts a mean vector and covariance matrix from a one-dimensional array of length 2n + (2n)^2, assumed to be constructed via compress_mean_covariance.

The covariance matrix is reshaped from the flattened data and then symmetrized as (A + A.conj().T) / 2 before validation.

Parameters:

mean_covariance_vector (numpy.ndarray) – The flattened representation of the Gaussian state.

Returns:

The reconstructed mean vector of shape (2n,) and covariance matrix of shape (2n, 2n).

Return type:

tuple of (numpy.ndarray, numpy.ndarray)

Raises:

  • TypeError – If the input is not a real-valued NumPy vector.

  • ValueError – If the length is incompatible with 2n + 4n^2, or if the reconstructed mean vector and covariance matrix are not valid.

Notes

The covariance matrix is assumed to be flattened in row-major (C) order. The reconstruction enforces symmetry via projection, so the operation is not an exact inverse of compress_mean_covariance if the input matrix was not symmetric.

gaussian_systems.conventions.index_list(n, indices)

Construct phase-space indices for a subsystem in x-then-p ordering.

This function returns the indices corresponding to a subsystem defined by indices in an n-mode phase-space vector ordered as $(x_1, …, x_n, p_1, …, p_n)$.

For a subsystem with modes $(i_1, …, i_k)$, the returned indices correspond to $(x_{i_1}, …, x_{i_k}, p_{i_1}, …, p_{i_k})$, preserving the input ordering.

Parameters:

  • n (Integral) – The total number of modes in the system. Must be strictly positive.

  • indices (tuple of Integral) – The mode indices (1-based) defining the subsystem.

Returns:

The zero-based indices of the subsystem in phase space, ordered as x-components followed by p-components.

Return type:

numpy.ndarray of int

Raises:

  • TypeError – If inputs have invalid types.

  • ValueError – If any index is outside the range [1, n].

Notes

The returned ordering is (x-subsystem, p-subsystem), not interleaved. No sorting or deduplication of indices is performed.

gaussian_systems.conventions.is_physical_covariance_matrix(covariance_matrix, tol=1e-08)

Test whether a covariance matrix is physically admissible for a Gaussian state.

This function checks whether a covariance matrix satisfies the Gaussian quantum uncertainty relation

V + (i/2) Ω >= 0,

where Ω is the canonical symplectic form for an n-mode system in the x-then-p phase-space ordering

(x_1, …, x_n, p_1, …, p_n).

Parameters:

  • covariance_matrix (numpy.ndarray) – The covariance matrix of shape (2n, 2n).

  • tol (Real, optional) – Numerical tolerance used in the positivity test. Eigenvalues greater than or equal to -tol are accepted as non-negative. Default is 1e-8.

Returns:

True if the covariance matrix satisfies the quantum uncertainty relation within tolerance, and False otherwise.

Return type:

bool

Raises:

  • TypeError – If covariance_matrix is not a NumPy array or if tol is not a real scalar.

  • ValueError – If covariance_matrix is not a valid covariance matrix or if tol is not strictly positive.

Notes

The covariance matrix is symmetrized internally as (V + V.conj().T) / 2 before testing physicality. This function checks quantum admissibility, not just classical positive semidefiniteness.

gaussian_systems.conventions.mean_subsystem(mean_vector, indices=None)

Extract a subsystem mean vector in x-then-p ordering.

This function returns the mean vector corresponding to a subset of modes from an n-mode Gaussian state. The input mean vector is assumed to be ordered as (x_1, …, x_n, p_1, …, p_n).

For a subsystem with modes (i_1, …, i_k), the returned vector is ordered as (x_{i_1}, …, x_{i_k}, p_{i_1}, …, p_{i_k}).

Parameters:

  • mean_vector (numpy.ndarray) – The full mean vector of shape (2n,).

  • indices (tuple of Integral or None) – The mode indices (1-based) defining the subsystem.

Returns:

The subsystem mean vector of shape (2k,). If indices is None, returns a copy of the original mean vector

Return type:

numpy.ndarray

Raises:

  • TypeError – If mean_vector is not a NumPy array.

  • ValueError – If mean_vector is not valid or if indices are out of bounds.

Notes

The returned ordering is (x-subsystem, p-subsystem), not interleaved.

gaussian_systems.conventions.require_physical_covariance(covariance_matrix)

Enforce that a covariance matrix is physically admissible.

This function validates that covariance_matrix is both: 1. a valid classical covariance matrix (real, symmetric, positive semidefinite), 2. a valid quantum Gaussian covariance matrix satisfying the Heisenberg uncertainty relation.

Parameters:

covariance_matrix (numpy.ndarray) – The covariance matrix of shape (2n, 2n).

Raises:

  • TypeError – If covariance_matrix is not a NumPy array.

  • ValueError – If covariance_matrix is not a valid covariance matrix or fails the Heisenberg uncertainty condition.

Return type:

None

Notes

This function performs a strict validation and raises an exception on failure. It uses the default numerical tolerance of is_physical_covariance_matrix.

gaussian_systems.conventions.rotation_matrix(theta)

Construct a single-mode phase-space rotation matrix.

This function returns the 2×2 rotation matrix acting on a single mode in phase space. In the (x, p) quadrature basis, the transformation is

(x, p) → R(theta) (x, p),

where
R(theta) = [[cos(theta), -sin(theta)],

[sin(theta), cos(theta)]].

Parameters:

theta (Real) – The rotation angle in radians. Must be finite.

Returns:

The rotation matrix.

Return type:

numpy.ndarray of shape (2, 2)

Raises:

  • TypeError – If theta is not a real scalar.

  • ValueError – If theta is not finite.

Notes

This matrix represents a symplectic transformation for a single mode in the x-then-p convention.

gaussian_systems.conventions.symmetrize_matrix(matrix)

Return the symmetric (Hermitian) part of a matrix.

This function maps a square matrix matrix to its self-adjoint part: (matrix + matrix.conj().T) / 2.

For real-valued matrices, this reduces to standard symmetrization: (A + A.T) / 2. For complex-valued matrices, this produces a Hermitian matrix.

Parameters:

matrix (numpy.ndarray) – The input matrix. Must be square and contain only finite values.

Returns:

The symmetric (Hermitian) part of the input matrix.

Return type:

numpy.ndarray

Raises:

  • TypeError – If matrix is not a NumPy array.

  • ValueError – If matrix is not square or contains non-finite values.

Notes

The output is always self-adjoint (Hermitian).

gaussian_systems.conventions.symplectic_matrix(n)

Construct the canonical symplectic form for an n-mode system.

This function returns the 2n × 2n symplectic matrix Ω defined by

Ω = [[0, I],

[-I, 0]],

where I is the n × n identity matrix.

In the x-then-p phase-space ordering (x_1, …, x_n, p_1, …, p_n), this matrix defines the canonical commutation relations:

[R_i, R_j] = i Ω_{ij},

where R = (x_1, …, x_n, p_1, …, p_n).

Parameters:

n (Integral) – The number of modes. Must be strictly positive.

Returns:

The symplectic form Ω.

Return type:

numpy.ndarray of shape (2n, 2n)

Raises:

  • TypeError – If n is not an integer.

  • ValueError – If n is not strictly positive.

Notes

The matrix is constructed as kron([[0, 1], [-1, 0]], I_n).