From 308d723b36a46af1c72b451c839d66f0b0f24fad Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 16 Apr 2026 09:59:39 -0600 Subject: [PATCH 01/23] add helper methods for spin & mot. DOFs to basis, add W method, state --- src/ionsim/basis.py | 20 +++++++++++++++++--- src/ionsim/state.py | 32 ++++++++++++++++++++++++++++---- 2 files changed, 45 insertions(+), 7 deletions(-) diff --git a/src/ionsim/basis.py b/src/ionsim/basis.py index b042f4b..88adbb3 100644 --- a/src/ionsim/basis.py +++ b/src/ionsim/basis.py @@ -192,6 +192,9 @@ def spin_DOFs(self): spins = [DOF for DOF in self.degrees_of_freedom if isinstance(DOF, AtomicSpin)] return spins + @property + def motional_modes(self): + return [dof for dof in self.degrees_of_freedom if dof not in self.spin_DOFs] @dataclass(frozen=True, eq=False) class ZPauliBasis(StandardBasis): @@ -221,6 +224,11 @@ def vectors(self): pairs = list(itertools.product(*[[plus, minus] for dof in self.degrees_of_freedom])) return [np.kron(*pair) for pair in pairs] + @property + def spin_DOFs(self): + """ Returns list of spin degrees of freedom """ + return self.degrees_of_freedom + @dataclass(frozen=True, eq=False) class YPauliBasis(Basis): """A basis in which the basis vectors correspond to the (plus/minus) eigenstates of the x-Pauli spin matrix.""" @@ -243,11 +251,16 @@ def vectors(self): @dataclass(frozen=True, eq=False) class XPauliAndFockBasis(Basis): """A basis in which the basis vectors correspond to the (plus/minus) eigenstates of the x-Pauli spin matrix and Fock states.""" - atomic_spins: list[AtomicSpin] + + @property + def spin_DOFs(self): + """ Returns list of spin degrees of freedom or empty list if none. """ + spins = [DOF for DOF in self.degrees_of_freedom if isinstance(DOF, AtomicSpin)] + return spins @property def motional_modes(self): - return [dof for dof in self.degrees_of_freedom if dof not in self.atomic_spins] + return [dof for dof in self.degrees_of_freedom if dof not in self.spin_DOFs] @property def vectors(self): @@ -256,9 +269,10 @@ def vectors(self): minus = 1/np.sqrt(2)*np.array([1, -1]) groups = list(itertools.product( *[ - [plus, minus] if dof in self.atomic_spins else + [plus, minus] if dof in self.spin_DOFs else [np.eye(len(dof.energy_levels))[i] for i in range(len(dof.energy_levels))] for dof in self.degrees_of_freedom ] )) return [ft.reduce(np.kron, group) for group in groups] + diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 0771d0a..3c2573f 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -4,6 +4,7 @@ from ionsim.ionsim_error import IonSimError from ionsim.hamiltonian import Hamiltonian from ionsim.dissipator import Dissipator, Lindbladian +from ionsim.named_operators import Fock import numpy as np # from typing import Any @@ -159,7 +160,7 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio """Compute the coherent displacement (expectation value of the lowering operator) for each spin state.""" assert(len(self.basis.degrees_of_freedom) == len(spin_dofs) + 1) # TODO: trace out other degrees of freedom spin_basis = StandardBasis(spin_dofs) - lowering = lowering_motion(len(motional_dof.energy_levels)) + lowering = Fock.lowering(len(motional_dof.energy_levels)) diplacements = [] for vector in spin_basis.vectors: spin_proj = spin_basis.compute_projector_matrix(vector) @@ -167,7 +168,30 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio diplacements.append(displacement) return diplacements - # def transform_to_spin_eigenbasis(self): + def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): + """ Computes W(x,p) the Wigner distribution for each motional mode in the basis; assumes a Fock basis for each mode. + - requires a specification of the x and p grids as 1-dimensional arrays + """ + from qutip import Qobj, wigner + wigner_distributions = [] + + # Retrieve and trace out spin DOFs from the density matrix + spin_DOFs = self.basis.spin_DOFs + + motional_state = self + for spin in spin_DOFs: + motional_state = motional_state.trace_out_degree_of_freedom(spin) -def lowering_motion(fock_dimension: int): - return np.diag([np.sqrt(n+1) for n in range(fock_dimension-1)], k=1) + # For each mode, trace out all other modes + for i, mode_i in enumerate(self.basis.motional_modes): + mode_state = motional_state # reset the state + for j, mode_j in enumerate(self.basis.motional_modes): + if i == j: + continue + mode_state = motional_state.trace_out_degree_of_freedom(mode_j) + N_fock = len(mode_i.energy_levels) + wigner_distributions.append(wigner(Qobj(mode_state.density_matrix, dims=[[N_fock], [N_fock]]), x_grid, p_grid)) + + return wigner_distributions + + # def transform_to_spin_eigenbasis(self): From 94504d9ce807685818aab67f68f18944034ea311 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 16 Apr 2026 10:36:25 -0600 Subject: [PATCH 02/23] add test for Wigner distribution calcs --- tests/test_state.py | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/tests/test_state.py b/tests/test_state.py index 468065f..4422e6c 100644 --- a/tests/test_state.py +++ b/tests/test_state.py @@ -3,7 +3,7 @@ import numpy as np from ionsim.state import State -from ionsim.degree_of_freedom import AtomicSpin +from ionsim.degree_of_freedom import AtomicSpin, MotionalMode from ionsim.basis import StandardBasis, XPauliBasis from ionsim.named_operators import Pauli from ionsim.testing import assert_array_close @@ -19,6 +19,13 @@ def setUp(self): self.eig_x = np.array([[1, 0], [0, -1]]) self.state = State.from_density_matrix(self.basis_x, np.kron(self.eig_x, Pauli.I)) + # Set up spin-motional state with 2 spins and 1 motional mode + self.mode = MotionalMode.from_frequency(frequency=2*np.pi * 3e6, fock_dimension=5) + self.spin_motion_basis = StandardBasis([self.spin_a, self.spin_b, self.mode]) + state_coefficients = np.zeros(len(self.spin_motion_basis.states)) + state_coefficients[0] = 1. + self.spin_motion_state = State.from_coefficients(self.spin_motion_basis, list(state_coefficients)) + def test_density_matrix_in_new_basis(self): """Test the density matrix in the new basis.""" rho_p = State.from_state(self.basis, self.state).density_matrix @@ -27,5 +34,27 @@ def test_density_matrix_in_new_basis(self): # Check that the density matrix matches the expected value assert_array_close(rho_p, expected_rho_p) + def test_wigner_function(self): + """ Test Wigner function computation from a motional density matrix in the Fock state basis""" + domain_limit = 2. + x_grid = np.linspace(-domain_limit, domain_limit, 25) + p_grid = np.linspace(-domain_limit, domain_limit, 25) + + W_distribution = self.spin_motion_state.compute_wigner_distribution(x_grid, p_grid)[0] + + W_max = np.max(W_distribution) + # Test just one slice of the array and half of its contents (U(1) symmetric about origin)) + slice_indx = 12 + W_expected_at_slice = np.array([0.00583005, 0.0110443, 0.0197914, 0.03354962, 0.05379861, 0.08160694, 0.11709966, 0.15894861, 0.20409406, 0.24789999, 0.2848362, 0.30958962, 0.31830989]) + + norm = np.trapezoid(np.trapezoid(W_distribution, x_grid, axis=0), p_grid) + expected_norm = 0.9902872798196112 + self.assertAlmostEqual(norm, expected_norm, places=7) + expected_W_max = 0.31830988618379075 + self.assertAlmostEqual(W_max, expected_W_max, places=7) + + assert_array_close(W_distribution[slice_indx, 0:13], W_expected_at_slice, rtol=1e-04, atol=1E-7) + + if __name__ == '__main__': unittest.main() From 150d8faac76346e729310077925609e896c42131 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 16 Apr 2026 10:37:18 -0600 Subject: [PATCH 03/23] finalize branch for MR --- tests/test_state.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/tests/test_state.py b/tests/test_state.py index 4422e6c..260a2c0 100644 --- a/tests/test_state.py +++ b/tests/test_state.py @@ -39,20 +39,17 @@ def test_wigner_function(self): domain_limit = 2. x_grid = np.linspace(-domain_limit, domain_limit, 25) p_grid = np.linspace(-domain_limit, domain_limit, 25) - W_distribution = self.spin_motion_state.compute_wigner_distribution(x_grid, p_grid)[0] - W_max = np.max(W_distribution) # Test just one slice of the array and half of its contents (U(1) symmetric about origin)) slice_indx = 12 W_expected_at_slice = np.array([0.00583005, 0.0110443, 0.0197914, 0.03354962, 0.05379861, 0.08160694, 0.11709966, 0.15894861, 0.20409406, 0.24789999, 0.2848362, 0.30958962, 0.31830989]) - norm = np.trapezoid(np.trapezoid(W_distribution, x_grid, axis=0), p_grid) expected_norm = 0.9902872798196112 self.assertAlmostEqual(norm, expected_norm, places=7) expected_W_max = 0.31830988618379075 + W_max = np.max(W_distribution) self.assertAlmostEqual(W_max, expected_W_max, places=7) - assert_array_close(W_distribution[slice_indx, 0:13], W_expected_at_slice, rtol=1e-04, atol=1E-7) From fbd2844cf1fbdd8c23f69f9f21e399a9de6016df Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 16 Apr 2026 10:44:28 -0600 Subject: [PATCH 04/23] + check that state is in Fock number state basis in Wigner calc --- src/ionsim/state.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 3c2573f..7b5648b 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -5,6 +5,7 @@ from ionsim.hamiltonian import Hamiltonian from ionsim.dissipator import Dissipator, Lindbladian from ionsim.named_operators import Fock +from ionsim.collective_motional_energy_level import CollectiveMotionalEnergyLevel import numpy as np # from typing import Any @@ -170,7 +171,10 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): """ Computes W(x,p) the Wigner distribution for each motional mode in the basis; assumes a Fock basis for each mode. - - requires a specification of the x and p grids as 1-dimensional arrays + - requires a specification of the x and p grids as 1-dimensional arrays, determines resolution of Wigner distributions. + - returns a list of Wigner distributions -> one for each mode. + - returns an empty list if there are no modes in the basis. + - assumes the motional modes are in the Fock number state basis """ from qutip import Qobj, wigner wigner_distributions = [] @@ -184,6 +188,9 @@ def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): # For each mode, trace out all other modes for i, mode_i in enumerate(self.basis.motional_modes): + if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): + raise IonSimError("Wigner distribution calculation assumes motional mode is in the Fock number state basis.") + mode_state = motional_state # reset the state for j, mode_j in enumerate(self.basis.motional_modes): if i == j: From 5bd246486b7b537f26ecf79c6b725b439983690e Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 16 Apr 2026 10:47:44 -0600 Subject: [PATCH 05/23] add qutip dependency to pyproject toml req's --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 1c0d792..a0fca57 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,7 +17,8 @@ dependencies = [ "icecream", "pyyaml", "pint", - "sympy" + "sympy", + "qutip" ] classifiers = [ "Development Status :: 2 - Beta", From fcb175e467663dd12045b8f0ab8234cb88ad4a19 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Tue, 21 Apr 2026 12:19:00 -0600 Subject: [PATCH 06/23] add type hinting to Wigner dist fxn and check CI/CD --- src/ionsim/state.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 7b5648b..778b8df 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -169,7 +169,7 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio diplacements.append(displacement) return diplacements - def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): + def compute_wigner_distribution(self, x_grid: Vector[float], p_grid: Vector[float]): """ Computes W(x,p) the Wigner distribution for each motional mode in the basis; assumes a Fock basis for each mode. - requires a specification of the x and p grids as 1-dimensional arrays, determines resolution of Wigner distributions. - returns a list of Wigner distributions -> one for each mode. From 4370549bbd2be4ad26c43d99edb1bc65be9482b7 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Tue, 21 Apr 2026 12:23:44 -0600 Subject: [PATCH 07/23] rm type hinting with Vector[float], functionality not supported --- src/ionsim/state.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 778b8df..7b5648b 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -169,7 +169,7 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio diplacements.append(displacement) return diplacements - def compute_wigner_distribution(self, x_grid: Vector[float], p_grid: Vector[float]): + def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): """ Computes W(x,p) the Wigner distribution for each motional mode in the basis; assumes a Fock basis for each mode. - requires a specification of the x and p grids as 1-dimensional arrays, determines resolution of Wigner distributions. - returns a list of Wigner distributions -> one for each mode. From fde0f719bbacd7dbc36336c127942140a89c7782 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Mon, 18 May 2026 22:08:47 -0600 Subject: [PATCH 08/23] rm unused lists from dissipator test --- tests/test_dissipators.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/test_dissipators.py b/tests/test_dissipators.py index 24592f1..778f972 100644 --- a/tests/test_dissipators.py +++ b/tests/test_dissipators.py @@ -133,10 +133,6 @@ def setUp(self): spin_identity = np.sum(spin_identities, axis=0) - # Tensor raising/lowering operators with spin identity for each mode that is heating - lowering_operators = [] - raising_operators = [] - # Basis is of the form spins x modes, so we tensor in the that order: N_thermal = case['N thermal'] lowering_operator = np.sqrt(decay_rate*(N_thermal+1)) * skron(spin_identity, mode_lowering_matrix) From ebe71da8c28fd3dc24c5aa625b7d68320c4c7e39 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 09:55:49 -0600 Subject: [PATCH 09/23] add X and P operators to Fock named operators class --- src/ionsim/named_operators.py | 8 +++++ src/ionsim/state.py | 60 +++++++++++++++++++++++++++++++---- 2 files changed, 62 insertions(+), 6 deletions(-) diff --git a/src/ionsim/named_operators.py b/src/ionsim/named_operators.py index 16d012e..c054768 100644 --- a/src/ionsim/named_operators.py +++ b/src/ionsim/named_operators.py @@ -49,6 +49,14 @@ def raising(cls, fock_dimension: int): def number(cls, fock_dimension: int): return np.diag(np.arange(fock_dimension)) + @classmethod # or "x" or "X"? + def position(cls, fock_dimension: int): + return np.sqrt(0.5)*(cls.raising(fock_dimension) + cls.lowering(fock_dimension)) + + @classmethod # or "p" or "P"? + def momentum(cls, fock_dimension: int): + return 1j*np.sqrt(0.5)*(cls.raising(fock_dimension) - cls.lowering(fock_dimension)) + @staticmethod def debye_waller_lowering(fock_dimension: int, lamb_dicke_parameter: float): """The lowering operator for a harmonic oscillator.""" diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 7b5648b..9fdfd56 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -169,6 +169,54 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio diplacements.append(displacement) return diplacements + + @property + def motional_state(self): + """ Motional portion of the state, if it exists. Returns None if no motional DOFs """ + if not self.basis.motional_modes: + # Raise an error if there are no motional modes in the basis + raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") + + # Trace out spin DOFs to obtain a purely motional state + motional_state = self + for spin in self.basis.spin_DOFs: + motional_state = motional_state.trace_out_degree_of_freedom(spin) + + return motional_state + + def compute_quadratures(self): + """ Returns quadrature expectation values and

