Creating Custom Observables

Creating Custom Observables with the Rydberg Basis ("g", "r")¶

This tutorial explains how to create custom observables for measuring expectation values during quantum simulations using emu_mps.

There are TWO main approaches to define operators for the Expectation observable:

1. MPO.from_operator_repr(): Use symbolic notation with basis state labels
2. Direct MPO construction: Build tensors manually with matrix elements

Both methods work with the Rydberg basis (with leakage you add the"x" and XY you can do it with "1" and "0" ) where:

• $|g \rangle$ (ground state) is represented as [1, 0]^T
• $|r \rangle$ (Rydberg state) is represented as [0, 1]^T

The matrix element notation "ab" (e.g., "rg", "gr") corresponds to $|a \rangle \langle b|$.

In [1]:

import torch
from math import pi

from pulser import Register
from pulser.pulse import Pulse
from pulser.sequence import Sequence
from pulser.devices import MockDevice

from emu_mps import MPO, MPSBackend, MPSConfig, Expectation

dtype = torch.complex128  # emu_mps base dtype for MPS and MPO

import torch from math import pi from pulser import Register from pulser.pulse import Pulse from pulser.sequence import Sequence from pulser.devices import MockDevice from emu_mps import MPO, MPSBackend, MPSConfig, Expectation dtype = torch.complex128 # emu_mps base dtype for MPS and MPO

Setup: Create a simple 2-atom constant pulse¶

In [2]:

natoms = 2
reg = Register.rectangle(1, natoms, spacing=10000.0, prefix="q")

seq = Sequence(reg, MockDevice)
seq.declare_channel("ch0", "rydberg_global")
duration = 500
pulse = Pulse.ConstantPulse(duration, 2 * pi, 0.0, 0.0)
seq.add(pulse, "ch0")

natoms = 2 reg = Register.rectangle(1, natoms, spacing=10000.0, prefix="q") seq = Sequence(reg, MockDevice) seq.declare_channel("ch0", "rydberg_global") duration = 500 pulse = Pulse.ConstantPulse(duration, 2 * pi, 0.0, 0.0) seq.add(pulse, "ch0")

Approach 1: Using MPO.from_operator_repr() (Symbolic Notation)¶

This is the recommended approach for most use cases. You define operators using symbolic notation based on the basis states.

The operations parameter has the structure:

operations = [(coefficient, [(operator_dict, target_qubits), ...])]

Where:

• coefficient: a scalar multiplying the entire term
• operator_dict: defines the single-qubit operator using basis labels
• target_qubits: list of qubit indices where the operator acts

The operator_dict uses the notation {"ab": value} meaning: value * |a><b|

For the Rydberg basis ("g", "r"):

• "gg" = $|g \rangle \langle g|$ (projection onto ground state)
• "rr" = $|r \rangle \langle r|$ (projection onto Rydberg state, i.e., n operator)
• "gr" = $|g \rangle \langle r|$ (lowering operator, sigma-)
• "rg" = $|r \rangle \langle g|$ (raising operator, sigma+)

Common operators in terms of these basis elements:

• Pauli X: {"rg": 1.0, "gr": 1.0}
• Pauli Y: {"rg": -1j, "gr": 1j}
• Pauli Z: {"gg": 1.0, "rr": -1.0}
• Number operator (n): {"rr": 1.0}
• Identity: {"gg": 1.0, "rr": 1.0}

Example 1.1: Pauli X on qubit 0¶

$X = |r \rangle \langle g| + |g \rangle \langle r|$

In [3]:

operations_X0 = [(1.0, [({"rg": 1.0, "gr": 1.0}, [0])])]
X0 = MPO.from_operator_repr(
    eigenstates=("g", "r"), n_qudits=natoms, operations=operations_X0
)
print("\n1.1 Pauli X on qubit 0:")
print(f"    operations = {operations_X0}")
print(f"    Shape of first factor: {X0.factors[0].shape}\n")
print(f"    MPO factors: {X0.factors}")

operations_X0 = [(1.0, [({"rg": 1.0, "gr": 1.0}, [0])])] X0 = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=operations_X0 ) print("\n1.1 Pauli X on qubit 0:") print(f" operations = {operations_X0}") print(f" Shape of first factor: {X0.factors[0].shape}\n") print(f" MPO factors: {X0.factors}")

