Adding relaxation noise channel (#675)

* Adding relaxation noise channel

* Improving noise channel descriptions

* Switch relaxation and dephasing

* Rephrase depolarizing noise description

pulser-core/pulser/backend/noise_model.py

@@ -20,7 +20,13 @@
 import numpy as np
 
 NOISE_TYPES = Literal[
-    "doppler", "amplitude", "SPAM", "dephasing", "depolarizing", "eff_noise"
+    "doppler",
+    "amplitude",
+    "SPAM",
+    "dephasing",
+    "relaxation",
+    "depolarizing",
+    "eff_noise",
 ]
 
 
@@ -37,11 +43,17 @@ class NoiseModel:
         noise_types: Noise types to include in the emulation. Available
             options:
 
+            - "relaxation": Noise due to a decay from the Rydberg to
+              the ground state (parametrized by `relaxation_rate`), commonly
+              characterized experimentally by the T1 time.
             - "dephasing": Random phase (Z) flip (parametrized
-              by `dephasing_rate`).
+              by `dephasing_rate`), commonly characterized experimentally
+              by the T2 time.
             - "depolarizing": Quantum noise where the state is
-              turned into a mixed state I/2 with rate
-              `depolarizing_rate`.
+              turned into a mixed state I/2 with rate `depolarizing_rate`.
+              While it does not describe a physical phenomenon, it is a
+              commonly used tool to test the system under a uniform
+              combination of phase flip (Z) and bit flip (X) errors.
             - "eff_noise": General effective noise channel defined by
               the set of collapse operators `eff_noise_opers`
               and the corresponding rates distribution
@@ -67,11 +79,14 @@ class NoiseModel:
             pulses.
         amp_sigma: Dictates the fluctuations in amplitude as a standard
             deviation of a normal distribution centered in 1.
-        dephasing_rate: The rate of a dephasing error occuring (in rad/µs).
-        depolarizing_rate: The rate (in rad/µs) at which a depolarizing
+        relaxation_rate: The rate of relaxation from the Rydberg to the
+            ground state (in 1/µs). Corresponds to 1/T1.
+        dephasing_rate: The rate of a dephasing error occuring (in 1/µs).
+            Corresponds to 1/T2.
+        depolarizing_rate: The rate (in 1/µs) at which a depolarizing
             error occurs.
         eff_noise_rates: The rate associated to each effective noise operator