for each motional mode in the state. + - fails if there are no motional DOFs + - returns lists of and

; each list element corresponds to 1 motional mode + - list order matches DOF order in the state's basis + """ + + # if not self.basis.motional_modes: + # # Raise an error if there are no motional modes in the basis + # raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") + # + # # Trace out spin DOFs if present + # motional_state = self + # for spin in self.basis.spin_DOFs: + # motional_state = motional_state.trace_out_degree_of_freedom(spin) + + x = [] + p = [] + for i, mode_i in enumerate(self.basis.motional_modes): + if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): + raise IonSimError("Quadrature calculation assumes motional mode is in the Fock number state basis.") + + mode_state = self.motional_state # reset the state + # Trace out other modes + for j, mode_j in enumerate(self.basis.motional_modes): + if i == j: + continue + mode_state = motional_state.trace_out_degree_of_freedom(mode_j) + + + + + def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): """ Computes W(x,p) the Wigner distribution for each motional mode in the basis; assumes a Fock basis for each mode. - requires a specification of the x and p grids as 1-dimensional arrays, determines resolution of Wigner distributions. @@ -180,18 +228,18 @@ def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): wigner_distributions = [] # Retrieve and trace out spin DOFs from the density matrix - spin_DOFs = self.basis.spin_DOFs - - motional_state = self - for spin in spin_DOFs: - motional_state = motional_state.trace_out_degree_of_freedom(spin) + # spin_DOFs = self.basis.spin_DOFs + # + # motional_state = self + # for spin in spin_DOFs: + # motional_state = motional_state.trace_out_degree_of_freedom(spin) # For each mode, trace out all other modes for i, mode_i in enumerate(self.basis.motional_modes): if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): raise IonSimError("Wigner distribution calculation assumes motional mode is in the Fock number state basis.") - mode_state = motional_state # reset the state + mode_state = self.motional_state # reset the state for j, mode_j in enumerate(self.basis.motional_modes): if i == j: continue From 5d72246ed309322a40b7249993c0ee6298ffc909 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 10:14:39 -0600 Subject: [PATCH 10/23] add observable computational methods with density matrix, add