1.1 Pauli X on qubit 0:
    operations = [(1.0, [({'rg': 1.0, 'gr': 1.0}, [0])])]
    Shape of first factor: torch.Size([1, 2, 2, 1])

    MPO factors: [tensor([[[[0.+0.j],
          [1.+0.j]],

         [[1.+0.j],
          [0.+0.j]]]], dtype=torch.complex128), tensor([[[[1.+0.j],
          [0.+0.j]],

         [[0.+0.j],
          [1.+0.j]]]], dtype=torch.complex128)]

Example 1.2: Tensor product $X_0 \otimes X_1$ (with coefficient 2)¶

This measures correlations in the $X$ basis

In [4]:

operations_XX = [(2.0, [({"rg": 1.0, "gr": 1.0}, [0, 1])])]
XX = MPO.from_operator_repr(
    eigenstates=("g", "r"), n_qudits=natoms, operations=operations_XX
)
print(r"1.2 Tensor product 2 * X0 x X1:") 
print(f"    operations = {operations_XX}")

operations_XX = [(2.0, [({"rg": 1.0, "gr": 1.0}, [0, 1])])] XX = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=operations_XX ) print(r"1.2 Tensor product 2 * X0 x X1:") print(f" operations = {operations_XX}")

1.2 Tensor product 2 * X0 x X1:
    operations = [(2.0, [({'rg': 1.0, 'gr': 1.0}, [0, 1])])]

Example 1.3: Number operator on qubit 1¶

$n_1 = |r \rangle \langle r|_1 = Id \otimes |r \rangle \langle r| $, where $Id$ is the identity 2 by 2

In [5]:

operations_n1 = [(1.0, [({"rr": 1.0}, [1])])]
n1 = MPO.from_operator_repr(
    eigenstates=("g", "r"), n_qudits=natoms, operations=operations_n1
)

operations_n1 = [(1.0, [({"rr": 1.0}, [1])])] n1 = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=operations_n1 )

Example 1.4: Sum of local operators: $X_0 + X_1$¶

Each term in the list is a separate operator that gets summed

In [6]:

operations_X0_plus_X1 = [
    (1.0, [({"rg": 1.0, "gr": 1.0}, [0])]),  # X on qubit 0
    (1.0, [({"rg": 1.0, "gr": 1.0}, [1])]),  # X on qubit 1
]
X0_plus_X1 = MPO.from_operator_repr(
    eigenstates=("g", "r"), n_qudits=natoms, operations=operations_X0_plus_X1
)

print("as MPO format: \n",X0_plus_X1)

operations_X0_plus_X1 = [ (1.0, [({"rg": 1.0, "gr": 1.0}, [0])]), # X on qubit 0 (1.0, [({"rg": 1.0, "gr": 1.0}, [1])]), # X on qubit 1 ] X0_plus_X1 = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=operations_X0_plus_X1 ) print("as MPO format: \n",X0_plus_X1)

as MPO format: 
 [tensor([[[[0.+0.j, 1.+0.j],
          [1.+0.j, 0.+0.j]],

         [[1.+0.j, 0.+0.j],
          [0.+0.j, 1.+0.j]]]], dtype=torch.complex128), tensor([[[[1.+0.j],
          [0.+0.j]],

         [[0.+0.j],
          [1.+0.j]]],


        [[[0.+0.j],
          [1.+0.j]],

         [[1.+0.j],
          [0.+0.j]]]], dtype=torch.complex128)]

Example 1.5: Pauli Z on qubit 0¶

$Z_0 = |g \rangle \langle g| - |r \rangle \langle r|$

In [7]:

operations_Z0 = [(1.0, [({"gg": 1.0, "rr": -1.0}, [0])])]
Z0 = MPO.from_operator_repr(
    eigenstates=("g", "r"), n_qudits=natoms, operations=operations_Z0
)

operations_Z0 = [(1.0, [({"gg": 1.0, "rr": -1.0}, [0])])] Z0 = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=operations_Z0 )

Example 1.6: Complex operator - different operators on different qubits¶