-            (in rad/µs).
+            (in 1/µs).
         eff_noise_opers: The operators for the effective noise model.
     """
 
@@ -84,6 +99,7 @@ class NoiseModel:
     temperature: float = 50.0
     laser_waist: float = 175.0
     amp_sigma: float = 5e-2
+    relaxation_rate: float = 0.01
     dephasing_rate: float = 0.05
     depolarizing_rate: float = 0.05
     eff_noise_rates: list[float] = field(default_factory=list)
@@ -92,6 +108,7 @@ class NoiseModel:
     def __post_init__(self) -> None:
         positive = {
             "dephasing_rate",
+            "relaxation_rate",
             "depolarizing_rate",
         }
         strict_positive = {

pulser-simulation/pulser_simulation/hamiltonian.py

@@ -119,6 +119,16 @@ def basis_check(noise_type: str) -> None:
             coeff = np.sqrt(config.dephasing_rate / 2)
             local_collapse_ops.append(coeff * qutip.sigmaz())
 
+        if "relaxation" in config.noise_types:
+            coeff = np.sqrt(config.relaxation_rate)
+            try:
+                local_collapse_ops.append(coeff * self.op_matrix["sigma_gr"])
+            except KeyError:
+                raise ValueError(
+                    "'relaxation' noise requires addressing of the"
+                    " 'ground-rydberg' basis."
+                )
+
         if "depolarizing" in config.noise_types:
             basis_check("depolarizing")
             coeff = np.sqrt(config.depolarizing_rate / 4)
@@ -175,7 +185,7 @@ def _extract_samples(self) -> None:
         """Populates samples dictionary with every pulse in the sequence."""
         local_noises = True
         if set(self.config.noise_types).issubset(
-            {"dephasing", "SPAM", "depolarizing", "eff_noise"}
+            {"dephasing", "relaxation", "SPAM", "depolarizing", "eff_noise"}
         ):
             local_noises = (
                 "SPAM" in self.config.noise_types

pulser-simulation/pulser_simulation/simconfig.py

@@ -17,15 +17,12 @@
 
 from dataclasses import dataclass, field
 from math import sqrt
-from typing import Any, Literal, Optional, Tuple, Type, TypeVar, Union, cast
+from typing import Any, Optional, Tuple, Type, TypeVar, Union, cast
 
 import qutip
 
-from pulser.backend.noise_model import NoiseModel
+from pulser.backend.noise_model import NOISE_TYPES, NoiseModel
 
-NOISE_TYPES = Literal[
-    "doppler", "amplitude", "SPAM", "dephasing", "depolarizing", "eff_noise"
-]
 MASS = 1.45e-25  # kg
 KB = 1.38e-23  # J/K
 KEFF = 8.7  # µm^-1
@@ -36,6 +33,7 @@
     "ising": {
         "amplitude",
         "dephasing",
+        "relaxation",
         "depolarizing",
         "doppler",
         "eff_noise",
@@ -72,6 +70,7 @@ class SimConfig:
             simulation. You may specify just one, or a tuple of the allowed
             noise types:
 
+            - "relaxation": Relaxation from the Rydberg to the ground state.
             - "dephasing": Random phase (Z) flip.
             - "depolarizing": Quantum noise where the state (rho) is
               turned into a mixed state I/2 at a rate gamma (in rad/µs).
@@ -109,6 +108,7 @@ class SimConfig:
     eta: float = 0.005
     epsilon: float = 0.01
     epsilon_prime: float = 0.05
+    relaxation_rate: float = 0.01
     dephasing_rate: float = 0.05
     depolarizing_rate: float = 0.05
     eff_noise_rates: list[float] = field(default_factory=list, repr=False)
@@ -129,6 +129,7 @@ def from_noise_model(cls: Type[T], noise_model: NoiseModel) -> T:
             epsilon=noise_model.p_false_pos,
             epsilon_prime=noise_model.p_false_neg,
             dephasing_rate=noise_model.dephasing_rate,
+            relaxation_rate=noise_model.relaxation_rate,
             depolarizing_rate=noise_model.depolarizing_rate,
             eff_noise_rates=noise_model.eff_noise_rates,
             eff_noise_opers=list(map(qutip.Qobj, noise_model.eff_noise_opers)),
@@ -147,6 +148,7 @@ def to_noise_model(self) -> NoiseModel:
             laser_waist=self.laser_waist,
             amp_sigma=self.amp_sigma,
             dephasing_rate=self.dephasing_rate,
+            relaxation_rate=self.relaxation_rate,
             depolarizing_rate=self.depolarizing_rate,
             eff_noise_rates=self.eff_noise_rates,
             eff_noise_opers=[op.full() for op in self.eff_noise_opers],
@@ -209,6 +211,8 @@ def __str__(self, solver_options: bool = False) -> str:
         if "amplitude" in self.noise:
             lines.append(f"Laser waist:           {self.laser_waist}μm")
             lines.append(f"Amplitude standard dev.:  {self.amp_sigma}")
+        if "relaxation" in self.noise:
+            lines.append(f"Relaxation rate: {self.relaxation_rate}")
         if "dephasing" in self.noise:
             lines.append(f"Dephasing rate: {self.dephasing_rate}")
         if "depolarizing" in self.noise:

pulser-simulation/pulser_simulation/simresults.py

@@ -52,7 +52,7 @@ def __init__(
             size: The number of atoms in the register.
             basis_name: The basis indicating the addressed atoms after
                 the pulse sequence ('ground-rydberg', 'digital' or 'all').
-            sim_times: Array of times (µs) when simulation results are
+            sim_times: Array of times (in µs) when simulation results are
                 returned.
         """
         self._dim = 3 if basis_name == "all" else 2