--- src/ionsim/state.py | 82 +++++++++++++++++++++++++++++++-------------- 1 file changed, 56 insertions(+), 26 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 9fdfd56..536de7d 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -25,6 +25,20 @@ class State: def supervector(self) -> Vector: return self.basis.compute_supervector_from_density_matrix(self.density_matrix) + @property + def motional_state(self): + """ Motional portion of the state, if it exists. Returns None if no motional DOFs """ + if not self.basis.motional_modes: + # Raise an error if there are no motional modes in the basis + raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") + + # Trace out spin DOFs to obtain a purely motional state + motional_state = self + for spin in self.basis.spin_DOFs: + motional_state = motional_state.trace_out_degree_of_freedom(spin) + + return motional_state + @classmethod def from_coefficients(cls, basis: Basis, coefficients: list[float]): """Build a state from a list of basis-state coefficients.""" @@ -169,39 +183,46 @@ def compute_coherent_displacements(self, spin_dofs: list[DegreeOfFreedom], motio diplacements.append(displacement) return diplacements + ### Added methods by ECM in this branch:### + def compute_matrix_observable_expectation(self, observable_operator: Matrix) -> Matrix: + """ Compute the expectation value of an operator observable, represented by . - @property - def motional_state(self): - """ Motional portion of the state, if it exists. Returns None if no motional DOFs """ - if not self.basis.motional_modes: - # Raise an error if there are no motional modes in the basis - raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") + - Computes via Tr[O rho], where rho is the density matrix. + - O and rho must be in compatable bases. - # Trace out spin DOFs to obtain a purely motional state - motional_state = self - for spin in self.basis.spin_DOFs: - motional_state = motional_state.trace_out_degree_of_freedom(spin) - - return motional_state + """ + # Check that observable and density matrix are compatible + basis_size = len(self.basis.vectors) + if observable_operator.shape != (basis_size, basis_size): + raise IonSimError("Observable operator must correspond with a matrix of shape: {(basis_size, basis_size)}.") + + return np.trace(observable_operator.dot(self.density_matrix)) + + def compute_operator_observable_expectation(self, observable_operator: Operator) -> Matrix: + """ Compute the expectation value of an operator observable, represented by . + + - Computes via Tr[O rho], where rho is the density matrix. + - O and rho must be in compatable bases. + + """ + # Check that observable and density matrix are compatible + if isinstance(observable_operator, CouplingOperator) and observable_operator.oscillation_rate != 0.): + raise IonSimeError(f"Observable operator must be a static operator, but oscillation rate is {observable_operator.oscillation_rate}.") - def compute_quadratures(self): + return self.compute_matrix_observable_expectation(observable_operator.static_matrix) + + def compute_quadratures(self, include_variance: bool=False): """ Returns quadrature expectation values and

