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:

MPO.from_operator_repr(): Use symbolic notation with basis state labels
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

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")

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}")

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}")

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
)

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)

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
)

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
)

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)

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)

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>
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>
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><g|_1
proj_gg = MPO([proj_g_tensor, proj_g_tensor.clone()])

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

Example 2.4: Raising and lowering operators

In [12]:

# sigma+ = |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>
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()])

Example 2.5: Combining MPOs with addition and scalar multiplication

In [13]:

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

# Scalar multiplication: 2 * |rr><rr|
scaled_proj = 2.0 * proj_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
        )

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>
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>
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)

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>
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"   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])])]

Previous
Running noisy simulations Next
Overview