@@ -132,7 +132,7 @@ def sample_state(
         """Returns the result of multiple measurements at time t.
 
         Args:
-            t: Time at which the state is sampled.
+            t: Time at which the state is sampled (in µs).
             n_samples: Number of samples to return.
             t_tol: Tolerance for the difference between t and
                 closest time.
@@ -291,7 +291,7 @@ def get_state(self, t: float, t_tol: float = 1.0e-3) -> qutip.Qobj:
             way of computing expectation values of observables.
 
         Args:
-            t: Time (µs) at which to return the state.
+            t: Time (in µs) at which to return the state.
             t_tol: Tolerance for the difference between t and
                 closest time.
 
@@ -423,7 +423,7 @@ def get_state(
         """Get the state at time t of the simulation.
 
         Args:
-            t: Time (µs) at which to return the state.
+            t: Time (in µs) at which to return the state.
             reduce_to_basis: Reduces the full state vector
                 to the given basis ("ground-rydberg" or "digital"), if the
                 population of the states to be ignored is negligible. Doesn't
@@ -515,7 +515,7 @@ def sample_state(
         """Returns the result of multiple measurements at time t.
 
         Args:
-            t: Time at which the state is sampled.
+            t: Time (in µs) at which the state is sampled.
             n_samples: Number of samples to return.
             t_tol: Tolerance for the difference between t and
                 closest time.

pulser-simulation/pulser_simulation/simulation.py

@@ -263,6 +263,8 @@ def add_config(self, config: SimConfig) -> None:
             param_dict["amp_sigma"] = noise_model.amp_sigma
         if "dephasing" in diff_noise_set:
             param_dict["dephasing_rate"] = noise_model.dephasing_rate
+        if "relaxation" in diff_noise_set:
+            param_dict["relaxation_rate"] = noise_model.relaxation_rate
         if "depolarizing" in diff_noise_set:
             param_dict["depolarizing_rate"] = noise_model.depolarizing_rate
         if "eff_noise" in diff_noise_set:
@@ -542,6 +544,7 @@ def _run_solver() -> CoherentResults:
 
             if (
                 "dephasing" in self.config.noise
+                or "relaxation" in self.config.noise
                 or "depolarizing" in self.config.noise
                 or "eff_noise" in self.config.noise
             ):
@@ -581,7 +584,7 @@ def _run_solver() -> CoherentResults:
 
         # Check if noises ask for averaging over multiple runs:
         if set(self.config.noise).issubset(
-            {"dephasing", "SPAM", "depolarizing", "eff_noise"}
+            {"dephasing", "relaxation", "SPAM", "depolarizing", "eff_noise"}
         ):
             # If there is "SPAM", the preparation errors must be zero
             if "SPAM" not in self.config.noise or self.config.eta == 0:

tests/test_backend.py

@@ -110,19 +110,15 @@ def test_init_strict_pos(self, param):
         "param",
         [
             "dephasing_rate",
+            "relaxation_rate",
             "depolarizing_rate",
         ],
     )
     def test_init_rate_like(self, param, value):
         if value < 0:
-            param_mess = (
-                "depolarizing_rate"
-                if "depolarizing" in param
-                else "dephasing_rate"
-            )
             with pytest.raises(
                 ValueError,
-                match=f"'{param_mess}' must be None or greater "
+                match=f"'{param}' must be None or greater "
                 f"than or equal to zero, not {value}.",
             ):
                 NoiseModel(**{param: value})
@@ -132,6 +128,8 @@ def test_init_rate_like(self, param, value):
                 assert noise_model.depolarizing_rate == value
             elif "dephasing" in param:
                 assert noise_model.dephasing_rate == value
+            elif "relaxation" in param:
+                assert noise_model.relaxation_rate == value
 
     @pytest.mark.parametrize("value", [-1e-9, 1.0001])
     @pytest.mark.parametrize(

tests/test_simconfig.py

@@ -49,11 +49,13 @@ def test_init():
         and "100" in str_config
         and "Solver Options" in str_config
     )
-    config = SimConfig(noise="depolarizing")
+    config = SimConfig(noise=("depolarizing", "relaxation"))
     assert config.temperature == 5e-5
     assert config.to_noise_model().temperature == 50
     str_config = config.__str__(True)
-    assert "depolarizing" in str_config
+    assert "depolarizing" in str_config and "relaxation" in str_config
+    assert f"Depolarizing rate: {config.depolarizing_rate}" in str_config
+    assert f"Relaxation rate: {config.relaxation_rate}" in str_config
     config = SimConfig(
         noise="eff_noise",
         eff_noise_opers=[qeye(2), sigmax()],

tests/test_simulation.py

@@ -746,9 +746,10 @@ def test_noise_with_zero_epsilons(seq, matrices):
     "noise, result, n_collapse_ops",
     [
         ("dephasing", {"0": 595, "1": 405}, 1),
+        ("relaxation", {"0": 595, "1": 405}, 1),
         ("eff_noise", {"0": 595, "1": 405}, 1),
         ("depolarizing", {"0": 587, "1": 413}, 3),
-        (("dephasing", "depolarizing"), {"0": 587, "1": 413}, 4),
+        (("dephasing", "depolarizing", "relaxation"), {"0": 587, "1": 413}, 5),
         (("eff_noise", "dephasing"), {"0": 595, "1": 405}, 2),
     ],
 )
@@ -778,6 +779,25 @@ def test_noises_rydberg(matrices, noise, result, n_collapse_ops):
     assert np.trace(trace_2) < 1 and not np.isclose(np.trace(trace_2), 1)
 
 
+def test_relaxation_noise():
+    seq = Sequence(Register({"q0": (0, 0)}), MockDevice)
+    seq.declare_channel("ryd", "rydberg_global")
+    seq.add(Pulse.ConstantDetuning(BlackmanWaveform(1000, np.pi), 0, 0), "ryd")
+    seq.delay(10000, "ryd")
+
+    sim = QutipEmulator.from_sequence(seq)
+    sim.add_config(SimConfig(noise="relaxation", relaxation_rate=0.1))
+    res = sim.run()
+    start_samples = res.sample_state(1)
+    ryd_pop = start_samples["1"]
+    assert ryd_pop > start_samples.get("0", 0)
+    # The Rydberg state population gradually decays
+    for t_ in range(2, 10):
+        new_ryd_pop = res.sample_state(t_)["1"]
+        assert new_ryd_pop < ryd_pop
+        ryd_pop = new_ryd_pop
+
+
 depo_res = {
     "111": 821,
     "110": 61,
@@ -822,6 +842,13 @@ def test_noises_digital(matrices, noise, result, n_collapse_ops, seq_digital):
             eff_noise_rates=[0.025],
         ),
     )
+
+    with pytest.raises(
+        ValueError,
+        match="'relaxation' noise requires addressing of the 'ground-rydberg'",
+    ):
+        sim.set_config(SimConfig(noise="relaxation"))
+
     res = sim.run()
     res_samples = res.sample_final_state()
     assert res_samples == Counter(result)