for each motional mode in the state. - fails if there are no motional DOFs - returns lists of and

; each list element corresponds to 1 motional mode - list order matches DOF order in the state's basis """ - - # if not self.basis.motional_modes: - # # Raise an error if there are no motional modes in the basis - # raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") - # - # # Trace out spin DOFs if present - # motional_state = self - # for spin in self.basis.spin_DOFs: - # motional_state = motional_state.trace_out_degree_of_freedom(spin) - x = [] p = [] + if include_variance: + x2 = [] + p2 = [] + for i, mode_i in enumerate(self.basis.motional_modes): if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): raise IonSimError("Quadrature calculation assumes motional mode is in the Fock number state basis.") @@ -212,10 +233,19 @@ def compute_quadratures(self): if i == j: continue mode_state = motional_state.trace_out_degree_of_freedom(mode_j) - - - + Fock_dim = len(mode_i.energy_levels) + # Compute expectation values: + x[i] = mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim)) + p[i] = mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim)) + if include_variance: + x2[i] = mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim)) + p2[i] = mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim)) + if include_variance: + return x, p, x2, p2 + else: + return x, p + def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): """ Computes W(x,p) the Wigner distribution for each motional mode in the basis; assumes a Fock basis for each mode. From 62cfe470d85d6f89d9eb44ccf84249e2da3f4340 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 10:44:25 -0600 Subject: [PATCH 11/23] test x,p and wigner distribution methods + factoring --- src/ionsim/state.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 536de7d..592934b 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -33,11 +33,11 @@ def motional_state(self): raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") # Trace out spin DOFs to obtain a purely motional state - motional_state = self + state = self for spin in self.basis.spin_DOFs: - motional_state = motional_state.trace_out_degree_of_freedom(spin) + state = state.trace_out_degree_of_freedom(spin) - return motional_state + return state @classmethod def from_coefficients(cls, basis: Basis, coefficients: list[float]): @@ -206,7 +206,7 @@ def compute_operator_observable_expectation(self, observable_operator: Operator) """ # Check that observable and density matrix are compatible - if isinstance(observable_operator, CouplingOperator) and observable_operator.oscillation_rate != 0.): + if isinstance(observable_operator, CouplingOperator) and observable_operator.oscillation_rate != 0.: raise IonSimeError(f"Observable operator must be a static operator, but oscillation rate is {observable_operator.oscillation_rate}.") return self.compute_matrix_observable_expectation(observable_operator.static_matrix) @@ -232,14 +232,14 @@ def compute_quadratures(self, include_variance: bool=False): for j, mode_j in enumerate(self.basis.motional_modes): if i == j: continue - mode_state = motional_state.trace_out_degree_of_freedom(mode_j) + mode_state = mode_state.trace_out_degree_of_freedom(mode_j) Fock_dim = len(mode_i.energy_levels) # Compute expectation values: - x[i] = mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim)) - p[i] = mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim)) + x.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim))) + p.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim))) if include_variance: - x2[i] = mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim)) - p2[i] = mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim)) + x2.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim))) + p2.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim))) if include_variance: return x, p, x2, p2 @@ -273,7 +273,7 @@ def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): for j, mode_j in enumerate(self.basis.motional_modes): if i == j: continue - mode_state = motional_state.trace_out_degree_of_freedom(mode_j) + mode_state = mode_state.trace_out_degree_of_freedom(mode_j) N_fock = len(mode_i.energy_levels) wigner_distributions.append(wigner(Qobj(mode_state.density_matrix, dims=[[N_fock], [N_fock]]), x_grid, p_grid)) From c929f8ecad5e9c27b41826dd93b0ff6379c4c14f Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 10:56:14 -0600 Subject: [PATCH 12/23] add Operator import --- src/ionsim/state.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 592934b..8c0a9c3 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -1,4 +1,5 @@ from ionsim.basis import Basis, StandardBasis +from ionsim.operator import Operator, CouplingOperator from ionsim.degree_of_freedom import DegreeOfFreedom from ionsim.custom_types import Vector, Matrix from ionsim.ionsim_error import IonSimError From 97a969e022340b73248c67a8f6ea8271c36e0c77 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 11:24:03 -0600 Subject: [PATCH 13/23] add tests for quadrature computation --- tests/test_state.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/test_state.py b/tests/test_state.py index 260a2c0..94227cb 100644 --- a/tests/test_state.py +++ b/tests/test_state.py @@ -53,5 +53,19 @@ def test_wigner_function(self): assert_array_close(W_distribution[slice_indx, 0:13], W_expected_at_slice, rtol=1e-04, atol=1E-7) + def test_quadrature_computation(self): + """ Yest computation of x, p, and x^2, and p^2 expectations """ + include_variance = True + x, p, x2, p2 = self.spin_motion_state.compute_quadratures(include_variance) + + x_expected = 0. + p_expected = 0. + x2_expected = 0.5 + p2_expected = 0.5 + self.assertAlmostEqual(x[0].real, x_expected, places=6) + self.assertAlmostEqual(p[0].real, p_expected, places=6) + self.assertAlmostEqual(x2[0].real, x2_expected, places=6) + self.assertAlmostEqual(p2[0].real, p2_expected, places=6) + if __name__ == '__main__': unittest.main() From 074ba5a983d4743d3bbc59f547d56743745ec4cf Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 11:25:19 -0600 Subject: [PATCH 14/23] rm commented code that was factored --- src/ionsim/state.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 8c0a9c3..715195d 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -258,13 +258,6 @@ def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): from qutip import Qobj, wigner wigner_distributions = [] - # Retrieve and trace out spin DOFs from the density matrix - # spin_DOFs = self.basis.spin_DOFs - # - # motional_state = self - # for spin in spin_DOFs: - # motional_state = motional_state.trace_out_degree_of_freedom(spin) - # For each mode, trace out all other modes for i, mode_i in enumerate(self.basis.motional_modes): if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): From bd6693b97a7bc95ff6ccfa7ad28c9741b6de8d2a Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 11:49:13 -0600 Subject: [PATCH 15/23] consolidate checking for motional modes in basis --- src/ionsim/state.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 715195d..2b2a5ff 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -9,7 +9,6 @@ from ionsim.collective_motional_energy_level import CollectiveMotionalEnergyLevel import numpy as np -# from typing import Any from dataclasses import dataclass from numpy.linalg import multi_dot @@ -30,8 +29,7 @@ def supervector(self) -> Vector: def motional_state(self): """ Motional portion of the state, if it exists. Returns None if no motional DOFs """ if not self.basis.motional_modes: - # Raise an error if there are no motional modes in the basis - raise IonSimError("No motional modes present in the basis. Cannot compute quadratures.") + return None # Trace out spin DOFs to obtain a purely motional state state = self @@ -234,6 +232,7 @@ def compute_quadratures(self, include_variance: bool=False): if i == j: continue mode_state = mode_state.trace_out_degree_of_freedom(mode_j) + assert mode_state is not None Fock_dim = len(mode_i.energy_levels) # Compute expectation values: x.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim))) @@ -264,6 +263,8 @@ def compute_wigner_distribution(self, x_grid: Vector, p_grid: Vector): raise IonSimError("Wigner distribution calculation assumes motional mode is in the Fock number state basis.") mode_state = self.motional_state # reset the state + assert mode_state is not None + for j, mode_j in enumerate(self.basis.motional_modes): if i == j: continue From 020d9b2603a7236a853fc2cb6030ffc4bb2b4d19 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 11:57:09 -0600 Subject: [PATCH 16/23] rename variable name for concept accuracy --- src/ionsim/state.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 2b2a5ff..7e8080f 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -210,15 +210,16 @@ def compute_operator_observable_expectation(self, observable_operator: Operator) return self.compute_matrix_observable_expectation(observable_operator.static_matrix) - def compute_quadratures(self, include_variance: bool=False): + def compute_quadratures(self, enable_squared_quadratures: bool=False): """ Returns quadrature expectation values and