$X_0 \otimes Z_1$ ($X$ on qubit 0, $Z$ on qubit 1)

In [8]:

operations_X0_Z1 = [
    (
        1.0,
        [
            ({"rg": 1.0, "gr": 1.0}, [0]),  # X on qubit 0
            ({"gg": 1.0, "rr": -1.0}, [1]),  # Z on qubit 1
        ],
    )
]
X0_Z1 = MPO.from_operator_repr(
    eigenstates=("g", "r"), n_qudits=natoms, operations=operations_X0_Z1
)

operations_X0_Z1 = [ ( 1.0, [ ({"rg": 1.0, "gr": 1.0}, [0]), # X on qubit 0 ({"gg": 1.0, "rr": -1.0}, [1]), # Z on qubit 1 ], ) ] X0_Z1 = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=operations_X0_Z1 )

Approach 2: Direct MPO Construction (Manual Tensor Definition)¶

For more control, you can build the MPO tensors directly.

Each MPO factor is a 4D tensor with shape:

(left_bond, output, input, right_bond)

For single-site operators without entanglement: shape (1, d, d, 1), d=2 (in case of leakage d=3 and in case of XY d = 2).

The matrix is embedded as:

tensor[0, :, :, 0] = matrix

In the Rydberg basis, the matrix indices are:

• index 0 corresponds to |g>
• index 1 corresponds to |r>

So a matrix [[a, b], [c, d]] represents:

a*|g><g| + b*|g><r| + c*|r><g| + d*|r><r|

Example 2.1: Number operator $|r \rangle \langle r|$ as a tensor¶

Matrix: [[0, 0], [0, 1]] because $|r \rangle \langle r|$ has 1 at position (1,1)

In [9]:

n_matrix = torch.tensor([[0.0, 0.0], [0.0, 1.0]], dtype=dtype)
n_tensor = n_matrix.reshape(1, 2, 2, 1)

# Identity for qubits where we don't want to measure
identity_matrix = torch.eye(2, dtype=dtype)
identity_tensor = identity_matrix.reshape(1, 2, 2, 1)

# n0 ⊗ Id (number operator on qubit 0, identity on qubit 1)
n0_mpo = MPO([n_tensor, identity_tensor])
print("\n2.1 Number operator on qubit 0 (direct construction):")
print(f"    n_matrix =\n{n_matrix}\n")
print("     MPO factors",n0_mpo.factors)

n_matrix = torch.tensor([[0.0, 0.0], [0.0, 1.0]], dtype=dtype) n_tensor = n_matrix.reshape(1, 2, 2, 1) # Identity for qubits where we don't want to measure identity_matrix = torch.eye(2, dtype=dtype) identity_tensor = identity_matrix.reshape(1, 2, 2, 1) # n0 ⊗ Id (number operator on qubit 0, identity on qubit 1) n0_mpo = MPO([n_tensor, identity_tensor]) print("\n2.1 Number operator on qubit 0 (direct construction):") print(f" n_matrix =\n{n_matrix}\n") print(" MPO factors",n0_mpo.factors)

2.1 Number operator on qubit 0 (direct construction):
    n_matrix =
tensor([[0.+0.j, 0.+0.j],
        [0.+0.j, 1.+0.j]], dtype=torch.complex128)

     MPO factors [tensor([[[[0.+0.j],
          [0.+0.j]],

         [[0.+0.j],
          [1.+0.j]]]], dtype=torch.complex128), tensor([[[[1.+0.j],
          [0.+0.j]],

         [[0.+0.j],
          [1.+0.j]]]], dtype=torch.complex128)]

Example 2.2: Pauli X operator¶

$X = |r \rangle \langle g| + |g \rangle \langle r|$ = [[0, 1], [1, 0]]

In [10]:

X_matrix = torch.tensor([[0.0, 1.0], [1.0, 0.0]], dtype=dtype)
X_tensor = X_matrix.reshape(1, 2, 2, 1)

# X0 ⊗ X1
XX_direct = MPO([X_tensor, X_tensor.clone()])
print("XX_direct MPO factors:\n", XX_direct.factors)

X_matrix = torch.tensor([[0.0, 1.0], [1.0, 0.0]], dtype=dtype) X_tensor = X_matrix.reshape(1, 2, 2, 1) # X0 ⊗ X1 XX_direct = MPO([X_tensor, X_tensor.clone()]) print("XX_direct MPO factors:\n", XX_direct.factors)

XX_direct MPO factors:
 [tensor([[[[0.+0.j],
          [1.+0.j]],

         [[1.+0.j],
          [0.+0.j]]]], dtype=torch.complex128), tensor([[[[0.+0.j],
          [1.+0.j]],

         [[1.+0.j],
          [0.+0.j]]]], dtype=torch.complex128)]

Example 2.3: Projection operators¶

In [11]:

# |g><g| projects onto ground state
proj_g_matrix = torch.tensor([[1.0, 0.0], [0.0, 0.0]], dtype=dtype)
proj_g_tensor = proj_g_matrix.reshape(1, 2, 2, 1)

# |r><r| projects onto Rydberg state
proj_r_matrix = torch.tensor([[0.0, 0.0], [0.0, 1.0]], dtype=dtype)
proj_r_tensor = proj_r_matrix.reshape(1, 2, 2, 1)

# |gg><gg| = |g><g|_0 ⊗ |g><g|_1
proj_gg = MPO([proj_g_tensor, proj_g_tensor.clone()])

# |rr><rr| = |r><r|_0 ⊗ |r><r|_1
proj_rr = MPO([proj_r_tensor, proj_r_tensor.clone()])

# |g>

Example 2.4: Raising and lowering operators¶

In [12]:

# sigma+ = |r><g| (raises from g to r)
sigma_plus_matrix = torch.tensor([[0.0, 0.0], [1.0, 0.0]], dtype=dtype)
sigma_plus_matrix = sigma_plus_matrix.reshape(1, 2, 2, 1)
# sigma- = |g><r| (lowers from r to g)
sigma_minus_matrix = torch.tensor([[0.0, 1.0], [0.0, 0.0]], dtype=dtype)
sigma_minus_matrix = sigma_minus_matrix.reshape(1, 2, 2, 1)

raising_x_raising_MPO = MPO([sigma_plus_matrix, sigma_plus_matrix.clone()])

# sigma+ = |r>

Example 2.5: Combining MPOs with addition and scalar multiplication¶

In [13]:

# You can add MPOs and multiply by scalars
# Example: (|gg><gg| + |rr><rr|) measures population in |gg> and |rr>
combined_proj = proj_gg + proj_rr

# Scalar multiplication: 2 * |rr><rr|
scaled_proj = 2.0 * proj_rr

# You can add MPOs and multiply by scalars # Example: (|gg> and |rr> combined_proj = proj_gg + proj_rr # Scalar multiplication: 2 * |rr>

Verification: Comparing Both Approaches¶

Let's verify that both approaches produce the same MPO tensors. This is important to ensure consistency between symbolic and direct methods.

In [14]:

def assert_mpo_equal(mpo1: MPO, mpo2: MPO, atol: float = 1e-10) -> None:
    """Assert that two MPOs have identical factors (up to numerical tolerance)."""
    assert len(mpo1.factors) == len(mpo2.factors)
    for _, (f1, f2) in enumerate(zip(mpo1.factors, mpo2.factors)):
        assert (
            f1.shape == f2.shape
        )
        assert torch.allclose(
            f1, f2, atol=atol
        )

def assert_mpo_equal(mpo1: MPO, mpo2: MPO, atol: float = 1e-10) -> None: """Assert that two MPOs have identical factors (up to numerical tolerance).""" assert len(mpo1.factors) == len(mpo2.factors) for _, (f1, f2) in enumerate(zip(mpo1.factors, mpo2.factors)): assert ( f1.shape == f2.shape ) assert torch.allclose( f1, f2, atol=atol )

In [15]:

# Build direct MPOs to compare with symbolic ones

# 1. Compare X0 (Pauli X on qubit 0)
X0_direct = MPO([X_tensor.clone(), identity_tensor.clone()])
assert_mpo_equal(X0, X0_direct)

# 2. Compare XX (X0 ⊗ X1 with coefficient 2)
XX_direct_scaled = 2.0 * MPO([X_tensor.clone(), X_tensor.clone()])
assert_mpo_equal(XX, XX_direct_scaled)

# 3. Compare n1 (number operator on qubit 1)
n1_direct = MPO([identity_tensor.clone(), n_tensor.clone()])
assert_mpo_equal(n1, n1_direct)

# 4. Compare Z0 (Pauli Z on qubit 0)
Z_matrix = torch.tensor([[1.0, 0.0], [0.0, -1.0]], dtype=dtype)
Z_tensor = Z_matrix.reshape(1, 2, 2, 1)
Z0_direct = MPO([Z_tensor.clone(), identity_tensor.clone()])
assert_mpo_equal(Z0, Z0_direct)

# 5. Compare X0_Z1 (X on qubit 0, Z on qubit 1)
X0_Z1_direct = MPO([X_tensor.clone(), Z_tensor.clone()])
assert_mpo_equal(X0_Z1, X0_Z1_direct)

# 6. Compare |rr><rr| projection
proj_rr_symbolic = MPO.from_operator_repr(
    eigenstates=("g", "r"),
    n_qudits=natoms,
    operations=[(1.0, [({"rr": 1.0}, [0, 1])])],
)
assert_mpo_equal(proj_rr_symbolic, proj_rr)

# 7. Compare |gg><gg| projection
proj_gg_symbolic = MPO.from_operator_repr(
    eigenstates=("g", "r"),
    n_qudits=natoms,
    operations=[(1.0, [({"gg": 1.0}, [0, 1])])],
)
assert_mpo_equal(proj_gg_symbolic, proj_gg)

# Build direct MPOs to compare with symbolic ones # 1. Compare X0 (Pauli X on qubit 0) X0_direct = MPO([X_tensor.clone(), identity_tensor.clone()]) assert_mpo_equal(X0, X0_direct) # 2. Compare XX (X0 ⊗ X1 with coefficient 2) XX_direct_scaled = 2.0 * MPO([X_tensor.clone(), X_tensor.clone()]) assert_mpo_equal(XX, XX_direct_scaled) # 3. Compare n1 (number operator on qubit 1) n1_direct = MPO([identity_tensor.clone(), n_tensor.clone()]) assert_mpo_equal(n1, n1_direct) # 4. Compare Z0 (Pauli Z on qubit 0) Z_matrix = torch.tensor([[1.0, 0.0], [0.0, -1.0]], dtype=dtype) Z_tensor = Z_matrix.reshape(1, 2, 2, 1) Z0_direct = MPO([Z_tensor.clone(), identity_tensor.clone()]) assert_mpo_equal(Z0, Z0_direct) # 5. Compare X0_Z1 (X on qubit 0, Z on qubit 1) X0_Z1_direct = MPO([X_tensor.clone(), Z_tensor.clone()]) assert_mpo_equal(X0_Z1, X0_Z1_direct) # 6. Compare |rr>

Running a simulation with custom observables¶

In [16]:

# Create several operators to measure
# 1. X0 ⊗ X1 using symbolic notation
XX_symbolic = MPO.from_operator_repr(
    eigenstates=("g", "r"),
    n_qudits=natoms,
    operations=[(1.0, [({"rg": 1.0, "gr": 1.0}, [0, 1])])],
)

# 2. Total number operator n0 + n1 (counts total Rydberg excitations)
total_n = MPO.from_operator_repr(
    eigenstates=("g", "r"),
    n_qudits=natoms,
    operations=[
        (1.0, [({"rr": 1.0}, [0])]),
        (1.0, [({"rr": 1.0}, [1])]),
    ],
)

# 3. |rr><rr| using direct construction (prob. of both atoms in Rydberg)
proj_rr_direct = MPO([proj_r_tensor, proj_r_tensor.clone()])

# Configure simulation with observables
dt = 10
eval_times = [0.5, 1.0]  # Evaluate at 50% and 100% of the simulation

config = MPSConfig(
    num_gpus_to_use=0,
    dt=dt,
    observables=[
        Expectation(operator=XX_symbolic, evaluation_times=eval_times, tag_suffix="XX"),
        Expectation(operator=total_n, evaluation_times=eval_times, tag_suffix="total_n"),
        Expectation(
            operator=proj_rr_direct,
            evaluation_times=eval_times,
            tag_suffix="rr_population",
        ),
    ],
    optimize_qubit_ordering=False,  # Required for custom Expectation observables
)

# Run simulation
simul = MPSBackend(seq, config=config)
results = simul.run()

# Print results
print("\nResults:")
print(f"  <X0 X1> at eval times: {results.expectation_XX}")
print(f"  <n0 + n1> (total excitations): {results.expectation_total_n}")
print(f"  |<rr|psi>|^2 (both Rydberg): {results.expectation_rr_population}")

# Create several operators to measure # 1. X0 ⊗ X1 using symbolic notation XX_symbolic = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=[(1.0, [({"rg": 1.0, "gr": 1.0}, [0, 1])])], ) # 2. Total number operator n0 + n1 (counts total Rydberg excitations) total_n = MPO.from_operator_repr( eigenstates=("g", "r"), n_qudits=natoms, operations=[ (1.0, [({"rr": 1.0}, [0])]), (1.0, [({"rr": 1.0}, [1])]), ], ) # 3. |rr> at eval times: {results.expectation_XX}") print(f" (total excitations): {results.expectation_total_n}") print(f" ||^2 (both Rydberg): {results.expectation_rr_population}")

emu-mps allows only {'energy', 'bitstrings', 'energy_second_moment', 'statistics', 'energy_variance', 'correlation_matrix', 'occupation'} observables with `optimize_qubit_ordering = True`. you provided unsupported {'expectation'} using `optimize_qubit_ordering = False` instead.
emu-mps allows only {'energy', 'bitstrings', 'energy_second_moment', 'statistics', 'energy_variance', 'correlation_matrix', 'occupation'} observables with `optimize_qubit_ordering = True`. you provided unsupported {'expectation'} using `optimize_qubit_ordering = False` instead.

Results:
  <X0 X1> at eval times: [tensor(0.+0.j, dtype=torch.complex128), tensor(1.4677e-39+0.j, dtype=torch.complex128)]
  <n0 + n1> (total excitations): [tensor(1.0000+5.7789e-36j, dtype=torch.complex128), tensor(2.+0.j, dtype=torch.complex128)]
  |<rr|psi>|^2 (both Rydberg): [tensor(0.2500+0.j, dtype=torch.complex128), tensor(1.+0.j, dtype=torch.complex128)]

Summary: Quick Reference¶

Operator Dictionary for Rydberg Basis¶

Basis: eigenstates=("g", "r")

• |g> = ground state = [1, 0]^T (index 0)
• |r> = Rydberg state = [0, 1]^T (index 1)

Symbolic notation (for MPO.from_operator_repr):

• "gg" = $|g \rangle \langle g|$ = [[1,0],[0,0]] (ground state projector)
• "rr" = $|r \rangle \langle r|$ = [[0,0],[0,1]] (Rydberg projector / number operator)
• "gr" = $|g \rangle \langle r|$ = [[0,1],[0,0]] (lowering operator sigma-)
• "rg" = $|r \rangle \langle g|$ = [[0,0],[1,0]] (raising operator sigma+)

Common operators:

• Identity ($Id$): {"gg": 1.0, "rr": 1.0}
• Pauli X: {"rg": 1.0, "gr": 1.0}
• Pauli Y: {"rg": -1j, "gr": 1j}
• Pauli Z: {"gg": 1.0, "rr": -1.0}
• Number (n): {"rr": 1.0}
• sigma+: {"rg": 1.0}
• sigma-: {"gr": 1.0}

Operations structure:

operations = [(coeff, [(op_dict, [qubit_indices]), ...]), ...]

Example - X on qubit 0, Z on qubit 1:

operations = [(1.0, [({"rg": 1.0, "gr": 1.0}, [0]),
                     ({"gg": 1.0, "rr": -1.0}, [1])])]

Example - Sum $X_0 + X_1$:

operations = [(1.0, [({"rg": 1.0, "gr": 1.0}, [0])]),
              (1.0, [({"rg": 1.0, "gr": 1.0}, [1])])]