for each motional mode in the state. - fails if there are no motional DOFs - returns lists of and

; each list element corresponds to 1 motional mode - list order matches DOF order in the state's basis + - option to return and , the expectation of the squared quadrature operators. """ x = [] p = [] - if include_variance: + if enable_squared_quadratures: x2 = [] p2 = [] @@ -237,11 +238,11 @@ def compute_quadratures(self, include_variance: bool=False): # Compute expectation values: x.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim))) p.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim))) - if include_variance: + if enable_squared_quadratures: x2.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim))) p2.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim))) - if include_variance: + if enable_squared_quadratures: return x, p, x2, p2 else: return x, p From c968855d53313f6eae1cf3c7c742e4c65d9eeffe Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 11:57:45 -0600 Subject: [PATCH 17/23] update test --- tests/test_state.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_state.py b/tests/test_state.py index 94227cb..c321fc1 100644 --- a/tests/test_state.py +++ b/tests/test_state.py @@ -55,8 +55,8 @@ def test_wigner_function(self): def test_quadrature_computation(self): """ Yest computation of x, p, and x^2, and p^2 expectations """ - include_variance = True - x, p, x2, p2 = self.spin_motion_state.compute_quadratures(include_variance) + include_squared_quantities = True + x, p, x2, p2 = self.spin_motion_state.compute_quadratures(include_squared_quantities) x_expected = 0. p_expected = 0. From c1b58285c111b65ef3fb4c5eb7a65b23335b627f Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 20 May 2026 18:01:53 -0600 Subject: [PATCH 18/23] add note to look into computation times for partial trace --- src/ionsim/state.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 7e8080f..8643a9c 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -223,6 +223,7 @@ def compute_quadratures(self, enable_squared_quadratures: bool=False): x2 = [] p2 = [] + # TODO: Check this: Wes says its cheaper to compute observable in the full space; for i, mode_i in enumerate(self.basis.motional_modes): if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): raise IonSimError("Quadrature calculation assumes motional mode is in the Fock number state basis.") From 3b0bce27bcee1c87db98aa71dc169d4825f63413 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Fri, 22 May 2026 10:55:32 -0600 Subject: [PATCH 19/23] add alternative calculation for x,p in full space --- src/ionsim/state.py | 70 +++++++++++++++++++++++++++++++++------------ 1 file changed, 51 insertions(+), 19 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 8643a9c..160ed9f 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -223,26 +223,58 @@ def compute_quadratures(self, enable_squared_quadratures: bool=False): x2 = [] p2 = [] - # TODO: Check this: Wes says its cheaper to compute observable in the full space; - for i, mode_i in enumerate(self.basis.motional_modes): - if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): - raise IonSimError("Quadrature calculation assumes motional mode is in the Fock number state basis.") - - mode_state = self.motional_state # reset the state - # Trace out other modes - for j, mode_j in enumerate(self.basis.motional_modes): - if i == j: - continue - mode_state = mode_state.trace_out_degree_of_freedom(mode_j) - assert mode_state is not None - Fock_dim = len(mode_i.energy_levels) - # Compute expectation values: - x.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim))) - p.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim))) - if enable_squared_quadratures: - x2.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim))) - p2.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim))) + use_partial_trace = False + #import time + #start = time.perf_counter() + # Either compute partial trace on the density matrix or elevate observable to full space + # It may be cheaper to compute the expectation value in the full space + if not use_partial_trace : + full_basis = self.basis + + # For each mode, build the full-space raising and lowering operators + for i, mode_i in enumerate(self.basis.motional_modes): + if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): + raise IonSimError("Quadrature calculation assumes motional mode is in the Fock number state basis.") + + # Enlarge the x and p operators for the current mode + Fock_dim = len(mode_i.energy_levels) + x_enlarged = full_basis.enlarge_one_dof_matrix(Fock.position(Fock_dim), mode_i) + p_enlarged = full_basis.enlarge_one_dof_matrix(Fock.momentum(Fock_dim), mode_i) + + # Compute expectation values: + x.append(self.compute_matrix_observable_expectation(x_enlarged)) + p.append(self.compute_matrix_observable_expectation(p_enlarged)) + if enable_squared_quadratures: + x2.append(self.compute_matrix_observable_expectation(x_enlarged @ x_enlarged)) + p2.append(self.compute_matrix_observable_expectation(p_enlarged @ p_enlarged)) + + else: + for i, mode_i in enumerate(self.basis.motional_modes): + if not isinstance(mode_i.energy_levels[0], CollectiveMotionalEnergyLevel): + raise IonSimError("Quadrature calculation assumes motional mode is in the Fock number state basis.") + + mode_state = self.motional_state # reset the state + # Trace out other modes + for j, mode_j in enumerate(self.basis.motional_modes): + if i == j: + continue + mode_state = mode_state.trace_out_degree_of_freedom(mode_j) + assert mode_state is not None + Fock_dim = len(mode_i.energy_levels) + # Compute expectation values: + x.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim))) + p.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim))) + if enable_squared_quadratures: + x2.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim))) + p2.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim))) + + #end = time.perf_counter() + # if not use_partial_trace : + # print(f"Promoting observable to full space") + # else: + # print(f"Using partial traces") + # print(f"Computing x,p took: {end-start} [s]") if enable_squared_quadratures: return x, p, x2, p2 else: From d3a683281ca0eedfcda9ec9fadd29d4380c9c102 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Tue, 26 May 2026 14:45:42 -0600 Subject: [PATCH 20/23] rm comments for x,p; wigner calculations --- src/ionsim/state.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/ionsim/state.py b/src/ionsim/state.py index 160ed9f..949076e 100644 --- a/src/ionsim/state.py +++ b/src/ionsim/state.py @@ -225,8 +225,6 @@ def compute_quadratures(self, enable_squared_quadratures: bool=False): use_partial_trace = False - #import time - #start = time.perf_counter() # Either compute partial trace on the density matrix or elevate observable to full space # It may be cheaper to compute the expectation value in the full space if not use_partial_trace : @@ -269,12 +267,6 @@ def compute_quadratures(self, enable_squared_quadratures: bool=False): x2.append(mode_state.compute_matrix_observable_expectation(Fock.position(Fock_dim) @ Fock.position(Fock_dim))) p2.append(mode_state.compute_matrix_observable_expectation(Fock.momentum(Fock_dim) @ Fock.momentum(Fock_dim))) - #end = time.perf_counter() - # if not use_partial_trace : - # print(f"Promoting observable to full space") - # else: - # print(f"Using partial traces") - # print(f"Computing x,p took: {end-start} [s]") if enable_squared_quadratures: return x, p, x2, p2 else: From f962165c954e07586e87fa3b404b7e13523fad22 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Wed, 1 Jul 2026 11:56:35 -0600 Subject: [PATCH 21/23] update MS gate example to use basis property method --- examples/example_simulated_multimode_MS_gate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/example_simulated_multimode_MS_gate.py b/examples/example_simulated_multimode_MS_gate.py index b898c6e..a125e50 100644 --- a/examples/example_simulated_multimode_MS_gate.py +++ b/examples/example_simulated_multimode_MS_gate.py @@ -163,7 +163,7 @@ def main(): end = time.perf_counter() ic(f'Propagating state took {end - start} s.') - basis_xx = ism.XPauliAndFockBasis([*spins, *modes], spins) + basis_xx = ism.XPauliAndFockBasis([*spins, *modes]) psis_xx = [ism.State.from_state(basis_xx, psi) for psi in psis] alphas = np.array([psi.compute_coherent_displacements(spins, modes[0]) for psi in psis_xx]) @@ -218,4 +218,4 @@ def main(): if __name__ == '__main__': - main() \ No newline at end of file + main() From ab382b5c8c1bf19afe975163f99161327c4942cd Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 23 Jul 2026 11:53:37 -0600 Subject: [PATCH 22/23] fix MS example and rm merge conflict comments --- examples/example_simulated_multimode_MS_gate.py | 2 +- src/ionsim/basis.py | 17 ----------------- 2 files changed, 1 insertion(+), 18 deletions(-) diff --git a/examples/example_simulated_multimode_MS_gate.py b/examples/example_simulated_multimode_MS_gate.py index a1e0c11..805638f 100644 --- a/examples/example_simulated_multimode_MS_gate.py +++ b/examples/example_simulated_multimode_MS_gate.py @@ -163,7 +163,7 @@ def main(): end = time.perf_counter() ic(f'Propagating state took {end - start} s.') - basis_xx = ism.XPauliAndFockBasis([*spins, *modes]) + basis_xx = ism.XPauliAndFockBasis([*spins, *modes], spins) psis_xx = [ism.State.from_state(basis_xx, psi) for psi in psis] alphas = np.array([psi.compute_coherent_displacements(spins, modes[0]) for psi in psis_xx]) diff --git a/src/ionsim/basis.py b/src/ionsim/basis.py index 52299b0..4363198 100644 --- a/src/ionsim/basis.py +++ b/src/ionsim/basis.py @@ -259,24 +259,11 @@ def vectors(self): @dataclass(frozen=True, eq=False) class XPauliAndFockBasis(Basis): """A basis in which the basis vectors correspond to the (plus/minus) eigenstates of the x-Pauli spin matrix and Fock states.""" -#<<<<<<< HEAD - - # @property - # def spin_DOFs(self): - # """ Returns list of spin degrees of freedom or empty list if none. """ - # spins = [DOF for DOF in self.degrees_of_freedom if isinstance(DOF, AtomicSpin)] - # return spins - - # @property - # def motional_modes(self): - # return [dof for dof in self.degrees_of_freedom if dof not in self.spin_DOFs] -#======= atomic_structure_DOFs: list[AtomicStructure] @property def motional_modes(self): return [dof for dof in self.degrees_of_freedom if dof not in self.atomic_structure_DOFs] -#>>>>>>> main @property def vectors(self): @@ -285,11 +272,7 @@ def vectors(self): minus = 1/np.sqrt(2)*np.array([1, -1]) groups = list(itertools.product( *[ - #<<<<<<< HEAD - # [plus, minus] if dof in self.spin_DOFs else - #======= [plus, minus] if dof in self.atomic_structure_DOFs else -#>>>>>>> main [np.eye(len(dof.energy_levels))[i] for i in range(len(dof.energy_levels))] for dof in self.degrees_of_freedom ] From 577c8aaab754825f1d8493c4e5507016c00e6775 Mon Sep 17 00:00:00 2001 From: Ethan McGarrigle Date: Thu, 23 Jul 2026 12:14:30 -0600 Subject: [PATCH 23/23] fix syntax error in pyproj, test CI/CD --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 81d2d82..00f9fa5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -18,7 +18,7 @@ dependencies = [ "pyyaml", "pint", "sympy", - "qutip" + "qutip", "h5py" ] classifiers = [