diff --git a/qiskit_aer/noise/__init__.py b/qiskit_aer/noise/__init__.py
index e26aa70626..ccc3cf71ae 100644
--- a/qiskit_aer/noise/__init__.py
+++ b/qiskit_aer/noise/__init__.py
@@ -181,6 +181,7 @@
 
     NoiseModel
     QuantumError
+    PauliLindbladError
     ReadoutError
 
 
@@ -240,6 +241,7 @@
 from .noise_model import NoiseModel
 from .errors import QuantumError
 from .errors import PauliError
+from .errors import PauliLindbladError
 from .errors import ReadoutError
 
 # Error generating functions
diff --git a/qiskit_aer/noise/errors/__init__.py b/qiskit_aer/noise/errors/__init__.py
index 3472ac7c2f..b4ff38136a 100644
--- a/qiskit_aer/noise/errors/__init__.py
+++ b/qiskit_aer/noise/errors/__init__.py
@@ -17,6 +17,7 @@
 from .readout_error import ReadoutError
 from .quantum_error import QuantumError
 from .pauli_error import PauliError
+from .pauli_lindblad_error import PauliLindbladError
 from .standard_errors import kraus_error
 from .standard_errors import mixed_unitary_error
 from .standard_errors import coherent_unitary_error
diff --git a/qiskit_aer/noise/errors/pauli_lindblad_error.py b/qiskit_aer/noise/errors/pauli_lindblad_error.py
new file mode 100644
index 0000000000..136a1d53ce
--- /dev/null
+++ b/qiskit_aer/noise/errors/pauli_lindblad_error.py
@@ -0,0 +1,363 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018-2024.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Class for representing a Pauli noise channel generated by a Pauli Lindblad dissipator.
+"""
+
+from __future__ import annotations
+from collections.abc import Sequence
+import numpy as np
+
+from qiskit.quantum_info import Pauli, PauliList, SparsePauliOp, ScalarOp, SuperOp
+from qiskit.quantum_info.operators.mixins import TolerancesMixin
+from .base_quantum_error import BaseQuantumError
+from .quantum_error import QuantumError
+from .pauli_error import PauliError, sort_paulis
+from ..noiseerror import NoiseError
+
+
+class PauliLindbladError(BaseQuantumError, TolerancesMixin):
+    r"""A Pauli channel generated by a Pauli Lindblad dissipator.
+
+    This operator represents an N-qubit quantum error channel
+    :math:`E(ρ) = e^{\sum_j r_j D_{P_j}}(ρ)` generated by Pauli Lindblad dissipators
+    :math:`D_P(ρ) = P ρ P - ρ`, where :math:`P_j` are N-qubit :class:`~.Pauli` operators.
+
+    The list of Pauli generator terms are stored as a :class:`~.PauliList` and can be
+    accessed via the :attr:`generators` attribute. The array of dissipator rates
+    :math:`r_j` can be accessed via the :attr:`rates` attribute.
+
+    A Pauli lindblad error is equivalent to a :class:`.PauliError` and can be converted
+    using the :meth:`to_pauli_error` method. Though note that for a sparse generated
+    ``PauliLindbladError`` there may in general be exponentially many terms in the
+    converted :class:`.PauliError` operator. Because of this, this operator can be
+    used to more efficiently represent N-qubit Pauli channels for simulation if they
+    have few generator terms.
+
+    The equivalent Pauli error is constructed as a composition of single-Pauli channel terms
+
+    .. math::
+
+        e^{\sum_j r_j D_{P_j}} = \prod_j e^{r_j D_{P_j}}
+        = prod_j \left( (1 - p_j) S_I + p_j S_{P_j} \right)
+
+    where :math:`p_j = \frac12 - \frac12 e^{-2 r_j}`.
+
+    .. note::
+
+        This operator can also represent a non-physical (non-CPTP) channel if any of the
+        rates are negative. Non-physical operators cannot be converted to a
+        :class:`~.QuantumError` or used in an :class:`~.AerSimulator` simulation. You can
+        check if an operator is physical using the :meth:`is_cptp` method.
+    """
+
+    def __init__(
+        self,
+        generators: Sequence[Pauli],
+        rates: Sequence[float],
+    ):
+        """Initialize a Pauli-Lindblad dissipator model.
+
+        Args:
+            generators: A list of Pauli's corresponding to the Lindblad dissipator generator terms.
+            rates: The Pauli Lindblad dissipator generator rates.
+        """
+        self._generators = PauliList(generators)
+        self._rates = np.asarray(rates, dtype=float)
+        if self._rates.shape != (len(self._generators),):
+            raise NoiseError("Input generators and rates are different lengths.")
+        super().__init__(num_qubits=self._generators.num_qubits)
+
+    def __repr__(self):
+        return f"{type(self).__name__}({self.generators.to_labels()}, {self.rates.tolist()})"
+
+    def __eq__(self, other):
+        # Use BaseOperator eq to check type and shape
+        if not super().__eq__(other):
+            return False
+        lhs = self.simplify()
+        rhs = other.simplify()
+        if lhs.size != rhs.size:
+            return False
+        lpaulis, lrates = sort_paulis(lhs.generators, lhs.rates)
+        rpaulis, rrates = sort_paulis(rhs.generators, rhs.rates)
+        return np.allclose(lrates, rrates) and lpaulis == rpaulis
+
+    @property
+    def size(self):
+        """Return the number of error circuit."""
+        return len(self.generators)
+
+    @property
+    def generators(self) -> PauliList:
+        """Return the Pauli Lindblad dissipator generator terms"""
+        return self._generators
+
+    @property
+    def rates(self) -> np.ndarray:
+        """Return the Lindblad dissipator generator rates"""
+        return self._rates
+
+    @property
+    def settings(self):
+        """Settings for IBM RuntimeEncoder JSON encoding"""
+        return {
+            "generators": self.generators,
+            "rates": self.rates,
+        }
+
+    def ideal(self) -> bool:
+        """Return True if this error object is composed only of identity operations.
+        Note that the identity check is best effort and up to global phase."""
+        return np.allclose(self.rates, 0) or not (
+            np.any(self.generators.z) or np.any(self.generators.x)
+        )
+
+    def is_cptp(self, atol: float | None = None, rtol: float | None = None) -> bool:
+        """Return True if completely-positive trace-preserving (CPTP)."""
+        return self.is_cp(atol=atol, rtol=rtol)
+
+    def is_tp(self, atol: float | None = None, rtol: float | None = None) -> bool:
+        """Test if a channel is trace-preserving (TP)"""
+        # pylint: disable = unused-argument
+        # This error is TP by construction regardless of rates.
+        return True
+
+    def is_cp(self, atol: float | None = None, rtol: float | None = None) -> bool:
+        """Test if Choi-matrix is completely-positive (CP)"""
+        if atol is None:
+            atol = self.atol
+        if rtol is None:
+            rtol = self.rtol
+        neg_rates = self.rates[self.rates < 0]
+        return np.allclose(neg_rates, 0, atol=atol, rtol=rtol)
+
+    def tensor(self, other: PauliLindbladError) -> PauliLindbladError:
+        if not isinstance(other, PauliLindbladError):
+            raise NoiseError("other must be a PauliLindbladError")
+        zeros_l = np.zeros(self.num_qubits, dtype=bool)
+        zeros_r = np.zeros(other.num_qubits, dtype=bool)
+        gens_left = self.generators.tensor(Pauli((zeros_r, zeros_r)))
+        gens_right = other.generators.expand(Pauli((zeros_l, zeros_l)))
+        generators = gens_left + gens_right
+        rates = np.concatenate([self.rates, other.rates])
+        return PauliLindbladError(generators, rates)
+
+    def expand(self, other) -> PauliLindbladError:
+        if not isinstance(other, PauliLindbladError):
+            raise NoiseError("other must be a PauliLindbladError")
+        return other.tensor(self)
+
+    def compose(self, other, qargs=None, front=False) -> PauliLindbladError:
+        if qargs is None:
+            qargs = getattr(other, "qargs", None)
+        if not isinstance(other, PauliLindbladError):
+            raise NoiseError("other must be a PauliLindbladError")
+        # Validate composition dimensions and qargs match
+        self._op_shape.compose(other._op_shape, qargs, front)
+        # pylint: disable = unused-argument
+        padded = ScalarOp(self.num_qubits * (2,)).compose(
+            SparsePauliOp(other.generators, copy=False, ignore_pauli_phase=True), qargs
+        )
+        generators = self.generators + padded.paulis
+        rates = np.concatenate([self.rates, other.rates])
+        return PauliLindbladError(generators, rates)
+
+    def power(self, n: float) -> PauliLindbladError:
+        return PauliLindbladError(self.generators, n * self.rates)
+
+    def inverse(self) -> PauliLindbladError:
+        """Return the inverse (non-CPTP) channel"""
+        return PauliLindbladError(self.generators, -self.rates)
+
+    def simplify(
+        self, atol: float | None = None, rtol: float | None = None
+    ) -> "PauliLindbladError":
+        """Simplify PauliList by combining duplicates and removing zeros.
+
+        Args:
+            atol (float): Optional. Absolute tolerance for checking if
+                          coefficients are zero (Default: 1e-8).
+            rtol (float): Optional. relative tolerance for checking if
+                          coefficients are zero (Default: 1e-5).
+
+        Returns:
+            SparsePauliOp: the simplified SparsePauliOp operator.
+        """
+        if atol is None:
+            atol = self.atol
+        if rtol is None:
+            rtol = self.rtol
+        # Remove identity terms
+        non_iden = np.any(np.logical_or(self.generators.z, self.generators.x), axis=1)
+        simplified = SparsePauliOp(
+            self.generators[non_iden],
+            self.rates[non_iden],
+            copy=False,
+            ignore_pauli_phase=True,
+        ).simplify(atol=atol, rtol=rtol)
+        return PauliLindbladError(simplified.paulis, simplified.coeffs.real)
+
+    def subsystem_errors(self) -> list[tuple[PauliLindbladError, tuple[int, ...]]]:
+        """Return a list errors for the subsystem decomposed error.
+
+        .. note::
+
+            This uses a greedy algorithm to find the largest non-identity subsystem
+            Pauli, checks if its non identity terms are covered by any previously
+            selected Pauli's, and if not adds it to list of covering subsystems.
+            This is repeated until no generators remain.
+
+            In terms of the number of Pauli terms this has runtime
+            `O(num_terms * num_coverings)`,
+            which in the worst case is `O(num_terms ** 2)`.
+
+        Returns:
+            A list of pairs of component PauliLindbladErrors and subsystem indices for the
+            that decompose the current errors.
+        """
+        # Find non-identity paulis we wish to cover
+        paulis = self.generators
+        non_iden = np.logical_or(paulis.z, paulis.x)
+
+        # Mask to exclude generator terms that are trivial (identities)
+        non_trivial = np.any(non_iden, axis=1)
+
+        # Paulis that aren't covered yet
+        uncovered = np.arange(non_iden.shape[0])[non_trivial]
+
+        # Indices that cover each Pauli
+        # Note that trivial terms will be left at -1
+        covered = -1 * np.ones(non_iden.shape[0], dtype=int)
+        coverings = []
+
+        # In general finding optimal coverings is NP-hard (I think?)
+        # so we do a heuristic of just greedily finding the largest
+        # pauli that isn't covered, and checking if it is covered
+        # by any previous coverings, if not adding it to coverings
+        # and repeat until we run out of paulis
+        while uncovered.size:
+            imax = np.argmax(np.sum(non_iden[uncovered], axis=1))
+            rmax = uncovered[imax]
+            add_covering = True
+            for rcov in coverings:
+                if np.all(non_iden[rcov][non_iden[rmax]]):
+                    add_covering = False
+                    covered[rmax] = rcov
+                    break
+            if add_covering:
+                covered[rmax] = rmax
+                coverings.append(rmax)
+            uncovered = uncovered[uncovered != rmax]
+
+        # Extract subsystem errors and qinds of non-identity terms
+        sub_errors = []
+        for cover in coverings:
+            pinds = covered == cover
+            qinds = np.any(non_iden[pinds], axis=0)
+            sub_z = paulis.z[pinds][:, qinds]
+            sub_x = paulis.x[pinds][:, qinds]
+            sub_gens = PauliList.from_symplectic(sub_z, sub_x)
+            sub_err = PauliLindbladError(sub_gens, self.rates[pinds])
+            sub_qubits = tuple(np.where(qinds)[0])
+            sub_errors.append((sub_err, sub_qubits))
+
+        return sub_errors
+
+    def to_pauli_error(self, simplify: bool = True) -> PauliError:
+        """Convert to a PauliError operator.
+
+        .. note::
+
+            If this objects represents an non-CPTP inverse channel with negative
+            rates the returned "probabilities" will be a quasi-probability
+            distribution containing negative values.
+
+        Args:
+            simplify: If True call :meth:`~.PauliError.simplify` each single
+                Pauli channel composition to reduce the number of duplicate terms.
+
+        Returns:
+            The :class:`~.PauliError` of the current Pauli channel.
+        """
+        chan_z = np.zeros((1, self.num_qubits), dtype=bool)
+        chan_x = np.zeros_like(chan_z)
+        chan_p = np.ones(1, dtype=float)
+        for term_z, term_x, term_r in zip(self.generators.z, self.generators.x, self.rates):
+            term_p = 0.5 - 0.5 * np.exp(-2 * term_r)
+            chan_z = np.concatenate([chan_z, np.logical_xor(chan_z, term_z)], axis=0)
+            chan_x = np.concatenate([chan_x, chan_x ^ term_x])
+            chan_p = np.concatenate([(1 - term_p) * chan_p, term_p * chan_p])
+            if simplify:
+                error_op = PauliError(PauliList.from_symplectic(chan_z, chan_x), chan_p).simplify()
+                chan_z, chan_x, chan_p = (
+                    error_op.paulis.z,
+                    error_op.paulis.x,
+                    error_op.probabilities,
+                )
+        return PauliError(PauliList.from_symplectic(chan_z, chan_x), chan_p)
+
+    def to_quantum_error(self) -> QuantumError:
+        """Convert to a general QuantumError object."""
+        if not self.is_cptp():
+            raise NoiseError("Cannot convert non-CPTP PauliLindbladError to a QuantumError")
+        return self.to_pauli_error().to_quantum_error()
+
+    def to_quantumchannel(self) -> SuperOp:
+        """Convert to a dense N-qubit QuantumChannel"""
+        return self.to_pauli_error().to_quantumchannel()
+
+    def to_dict(self):
+        """Return the current error as a dictionary."""
+        # Assemble noise circuits for Aer simulator
+        qubits = list(range(self.num_qubits))
+        instructions = [
+            [{"name": "pauli", "params": [pauli.to_label()], "qubits": qubits}]
+            for pauli in self.generators
+        ]
+        # Construct error dict
+        error = {
+            "type": "plerror",
+            "id": self.id,
+            "operations": [],
+            "instructions": instructions,
+            "rates": self.rates.tolist(),
+        }
+        return error
+
+    @staticmethod
+    def from_dict(error: dict) -> PauliLindbladError:
+        """Implement current error from a dictionary."""
+        # check if dictionary
+        if not isinstance(error, dict):
+            raise NoiseError("error is not a dictionary")
+        # check expected keys "type, id, operations, instructions, probabilities"
+        if (
+            ("type" not in error)
+            or ("id" not in error)
+            or ("operations" not in error)
+            or ("instructions" not in error)
+            or ("rates" not in error)
+        ):
+            raise NoiseError("error dictionary not containing expected keys")
+        instructions = error["instructions"]
+        rates = error["rates"]
+        if len(instructions) != len(rates):
+            raise NoiseError("rates not matching with instructions")
+        # parse instructions and turn to noise_ops
+        generators = []
+        for inst in instructions:
+            if len(inst) != 1 or inst[0]["name"] != "pauli":
+                raise NoiseError("Invalid PauliLindbladError dict")
+            generators.append(inst[0]["params"][0])
+        return PauliLindbladError(generators, rates)
diff --git a/releasenotes/notes/pauli-lindblad-error-2fccae840ceec2aa.yaml b/releasenotes/notes/pauli-lindblad-error-2fccae840ceec2aa.yaml
new file mode 100644
index 0000000000..8271eaafc6
--- /dev/null
+++ b/releasenotes/notes/pauli-lindblad-error-2fccae840ceec2aa.yaml
@@ -0,0 +1,8 @@
+---
+features:
+  - |
+    Adds a new :class:`.PauliLindbladError` quantum error subclass.
+    This class represents a Pauli error quantum channel generated by 
+    Pauli-Lindblad dissipator terms. This allows for more efficient
+    representation and error sampling for large qubit Pauli channels
+    generated by sparse Pauli-lindblad terms.
diff --git a/src/noise/noise_model.hpp b/src/noise/noise_model.hpp
index 4063aaa58a..e377a5a1fd 100644
--- a/src/noise/noise_model.hpp
+++ b/src/noise/noise_model.hpp
@@ -1047,9 +1047,11 @@ void NoiseModel::load_from_json(const json_t &js) {
       throw std::invalid_argument(
           "Invalid noise_params JSON: \"error\" field is not a list");
     }
+    std::unordered_set<std::string> qerror_types = {"qerror", "plerror"};
     for (const auto &gate_js : JSON::get_value("errors", js)) {
       std::string type;
       JSON::get_value(type, "type", gate_js);
+      bool is_qerror = (qerror_types.find(type) != qerror_types.end());
       stringset_t
           ops; // want set so ops are unique, and we can pull out measure
       JSON::get_value(ops, "operations", gate_js);
@@ -1062,7 +1064,7 @@ void NoiseModel::load_from_json(const json_t &js) {
       // before the measure operation, rather than after like the other gates
       if (ops.find("measure") != ops.end() && type != "roerror") {
         ops.erase("measure"); // remove measure from set of ops
-        if (type != "qerror")
+        if (!is_qerror)
           throw std::invalid_argument("NoiseModel: Invalid noise type (" +
                                       type + ")");
         QuantumError error;
@@ -1071,7 +1073,7 @@ void NoiseModel::load_from_json(const json_t &js) {
         add_quantum_error(error, {"measure"}, gate_qubits, noise_qubits);
       }
       // Load the remaining ops as errors that come after op
-      if (type == "qerror") {
+      if (is_qerror) {
         QuantumError error;
         error.load_from_json(gate_js);
         add_quantum_error(error, ops, gate_qubits, noise_qubits);
diff --git a/src/noise/quantum_error.hpp b/src/noise/quantum_error.hpp
index b791f5bfe5..4dfa9324c7 100644
--- a/src/noise/quantum_error.hpp
+++ b/src/noise/quantum_error.hpp
@@ -17,6 +17,7 @@
 
 #include "framework/noise_utils.hpp"
 #include "framework/opset.hpp"
+#include "simulators/stabilizer/pauli.hpp"
 #include "simulators/superoperator/superoperator_state.hpp"
 
 namespace AER {
@@ -70,6 +71,11 @@ class QuantumError {
   void set_circuits(const std::vector<NoiseOps> &circuits,
                     const rvector_t &probs);
 
+  // Sets the generator Pauli sub-circuits and rates to be sampled from.
+  // The length of the circuits vector and probability vector must be equal.
+  void set_generators(const std::vector<NoiseOps> &circuits,
+                      const rvector_t &rates);
+
   // Construct a quantum error from a set of Kraus matrices
   // This will factor out any identity or unitary Kraus operators into
   // non-kraus subcircuits.
@@ -110,7 +116,10 @@ class QuantumError {
   // Probabilities, first entry is no-error (identity)
   rvector_t probabilities_;
 
-  // List of unitary error matrices
+  // Generator rates if Pauli lindblad generated error
+  rvector_t rates_;
+
+  // List of unitary error matrices or Pauli generators
   std::vector<NoiseOps> circuits_;
 
   // List of OpTypes contained in error circuits
@@ -126,6 +135,29 @@ class QuantumError {
 
   // flag for where errors should be applied relative to the sampled op
   bool errors_after_op_ = true;
+
+  // flag for whether circuits_ represent channel generators or terms
+  bool generator_circuits_ = false;
+
+  // Sample gate level noise from circuits
+  NoiseOps sample_noise_circuits(const reg_t &qubits, RngEngine &rng) const;
+
+  // Sample gate level noise from generators
+  NoiseOps sample_noise_generators(const reg_t &qubits, RngEngine &rng) const;
+
+  // Helper method to convert a Pauli Op to a Pauli::Pauli
+  Pauli::Pauli<BV::BinaryVector> get_op_pauli(const Operations::Op &op) const;
+
+  // Helper method to convert a I,X,Y,Z or Pauli Op to a Pauli op
+  Operations::Op get_pauli_op(const Operations::Op &op) const;
+
+  // Compute the superoperator representation of the quantum error from channel
+  // circuits
+  void compute_circuits_superoperator();
+
+  // Compute the superoperator representation of the quantum error from
+  // generator circuits
+  void compute_generators_superoperator();
 };
 
 //-------------------------------------------------------------------------
@@ -157,24 +189,11 @@ QuantumError::NoiseOps QuantumError::sample_noise(const reg_t &qubits,
     return NoiseOps({op});
   }
   default: {
-    auto r = rng.rand_int(probabilities_);
-    // Check for invalid arguments
-    if (r + 1 > circuits_.size()) {
-      throw std::invalid_argument("QuantumError: probability outcome (" +
-                                  std::to_string(r) +
-                                  ")"
-                                  " is greater than number of circuits (" +
-                                  std::to_string(circuits_.size()) + ").");
-    }
-    NoiseOps noise_ops = circuits_[r];
-    // Add qubits to noise op commands;
-    for (auto &op : noise_ops) {
-      // Update qubits based on position in qubits list
-      for (auto &qubit : op.qubits) {
-        qubit = qubits[qubit];
-      }
+    if (generator_circuits_) {
+      return sample_noise_generators(qubits, rng);
+    } else {
+      return sample_noise_circuits(qubits, rng);
     }
-    return noise_ops;
   }
   }
 }
@@ -224,6 +243,47 @@ void QuantumError::set_circuits(const std::vector<NoiseOps> &circuits,
   set_num_qubits(num_qubits);
 }
 
+void QuantumError::set_generators(const std::vector<NoiseOps> &circuits,
+                                  const rvector_t &rates) {
+  if (rates.size() != circuits.size()) {
+    throw std::invalid_argument(
+        "QuantumError: invalid input, number of generator circuits (" +
+        std::to_string(circuits.size()) + ") and number of rates (" +
+        std::to_string(rates.size()) + ") are not equal.");
+  }
+  // Reset OpSet
+  generator_circuits_ = true;
+  opset_ = Operations::OpSet();
+  uint_t num_qubits = 0;
+
+  // Add elements with non-zero probability
+  for (size_t j = 0; j < rates.size(); j++) {
+    if (abs(rates[j]) > threshold_) {
+      if (rates[j] < 0) {
+        throw std::invalid_argument(
+            "QuantumError: cannot contain negative rates");
+      }
+      NoiseOps circuit;
+      for (const auto &op : circuits[j]) {
+        if (op.name == "i") {
+          // Ignore identities
+          continue;
+        }
+        auto pauli_op = get_pauli_op(op);
+        circuit.push_back(pauli_op);
+        opset_.insert(pauli_op);
+        // Check max qubit size
+        for (const auto &qubit : op.qubits) {
+          num_qubits = std::max(num_qubits, qubit + 1);
+        }
+      }
+      rates_.push_back(rates[j]);
+      circuits_.push_back(circuit);
+    }
+  }
+  set_num_qubits(num_qubits);
+}
+
 void QuantumError::set_from_kraus(const std::vector<cmatrix_t> &mats) {
   // Check input isn't empty
   if (mats.empty())
@@ -352,6 +412,85 @@ const std::vector<cmatrix_t> &QuantumError::kraus() const {
 }
 
 void QuantumError::compute_superoperator() {
+  if (generator_circuits_) {
+    compute_generators_superoperator();
+  } else {
+    compute_circuits_superoperator();
+  }
+}
+
+void QuantumError::compute_kraus() {
+  // Check superoperator representation is computed
+  if (superoperator_.empty()) {
+    compute_superoperator();
+  }
+  // Conver to Kraus
+  size_t dim = 1 << get_num_qubits();
+  canonical_kraus_ = Utils::superop2kraus(superoperator_, dim);
+}
+
+void QuantumError::load_from_json(const json_t &js) {
+  std::vector<NoiseOps> circuits;
+  JSON::get_value(circuits, "instructions", js);
+  if (JSON::check_key("rates", js)) {
+    rvector_t rates;
+    JSON::get_value(rates, "rates", js);
+    set_generators(circuits, rates);
+  } else {
+    rvector_t probs;
+    JSON::get_value(probs, "probabilities", js);
+    set_circuits(circuits, probs);
+  }
+}
+
+Pauli::Pauli<BV::BinaryVector>
+QuantumError::get_op_pauli(const Operations::Op &op) const {
+  // Initialize an empty Pauli
+  Pauli::Pauli<BV::BinaryVector> pauli(num_qubits_);
+  const auto size = op.qubits.size();
+  const auto pauli_str = op.string_params[0];
+  for (size_t i = 0; i < op.qubits.size(); ++i) {
+    const auto qubit = op.qubits[size - 1 - i];
+    switch (pauli_str[i]) {
+    case 'I':
+      break;
+    case 'X':
+      pauli.X.set1(qubit);
+      break;
+    case 'Y':
+      pauli.X.set1(qubit);
+      pauli.Z.set1(qubit);
+      break;
+    case 'Z':
+      pauli.Z.set1(qubit);
+      break;
+    default:
+      throw std::invalid_argument("invalid Pauli \'" +
+                                  std::to_string(pauli_str[i]) + "\'.");
+    }
+  }
+  return pauli;
+}
+
+Operations::Op QuantumError::get_pauli_op(const Operations::Op &op) const {
+  const std::string allowed_gates = "ixyz";
+  if (op.name == "pauli") {
+    return op;
+  }
+  if (op.name.size() == 1 &&
+      allowed_gates.find(op.name[0]) != std::string::npos) {
+    // Convert to a Pauli op
+    Operations::Op pauli_op;
+    pauli_op.type = Operations::OpType::gate;
+    pauli_op.name = "pauli";
+    pauli_op.qubits = op.qubits;
+    return pauli_op;
+  }
+  throw std::invalid_argument(
+      "QuantumError: generator errors can only contain Pauli ops");
+}
+
+void QuantumError::compute_circuits_superoperator() {
   // Initialize superoperator matrix to correct size
   size_t dim = 1ULL << (2 * get_num_qubits());
   superoperator_.initialize(dim, dim);
@@ -369,22 +508,85 @@ void QuantumError::compute_superoperator() {
   }
 }
 
-void QuantumError::compute_kraus() {
-  // Check superoperator representation is computed
-  if (superoperator_.empty()) {
-    compute_superoperator();
+void QuantumError::compute_generators_superoperator() {
+  // Initialize superoperator matrix to correct size
+  size_t dim = 1ULL << (2 * get_num_qubits());
+  const cmatrix_t iden = Linalg::Matrix::identity(dim);
+  superoperator_ = iden;
+
+  // We use the superoperator simulator state to do this
+  QubitSuperoperator::State<> superop;
+
+  for (size_t j = 0; j < circuits_.size(); j++) {
+    // Initialize identity superoperator
+    superop.initialize_qreg(get_num_qubits());
+    // Apply each gate in the circuit
+    // We don't need output data or RNG for this
+    ExperimentResult data;
+    RngEngine rng;
+    superop.apply_ops(circuits_[j].cbegin(), circuits_[j].cend(), data, rng);
+
+    // Now compute the superoperator for the generator term
+    auto p_iden = 0.5 + 0.5 * std::exp(-2 * rates_[j]);
+    const auto non_iden = superoperator_ * superop.move_to_matrix();
+    superoperator_ = p_iden * superoperator_ + (1 - p_iden) * non_iden;
   }
-  // Conver to Kraus
-  size_t dim = 1 << get_num_qubits();
-  canonical_kraus_ = Utils::superop2kraus(superoperator_, dim);
 }
 
-void QuantumError::load_from_json(const json_t &js) {
-  rvector_t probs;
-  JSON::get_value(probs, "probabilities", js);
-  std::vector<NoiseOps> circuits;
-  JSON::get_value(circuits, "instructions", js);
-  set_circuits(circuits, probs);
+QuantumError::NoiseOps
+QuantumError::sample_noise_circuits(const reg_t &qubits, RngEngine &rng) const {
+  auto r = rng.rand_int(probabilities_);
+  // Check for invalid arguments
+  if (r + 1 > circuits_.size()) {
+    throw std::invalid_argument("QuantumError: probability outcome (" +
+                                std::to_string(r) +
+                                ")"
+                                " is greater than number of circuits (" +
+                                std::to_string(circuits_.size()) + ").");
+  }
+  NoiseOps noise_ops = circuits_[r];
+  // Add qubits to noise op commands;
+  for (auto &op : noise_ops) {
+    // Update qubits based on position in qubits list
+    for (auto &qubit : op.qubits) {
+      qubit = qubits[qubit];
+    }
+  }
+  return noise_ops;
+}
+
+QuantumError::NoiseOps
+QuantumError::sample_noise_generators(const reg_t &qubits,
+                                      RngEngine &rng) const {
+  Pauli::Pauli<BV::BinaryVector> sampled_pauli(num_qubits_);
+  for (uint_t i = 0; i < rates_.size(); ++i) {
+    // Use the generator rate to calculate the probability of sampling the
+    // generator Pauli for the single term channel exp(r D[P]) = p S[I] + (1-p)
+    // * S[P] where p = 0.5 + 0.5 * exp(-2r)
+    auto p_iden = 0.5 + 0.5 * std::exp(-2 * rates_[i]);
+    if (rng.rand() > p_iden) {
+      // Convert each circuit Op to Pauli so we can more efficiently compose
+      // terms
+      for (auto op : circuits_[i]) {
+        auto term_pauli = get_op_pauli(op);
+        sampled_pauli.X += term_pauli.X;
+        sampled_pauli.Z += term_pauli.Z;
+      }
+    }
+  }
+
+  // If sampled Pauli is an identity we don't return error terms
+  if (sampled_pauli == Pauli::Pauli<BV::BinaryVector>(num_qubits_)) {
+    return NoiseOps();
+  }
+
+  // Convert back to an pauli Op for noise insertion
+  Operations::Op op;
+  op.type = Operations::OpType::gate;
+  op.name = "pauli";
+  op.qubits = qubits;
+  op.string_params = {sampled_pauli.str()};
+  return NoiseOps({op});
 }
 
 //-------------------------------------------------------------------------
diff --git a/src/simulators/circuit_executor.hpp b/src/simulators/circuit_executor.hpp
index 39ec6ce8fa..39f9a62cff 100644
--- a/src/simulators/circuit_executor.hpp
+++ b/src/simulators/circuit_executor.hpp
@@ -1155,7 +1155,7 @@ void Executor<state_t>::measure_sampler(InputIterator first_meas,
     Utils::apply_omp_parallel_for((npar > 1), 0, npar, copy_samples_lambda,
                                   npar);
 
-    for (int_t i = 0; i < npar; i++) {
+    for (uint_t i = 0; i < npar; i++) {
       result.combine(std::move(par_results[i]));
     }
   }
diff --git a/src/simulators/density_matrix/densitymatrix_executor.hpp b/src/simulators/density_matrix/densitymatrix_executor.hpp
index b6bdb67146..8a1081aa09 100644
--- a/src/simulators/density_matrix/densitymatrix_executor.hpp
+++ b/src/simulators/density_matrix/densitymatrix_executor.hpp
@@ -1329,10 +1329,10 @@ Executor<densmat_t>::sample_measure(const reg_t &qubits, uint_t shots,
   std::vector<SampleVector> all_samples(shots, SampleVector(qubits.size()));
 
   auto convert_to_bit_lambda = [this, &local_samples, &all_samples, shots,
-                                qubits, npar](int_t i) {
+                                qubits, npar](int_t k) {
     uint_t ishot, iend;
-    ishot = local_samples.size() * i / npar;
-    iend = local_samples.size() * (i + 1) / npar;
+    ishot = local_samples.size() * k / npar;
+    iend = local_samples.size() * (k + 1) / npar;
     for (; ishot < iend; ishot++) {
       SampleVector allbit_sample;
       allbit_sample.from_uint(local_samples[ishot], qubits.size());
diff --git a/src/simulators/density_matrix/densitymatrix_state.hpp b/src/simulators/density_matrix/densitymatrix_state.hpp
index cb54cb72a7..4300fed0db 100644
--- a/src/simulators/density_matrix/densitymatrix_state.hpp
+++ b/src/simulators/density_matrix/densitymatrix_state.hpp
@@ -1009,10 +1009,10 @@ std::vector<SampleVector> State<densmat_t>::sample_measure(const reg_t &qubits,
   std::vector<SampleVector> all_samples(shots, SampleVector(qubits.size()));
 
   auto convert_to_bit_lambda = [this, &allbit_samples, &all_samples, shots,
-                                qubits, npar](int_t i) {
+                                qubits, npar](int_t k) {
     uint_t ishot, iend;
-    ishot = shots * i / npar;
-    iend = shots * (i + 1) / npar;
+    ishot = shots * k / npar;
+    iend = shots * (k + 1) / npar;
     for (; ishot < iend; ishot++) {
       SampleVector allbit_sample;
       allbit_sample.from_uint(allbit_samples[ishot], qubits.size());
diff --git a/src/simulators/matrix_product_state/svd.cpp b/src/simulators/matrix_product_state/svd.cpp
index 150d97af4f..19c6e65cca 100644
--- a/src/simulators/matrix_product_state/svd.cpp
+++ b/src/simulators/matrix_product_state/svd.cpp
@@ -647,7 +647,7 @@ void lapack_csvd_wrapper(cmatrix_t &A, cmatrix_t &U, rvector_t &S,
   V = cmatrix_t::move_from_buffer(n, n, lapackV);
 
   S.clear();
-  for (int i = 0; i < min_dim; i++)
+  for (size_t i = 0; i < min_dim; i++)
     S.push_back(lapackS[i]);
 
   // Activated by default as requested in the PR
diff --git a/src/simulators/parallel_state_executor.hpp b/src/simulators/parallel_state_executor.hpp
index 0635da9ef5..030f92456a 100644
--- a/src/simulators/parallel_state_executor.hpp
+++ b/src/simulators/parallel_state_executor.hpp
@@ -706,7 +706,7 @@ void ParallelStateExecutor<state_t>::measure_sampler(InputIterator first_meas,
     Utils::apply_omp_parallel_for((npar > 1), 0, npar, copy_samples_lambda,
                                   npar);
 
-    for (int_t i = 0; i < npar; i++) {
+    for (uint_t i = 0; i < npar; i++) {
       result.combine(std::move(par_results[i]));
     }
   }
diff --git a/src/simulators/stabilizer/pauli.hpp b/src/simulators/stabilizer/pauli.hpp
index 39b3aae5fd..1a04f3c423 100644
--- a/src/simulators/stabilizer/pauli.hpp
+++ b/src/simulators/stabilizer/pauli.hpp
@@ -136,5 +136,11 @@ inline std::ostream &operator<<(std::ostream &out,
   return out;
 }
 
+template <typename bv_type>
+inline bool operator==(const AER::Pauli::Pauli<bv_type> &lhs,
+                       const AER::Pauli::Pauli<bv_type> &rhs) {
+  return lhs.X == rhs.X && lhs.Z == rhs.Z;
+}
+
 //------------------------------------------------------------------------------
 #endif
diff --git a/src/simulators/stabilizer/stabilizer_state.hpp b/src/simulators/stabilizer/stabilizer_state.hpp
index 12b80e981f..682ffee6ff 100644
--- a/src/simulators/stabilizer/stabilizer_state.hpp
+++ b/src/simulators/stabilizer/stabilizer_state.hpp
@@ -487,7 +487,7 @@ std::vector<SampleVector> State::sample_measure(const reg_t &qubits,
   // table
   auto qreg_cache = BaseState::qreg_;
   std::vector<SampleVector> samples(shots);
-  for (int_t ishot = 0; ishot < shots; ishot++) {
+  for (uint_t ishot = 0; ishot < shots; ishot++) {
     samples[ishot].from_vector(apply_measure_and_update(qubits, rng));
     BaseState::qreg_ = qreg_cache; // restore pre-measurement data from cache
   }
diff --git a/src/simulators/statevector/statevector_executor.hpp b/src/simulators/statevector/statevector_executor.hpp
index 26b00e801e..27193ff948 100644
--- a/src/simulators/statevector/statevector_executor.hpp
+++ b/src/simulators/statevector/statevector_executor.hpp
@@ -1247,10 +1247,10 @@ Executor<state_t>::sample_measure(const reg_t &qubits, uint_t shots,
   std::vector<SampleVector> all_samples(shots, SampleVector(qubits.size()));
 
   auto convert_to_bit_lambda = [this, &local_samples, &all_samples, shots,
-                                qubits, npar](int_t i) {
+                                qubits, npar](int_t k) {
     uint_t ishot, iend;
-    ishot = local_samples.size() * i / npar;
-    iend = local_samples.size() * (i + 1) / npar;
+    ishot = local_samples.size() * k / npar;
+    iend = local_samples.size() * (k + 1) / npar;
     for (; ishot < iend; ishot++) {
       SampleVector allbit_sample;
       allbit_sample.from_uint(local_samples[ishot], qubits.size());
diff --git a/src/simulators/statevector/statevector_state.hpp b/src/simulators/statevector/statevector_state.hpp
index 1981290256..7ea03ebdd3 100755
--- a/src/simulators/statevector/statevector_state.hpp
+++ b/src/simulators/statevector/statevector_state.hpp
@@ -1035,10 +1035,10 @@ std::vector<SampleVector> State<statevec_t>::sample_measure(const reg_t &qubits,
   std::vector<SampleVector> all_samples(shots, SampleVector(qubits.size()));
 
   auto convert_to_bit_lambda = [this, &allbit_samples, &all_samples, shots,
-                                qubits, npar](int_t i) {
+                                qubits, npar](int_t k) {
     uint_t ishot, iend;
-    ishot = shots * i / npar;
-    iend = shots * (i + 1) / npar;
+    ishot = shots * k / npar;
+    iend = shots * (k + 1) / npar;
     for (; ishot < iend; ishot++) {
       SampleVector allbit_sample;
       allbit_sample.from_uint(allbit_samples[ishot], qubits.size());
diff --git a/test/terra/backends/aer_simulator/test_noise.py b/test/terra/backends/aer_simulator/test_noise.py
index 3a7ebab34b..a4f223ae02 100644
--- a/test/terra/backends/aer_simulator/test_noise.py
+++ b/test/terra/backends/aer_simulator/test_noise.py
@@ -26,6 +26,7 @@
 from test.terra.backends.simulator_test_case import SimulatorTestCase, supported_methods
 from test.terra.reference import ref_kraus_noise
 from test.terra.reference import ref_pauli_noise
+from test.terra.reference import ref_pauli_lindblad_noise
 from test.terra.reference import ref_readout_noise
 from test.terra.reference import ref_reset_noise
 
@@ -140,6 +141,61 @@ def test_pauli_measure_noise(self, method, device, qerror_cls):
             self.assertSuccess(result)
             self.compare_counts(result, [circuit], [target], delta=0.05 * shots)
 
+    @supported_methods(ALL_METHODS)
+    def test_pauli_lindblad_gate_noise(self, method, device):
+        """Test simulation with Pauli gate error noise model."""
+        backend = self.backend(method=method, device=device)
+        shots = 1000
+        circuits = ref_pauli_lindblad_noise.pauli_lindblad_gate_error_circuits()
+        noise_models = ref_pauli_lindblad_noise.pauli_lindblad_gate_error_noise_models()
+        targets = ref_pauli_lindblad_noise.pauli_lindblad_gate_error_counts(shots)
+
+        for circuit, noise_model, target in zip(circuits, noise_models, targets):
+            backend.set_options(noise_model=noise_model)
+            result = backend.run(circuit, shots=shots).result()
+            self.assertSuccess(result)
+            self.compare_counts(result, [circuit], [target], delta=0.05 * shots)
+
+    @supported_methods(
+        [
+            "automatic",
+            "stabilizer",
+            "statevector",
+            "density_matrix",
+            "matrix_product_state",
+            "extended_stabilizer",
+            "tensor_network",
+        ]
+    )
+    def test_pauli_lindblad_reset_noise(self, method, device):
+        """Test simulation with Pauli reset error noise model."""
+        backend = self.backend(method=method, device=device)
+        shots = 1000
+        circuits = ref_pauli_lindblad_noise.pauli_lindblad_reset_error_circuits()
+        noise_models = ref_pauli_lindblad_noise.pauli_lindblad_reset_error_noise_models()
+        targets = ref_pauli_lindblad_noise.pauli_lindblad_reset_error_counts(shots)
+
+        for circuit, noise_model, target in zip(circuits, noise_models, targets):
+            backend.set_options(noise_model=noise_model)
+            result = backend.run(circuit, shots=shots).result()
+            self.assertSuccess(result)
+            self.compare_counts(result, [circuit], [target], delta=0.05 * shots)
+
+    @supported_methods(ALL_METHODS)
+    def test_pauli_lindblad_measure_noise(self, method, device):
+        """Test simulation with Pauli measure error noise model."""
+        backend = self.backend(method=method, device=device)
+        shots = 1000
+        circuits = ref_pauli_lindblad_noise.pauli_lindblad_measure_error_circuits()
+        noise_models = ref_pauli_lindblad_noise.pauli_lindblad_measure_error_noise_models()
+        targets = ref_pauli_lindblad_noise.pauli_lindblad_measure_error_counts(shots)
+
+        for circuit, noise_model, target in zip(circuits, noise_models, targets):
+            backend.set_options(noise_model=noise_model)
+            result = backend.run(circuit, shots=shots).result()
+            self.assertSuccess(result)
+            self.compare_counts(result, [circuit], [target], delta=0.05 * shots)
+
     @supported_methods(ALL_METHODS)
     def test_reset_gate_noise(self, method, device):
         """Test simulation with reset gate error noise model."""
diff --git a/test/terra/noise/test_noise_model.py b/test/terra/noise/test_noise_model.py
index 258f12579f..297e38d188 100644
--- a/test/terra/noise/test_noise_model.py
+++ b/test/terra/noise/test_noise_model.py
@@ -21,7 +21,7 @@
 from qiskit_aer.backends import AerSimulator
 from qiskit_aer.noise import NoiseModel
 from qiskit_aer.noise.device.models import _excited_population
-from qiskit_aer.noise.errors import PauliError
+from qiskit_aer.noise.errors import PauliError, PauliLindbladError
 from qiskit_aer.noise.errors.standard_errors import amplitude_damping_error
 from qiskit_aer.noise.errors.standard_errors import kraus_error
 from qiskit_aer.noise.errors.standard_errors import pauli_error
@@ -449,6 +449,77 @@ def test_pauli_error_equiv(self):
         self.assertTrue(result2.success)
         self.assertEqual(result1.get_counts(), result2.get_counts())
 
+    def test_pauli_lindblad_error_sampling(self):
+        num_qubits = 100
+        rate = 0.1
+        plerr = PauliLindbladError([num_qubits * "X"], [rate])
+        prob_err = 0.5 - 0.5 * np.exp(-2 * rate)
+
+        circ = QuantumCircuit(num_qubits)
+        circ.append(plerr, range(num_qubits))
+        circ.measure_all()
+
+        shots = 10000
+        backend = AerSimulator(method="stabilizer", seed_simulator=123)
+        result = backend.run(circ, shots=shots).result()
+        counts = result.get_counts()
+        self.assertEqual(len(counts), 2)
+        self.assertEqual(sum(counts.values()), shots)
+        self.assertIn(num_qubits * "0", counts)
+        self.assertIn(num_qubits * "1", counts)
+        np.testing.assert_allclose(counts[num_qubits * "1"] / shots, prob_err, atol=1e-3)
+
+    def test_pauli_lindblad_error_equiv(self):
+        circ = QuantumCircuit(2)
+        circ.h(0)
+        circ.cx(0, 1)
+        circ.measure_all()
+
+        perr1 = PauliLindbladError(["X", "Y", "Z"], [0.1, 0.2, 0.3])
+        perr2 = PauliLindbladError(["XX", "YY", "ZZ"], [0.1, 0.2, 0.3])
+
+        noise_model1 = NoiseModel()
+        noise_model1.add_all_qubit_quantum_error(perr1, ["h"])
+        noise_model1.add_all_qubit_quantum_error(perr2, ["cx"])
+
+        noise_model2 = NoiseModel()
+        noise_model2.add_all_qubit_quantum_error(perr1.to_quantum_error(), ["h"])
+        noise_model2.add_all_qubit_quantum_error(perr2.to_quantum_error(), ["cx"])
+
+        seed = 1234
+        result1 = (
+            AerSimulator()
+            .run(circ, shots=5000, noise_model=noise_model1, seed_simulator=seed)
+            .result()
+        )
+        result2 = (
+            AerSimulator()
+            .run(circ, shots=5000, noise_model=noise_model2, seed_simulator=seed)
+            .result()
+        )
+        self.assertTrue(result1.success)
+        self.assertTrue(result2.success)
+        self.assertEqual(result1.get_counts(), result2.get_counts())
+
+    def test_pauli_lindblad_error_sampling_equiv(self):
+        plerr = PauliLindbladError(["IX", "XI"], [0.5, 0.1])
+        circ1 = QuantumCircuit(2)
+        circ1.append(plerr, range(2))
+        circ1.measure_all()
+
+        circ2 = QuantumCircuit(2)
+        circ2.append(plerr.to_quantum_error(), range(2))
+        circ2.measure_all()
+
+        shots = 100000
+        backend = AerSimulator(seed_simulator=123)
+        result = backend.run([circ1, circ2], shots=shots).result()
+        counts1 = result.get_counts(0)
+        counts2 = result.get_counts(1)
+        probs1 = [counts1.get(i, 0) / shots for i in ["00", "01", "10", "11"]]
+        probs2 = [counts2.get(i, 0) / shots for i in ["00", "01", "10", "11"]]
+        np.testing.assert_allclose(probs1, probs2, atol=5e-2)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/terra/noise/test_pauli_lindblad_error.py b/test/terra/noise/test_pauli_lindblad_error.py
new file mode 100644
index 0000000000..1a2987870d
--- /dev/null
+++ b/test/terra/noise/test_pauli_lindblad_error.py
@@ -0,0 +1,332 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018-2024.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+PauliLindbladError class tests
+"""
+from test.terra.common import QiskitAerTestCase
+
+import unittest
+from itertools import combinations, product
+
+import ddt
+import numpy as np
+import scipy.linalg as la
+
+import qiskit.quantum_info as qi
+from qiskit_aer.noise import PauliLindbladError, PauliError, QuantumError
+from qiskit_aer.noise.noiseerror import NoiseError
+
+
+@ddt.ddt
+class TestPauliLindbladError(QiskitAerTestCase):
+    """Testing PauliLindbladError class"""
+
+    @ddt.data(str, qi.Pauli, qi.PauliList)
+    def test_init_with_paulis(self, cls):
+        """Test construction with different inputs."""
+        gens = ["I", "X"]
+        if cls == qi.Pauli:
+            gens = [qi.Pauli(i) for i in gens]
+        elif cls == qi.PauliList:
+            gens = qi.PauliList(gens)
+        perr = PauliLindbladError(gens, np.ones(len(gens)) / len(gens))
+        self.assertEqual(perr.size, 2)
+
+    def test_invalid_inputs(self):
+        """Test inputs with different lengths raises"""
+        with self.assertRaises(NoiseError):
+            PauliLindbladError(["X", "Z"], [1])
+        with self.assertRaises(NoiseError):
+            PauliLindbladError(["X", "Z"], np.eye(2))
+        with self.assertRaises(NoiseError):
+            PauliLindbladError(["X"], 1)
+
+    def test_negative_rates(self):
+        """Test initalizing with negative rate."""
+        perr = PauliLindbladError(["X", "Z"], [0.1, -0.1])
+        self.assertTrue(perr.is_tp())
+        self.assertFalse(perr.is_cp())
+        self.assertFalse(perr.is_cptp())
+        with self.assertRaises(NoiseError):
+            perr.to_quantum_error()
+
+    def test_attributes(self):
+        """Test basic attributes"""
+        rng = np.random.default_rng(555)
+        gens = qi.random_pauli_list(3, 8, phase=False, seed=rng)
+        rates = rng.random(len(gens))
+        perr = PauliLindbladError(gens, rates)
+        self.assertEqual(perr.size, len(gens))
+        self.assertEqual(perr.generators, gens)
+        np.testing.assert_allclose(perr.rates, rates)
+
+    @ddt.data(1, 10, 100, 1000)
+    def test_ideal_single_iden(self, num_qubits):
+        """Test ideal gates are identified correctly."""
+        self.assertTrue(PauliLindbladError([num_qubits * "I"], [0.1]).ideal())
+        self.assertTrue(PauliLindbladError([num_qubits * "I"], [-0.1]).ideal())  # non-CP
+
+    @ddt.data(1, 10, 100, 1000)
+    def test_ideal_single_zero_rates(self, num_qubits):
+        """Test ideal gates are identified correctly."""
+        gens = [num_qubits * s for s in ["X", "Y", "Z"]]
+        rates = [0, 0, 0]
+        self.assertTrue(PauliLindbladError(gens, rates).ideal())
+
+    @ddt.data(1, 2, 3)
+    def test_to_quantum_channel_iden(self, num_qubits):
+        """Test ideal gates are identified correctly."""
+        gens = [num_qubits * s for s in ["X", "Y", "Z"]]
+        rates = [0, 0, 0]
+        plerr = PauliLindbladError(gens, rates)
+        chan = plerr.to_quantumchannel()
+        self.assertEqual(chan, qi.SuperOp(np.eye(4**num_qubits)))
+
+    @ddt.idata(product(["pauli_error", "quantum_error", "quantumchannel"], [0, 0.1, 0.25, 0.5]))
+    @ddt.unpack
+    def test_conversion_1q_single(self, output, prob):
+        """Test conversion into quantum channel."""
+        if prob == 0.5:
+            rate = np.inf
+        else:
+            rate = -0.5 * np.log(1 - 2 * prob)
+        plerr = PauliLindbladError(["Y"], [rate])
+        target = PauliError(["I", "Y"], [1 - prob, prob])
+        if output == "pauli_error":
+            value = plerr.to_pauli_error()
+        elif output == "quantum_error":
+            value = plerr.to_quantum_error()
+            target = target.to_quantum_error()
+        elif output == "quantumchannel":
+            value = plerr.to_quantumchannel()
+            target = target.to_quantumchannel()
+        self.assertEqual(value, target)
+
+    @ddt.idata(product(["pauli_error", "quantum_error", "quantumchannel"], [0, 0.1, 0.25, 0.5]))
+    @ddt.unpack
+    def test_conversion_1q_multi(self, output, prob):
+        """Test conversion into quantum channel."""
+        if prob == 0.5:
+            rate = np.inf
+        else:
+            rate = -0.5 * np.log(1 - 2 * prob)
+        plerr = PauliLindbladError(["X", "Y", "Z"], 3 * [rate])
+        plerr = plerr.simplify()
+        target = PauliError(["I"], [1])
+        for pauli in ["X", "Y", "Z"]:
+            target = target.compose(PauliError(["I", pauli], [1 - prob, prob]))
+        target = target.simplify()
+        if output == "pauli_error":
+            value = plerr.to_pauli_error()
+        elif output == "quantum_error":
+            value = plerr.to_quantum_error()
+            target = target.to_quantum_error()
+        elif output == "quantumchannel":
+            value = plerr.to_quantumchannel()
+            target = target.to_quantumchannel()
+        self.assertEqual(value, target)
+
+    @ddt.data(1, 2, 3, 4)
+    def test_to_quantumchannel_rand(self, num_qubits):
+        """Test conversion into quantum channel."""
+        rng = np.random.default_rng(55)
+        gens = qi.random_pauli_list(num_qubits, 3 * num_qubits, phase=False, seed=rng)
+        rates = rng.random(len(gens))
+        plerr = PauliLindbladError(gens, rates)
+        target = plerr.to_pauli_error().to_quantumchannel()
+        self.assertEqual(plerr.to_quantumchannel(), target)
+
+    def test_dict_round_trip(self):
+        """Test to_dict and from_dict round trip."""
+        rng = np.random.default_rng(55)
+        gens = qi.random_pauli_list(5, 6, phase=False, seed=rng)
+        rates = rng.random(len(gens))
+        target = PauliLindbladError(gens, rates)
+        value = PauliLindbladError.from_dict(target.to_dict())
+        self.assertDictAlmostEqual(value, target)
+
+    def test_simplify_zero_rates(self):
+        """Test simplify removes zero rates"""
+        gens = ["XI", "IX", "YI", "IY"]
+        rates = [1e-1, 1e-2, 1e-5, 1e-9]
+        perr = PauliLindbladError(gens, rates)
+
+        value1 = perr.simplify()
+        target1 = PauliLindbladError(gens[:3], rates[:3])
+        self.assertEqual(value1, target1)
+
+        value2 = perr.simplify(atol=1e-4)
+        target2 = PauliLindbladError(gens[:2], rates[:2])
+        self.assertEqual(value2, target2)
+
+        value3 = perr.simplify(atol=0)
+        target3 = perr
+        self.assertEqual(value3, target3)
+
+    def test_simplify_duplicates(self):
+        """Test simplify combines duplicate terms"""
+        gens = ["XX", "ZZ"]
+        rates = [0.1, 0.2]
+        target = PauliLindbladError(gens, rates)
+        value = PauliLindbladError(4 * gens, np.array(4 * rates) / 4).simplify()
+        self.assertEqual(value, target)
+
+    def test_equal_diff_order(self):
+        rng = np.random.default_rng(33)
+        gens = qi.random_pauli_list(5, 10, phase=False, seed=rng)
+        rates = rng.random(len(gens))
+        perr1 = PauliLindbladError(gens, rates)
+        perr2 = PauliLindbladError(gens[::-1], rates[::-1])
+        self.assertEqual(perr1, perr2)
+
+    def test_not_equal_type(self):
+        perr = PauliLindbladError(["ZZ", "XX"], [0.2, 0.1])
+        qerr = perr.to_quantum_error()
+        chan = perr.to_quantumchannel()
+        self.assertNotEqual(perr, qerr)
+        self.assertNotEqual(perr, chan)
+
+    def test_not_equal_shape(self):
+        perr1 = PauliLindbladError(["YY", "XX"], [0.9, 0.1])
+        perr2 = PauliLindbladError(["YY", "XX", "ZZ"], [0.9, 0.05, 0.05])
+        self.assertNotEqual(perr1, perr2)
+
+    @ddt.data(1, 2, 3, 4)
+    def test_compose(self, qarg_qubits):
+        """Test compose two quantum errors."""
+        rng = np.random.default_rng(33)
+        gens1 = qi.random_pauli_list(4, 3, phase=False, seed=rng)
+        rates1 = rng.random(len(gens1))
+        perr1 = PauliLindbladError(gens1, rates1)
+        gens2 = qi.random_pauli_list(qarg_qubits, 3, phase=False, seed=rng)
+        rates2 = rng.random(len(gens2))
+        perr2 = PauliLindbladError(gens2, rates2)
+        target = perr1.to_quantumchannel().compose(perr2.to_quantumchannel(), range(qarg_qubits))
+        value = (perr1.compose(perr2, range(qarg_qubits))).to_quantumchannel()
+        self.assertEqual(value, target)
+
+    @ddt.idata(list(combinations(range(4), 1)) + list(combinations(range(4), 2)))
+    def test_compose_subsystem(self, qargs):
+        """Test compose with 1 and 2-qubit subsystem permutations"""
+        rng = np.random.default_rng(123)
+        gens1 = qi.random_pauli_list(4, 3, phase=False, seed=rng)
+        rates1 = rng.random(len(gens1))
+        perr1 = PauliLindbladError(gens1, rates1)
+        gens2 = qi.random_pauli_list(len(qargs), 3, phase=False, seed=rng)
+        rates2 = rng.random(len(gens2))
+        perr2 = PauliLindbladError(gens2, rates2)
+        target = perr1.to_quantumchannel().compose(perr2.to_quantumchannel(), qargs)
+        value = (perr1.compose(perr2, qargs)).to_quantumchannel()
+        self.assertEqual(value, target)
+
+    @ddt.data(1, 10, 100)
+    def test_dot_equals_compose(self, num_qubits):
+        """Test dot is equal to compose."""
+        rng = np.random.default_rng(999)
+
+        for _ in range(20):
+            gens1 = qi.random_pauli_list(num_qubits, 3, phase=False, seed=rng)
+            coeffs1 = rng.random(len(gens1))
+            perr1 = PauliLindbladError(gens1, coeffs1)
+            gens2 = qi.random_pauli_list(num_qubits, 5, phase=False, seed=rng)
+            coeffs2 = rng.random(len(gens2))
+            perr2 = PauliLindbladError(gens2, coeffs2)
+            self.assertEqual(perr1.dot(perr2), perr1.compose(perr2))
+
+    def test_tensor(self):
+        """Test tensor two quantum errors."""
+        rng = np.random.default_rng(99)
+        gens1 = qi.random_pauli_list(2, 4, phase=False, seed=rng)
+        rates1 = rng.random(len(gens1))
+        rates1 /= sum(rates1)
+        perr1 = PauliLindbladError(gens1, rates1)
+        gens2 = qi.random_pauli_list(2, 3, phase=False, seed=rng)
+        rates2 = rng.random(len(gens2))
+        rates2 /= sum(rates2)
+        perr2 = PauliLindbladError(gens2, rates2)
+        value = perr1.tensor(perr2).to_quantumchannel()
+        target = perr1.to_quantumchannel().tensor(perr2.to_quantumchannel())
+        self.assertEqual(value, target)
+
+    def test_expand(self):
+        """Test tensor two quantum errors."""
+        rng = np.random.default_rng(99)
+        gens1 = qi.random_pauli_list(2, 4, phase=False, seed=rng)
+        rates1 = rng.random(len(gens1))
+        perr1 = PauliLindbladError(gens1, rates1)
+        gens2 = qi.random_pauli_list(2, 3, phase=False, seed=rng)
+        rates2 = rng.random(len(gens2))
+        perr2 = PauliLindbladError(gens2, rates2)
+        value = perr1.expand(perr2).to_quantumchannel()
+        target = perr1.to_quantumchannel().expand(perr2.to_quantumchannel())
+        self.assertEqual(value, target)
+
+    @ddt.data(0, 1, -1, 0.5, 2)
+    def test_power(self, exponent):
+        """Test power method"""
+        rng = np.random.default_rng(seed=1234)
+        gens = qi.random_pauli_list(3, 5, phase=False)
+        rates = rng.random(5)
+        plerr = PauliLindbladError(gens, rates)
+        plpow = plerr.power(exponent)
+        supmat = qi.SuperOp(plerr).data
+        if exponent == 0.5:
+            suppow = la.sqrtm(supmat)
+        elif exponent == -1:
+            suppow = la.inv(supmat)
+        else:
+            suppow = la.fractional_matrix_power(supmat, exponent)
+        supop_pow = qi.SuperOp(suppow)
+        self.assertEqual(plpow, PauliLindbladError(gens, exponent * rates))
+        self.assertEqual(plpow.to_quantumchannel(), supop_pow)
+
+    def test_inverse(self):
+        """Test inverse"""
+        rng = np.random.default_rng(seed=1234)
+        gens = qi.random_pauli_list(3, 5, phase=False)
+        rates = rng.random(5)
+        plerr = PauliLindbladError(gens, rates)
+        inv = plerr.inverse()
+        supop = qi.SuperOp(plerr)
+        supop_inv = qi.SuperOp(np.linalg.inv(supop.data))
+        self.assertEqual(inv, PauliLindbladError(gens, -rates))
+        self.assertEqual(inv.to_quantumchannel(), supop_inv)
+
+    def test_subsystem_errors(self):
+        """Test subsystem errors method"""
+        rng = np.random.default_rng(seed=1234)
+        gens = qi.PauliList(["XX", "YY", "ZZ"])
+        plerr1 = PauliLindbladError(gens, rng.random(3))
+        plerr2 = PauliLindbladError(gens, rng.random(3))
+        plerr3 = PauliLindbladError(gens, rng.random(3))
+        plerr4 = PauliLindbladError(gens, rng.random(3))
+
+        error = plerr1.expand(PauliLindbladError(["III"], [0]))
+        error = error.compose(plerr2, [1, 2])
+        error = error.compose(plerr3, [2, 3])
+        error = error.compose(plerr4, [3, 4])
+
+        targets = {
+            (0, 1): plerr1,
+            (1, 2): plerr2,
+            (2, 3): plerr3,
+            (3, 4): plerr4,
+        }
+        sub_errors = error.subsystem_errors()
+        self.assertEqual(len(sub_errors), 4)
+        for suberr, subsys in sub_errors:
+            self.assertIn(subsys, targets)
+            self.assertEqual(suberr, targets[subsys])
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/terra/reference/ref_pauli_lindblad_noise.py b/test/terra/reference/ref_pauli_lindblad_noise.py
new file mode 100644
index 0000000000..83a05dc05b
--- /dev/null
+++ b/test/terra/reference/ref_pauli_lindblad_noise.py
@@ -0,0 +1,256 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018, 2019.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+AerSimulator readout error NoiseModel integration tests
+"""
+from math import log, inf
+from test.terra.utils.utils import list2dict
+
+from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
+from qiskit_aer.noise import NoiseModel
+from qiskit_aer.noise.errors import PauliLindbladError
+
+
+# ==========================================================================
+# Pauli Gate Errors
+# ==========================================================================
+
+
+def _pauli_rate(prob_error):
+    """Convert single error probability to generator rate"""
+    # NOTE: PauliLindbladError cannot produce a Pauli error with error prob > 50%
+    if prob_error > 0.5:
+        raise ValueError("Error probability cannot be > 50%")
+    if prob_error == 0.5:
+        return inf
+    if prob_error == 0:
+        return 0.0
+    return -0.5 * log(1 - 2 * prob_error)
+
+
+def pauli_lindblad_gate_error_circuits():
+    """Local PauliLindbladError gate error noise model circuits"""
+    circuits = []
+
+    qr = QuantumRegister(2, "qr")
+    cr = ClassicalRegister(2, "cr")
+
+    # 25% all-qubit Pauli error on "id" gates
+    circuit = QuantumCircuit(qr, cr)
+    circuit.id(qr)
+    circuit.barrier(qr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    # 25% all-qubit Pauli error on "id" gates on qubit-0
+    circuit = QuantumCircuit(qr, cr)
+    circuit.id(qr)
+    circuit.barrier(qr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    # 50% Pauli error on conditional gate that doesn't get applied
+    circuit = QuantumCircuit(qr, cr)
+    circuit.x(qr).c_if(cr, 1)
+    circuit.barrier(qr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    # 50% Pauli error on conditional gate that does get applied
+    circuit = QuantumCircuit(qr, cr)
+    circuit.x(qr).c_if(cr, 0)
+    circuit.barrier(qr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    return circuits
+
+
+def pauli_lindblad_gate_error_noise_models():
+    """Local Pauli gate error noise models"""
+    noise_models = []
+
+    # 25% all-qubit Pauli error on "id" gates
+    error = PauliLindbladError(["X"], [_pauli_rate(0.25)])
+    noise_model = NoiseModel()
+    noise_model.add_all_qubit_quantum_error(error, "id")
+    noise_models.append(noise_model)
+
+    # 25% all-qubit Pauli error on "id" gates on qubit-0
+    error = PauliLindbladError(["X"], [_pauli_rate(0.25)])
+    noise_model = NoiseModel()
+    noise_model.add_quantum_error(error, "id", [0])
+    noise_models.append(noise_model)
+
+    # 50% Pauli error on conditional gate that doesn't get applied
+    error = PauliLindbladError(["X"], [_pauli_rate(0.5)])
+    noise_model = NoiseModel()
+    noise_model.add_all_qubit_quantum_error(error, "x")
+    noise_models.append(noise_model)
+
+    # 50% Pauli error on conditional gate that does get applied
+    error = PauliLindbladError(["X"], [_pauli_rate(0.5)])
+    noise_model = NoiseModel()
+    noise_model.add_all_qubit_quantum_error(error, "x")
+    noise_models.append(noise_model)
+
+    return noise_models
+
+
+def pauli_lindblad_gate_error_counts(shots, hex_counts=True):
+    """Pauli gate error circuits reference counts"""
+    counts_lists = []
+
+    # 25% all-qubit Pauli error on "id" gates
+    counts = [9 * shots / 16, 3 * shots / 16, 3 * shots / 16, shots / 16]
+    counts_lists.append(counts)
+
+    # 25% all-qubit Pauli error on "id" gates on qubit-0
+    counts = [3 * shots / 4, shots / 4, 0, 0]
+    counts_lists.append(counts)
+
+    # 50% Pauli error on conditional gate that doesn't get applied
+    counts = [shots, 0, 0, 0]
+    counts_lists.append(counts)
+
+    # 50% Pauli error on conditional gate that does get applied
+    counts = 4 * [shots / 4]
+    counts_lists.append(counts)
+
+    return [list2dict(i, hex_counts) for i in counts_lists]
+
+
+# ==========================================================================
+# Pauli Measure Errors
+# ==========================================================================
+
+
+def pauli_lindblad_measure_error_circuits():
+    """Local Pauli measure error noise model circuits"""
+    circuits = []
+
+    qr = QuantumRegister(2, "qr")
+    cr = ClassicalRegister(2, "cr")
+
+    # 25% all-qubit Pauli error on measure
+    circuit = QuantumCircuit(qr, cr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    # 25% local Pauli error on measure of qubit 1
+    circuit = QuantumCircuit(qr, cr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    return circuits
+
+
+def pauli_lindblad_measure_error_noise_models():
+    """Local Pauli measure error noise models"""
+    noise_models = []
+
+    # 25% all-qubit Pauli error on measure
+    error = PauliLindbladError(["X"], [_pauli_rate(0.25)])
+    noise_model = NoiseModel()
+    noise_model.add_all_qubit_quantum_error(error, "measure")
+    noise_models.append(noise_model)
+
+    # 25% local Pauli error on measure of qubit 1
+    error = PauliLindbladError(["X"], [_pauli_rate(0.25)])
+    noise_model = NoiseModel()
+    noise_model.add_quantum_error(error, "measure", [1])
+    noise_models.append(noise_model)
+
+    return noise_models
+
+
+def pauli_lindblad_measure_error_counts(shots, hex_counts=True):
+    """Local Pauli measure error circuits reference counts"""
+    counts_lists = []
+
+    # 25% all-qubit Pauli error on measure
+    counts = [9 * shots / 16, 3 * shots / 16, 3 * shots / 16, shots / 16]
+    counts_lists.append(counts)
+
+    # 25% local Pauli error on measure of qubit 1
+    counts = [3 * shots / 4, 0, shots / 4, 0]
+    counts_lists.append(counts)
+
+    # Convert to counts dict
+    return [list2dict(i, hex_counts) for i in counts_lists]
+
+
+# ==========================================================================
+# Pauli Reset Errors
+# ==========================================================================
+
+
+def pauli_lindblad_reset_error_circuits():
+    """Local Pauli reset error noise model circuits"""
+    circuits = []
+
+    qr = QuantumRegister(2, "qr")
+    cr = ClassicalRegister(2, "cr")
+
+    # 25% all-qubit Pauli error on reset
+    circuit = QuantumCircuit(qr, cr)
+    circuit.barrier(qr)
+    circuit.reset(qr)
+    circuit.barrier(qr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    # 25% local Pauli error on reset of qubit 1
+    circuit = QuantumCircuit(qr, cr)
+    circuit.barrier(qr)
+    circuit.reset(qr)
+    circuit.barrier(qr)
+    circuit.measure(qr, cr)
+    circuits.append(circuit)
+
+    return circuits
+
+
+def pauli_lindblad_reset_error_noise_models():
+    """Local Pauli reset error noise models"""
+    noise_models = []
+
+    # 25% all-qubit Pauli error on reset
+    error = PauliLindbladError(["X"], [_pauli_rate(0.25)])
+    noise_model = NoiseModel()
+    noise_model.add_all_qubit_quantum_error(error, "reset")
+    noise_models.append(noise_model)
+
+    # 25% local Pauli error on reset of qubit 1
+    error = PauliLindbladError(["X"], [_pauli_rate(0.25)])
+    noise_model = NoiseModel()
+    noise_model.add_quantum_error(error, "reset", [1])
+    noise_models.append(noise_model)
+
+    return noise_models
+
+
+def pauli_lindblad_reset_error_counts(shots, hex_counts=True):
+    """Local Pauli reset error circuits reference counts"""
+    counts_lists = []
+
+    # 25% all-qubit Pauli error on reset
+    counts = [9 * shots / 16, 3 * shots / 16, 3 * shots / 16, shots / 16]
+    counts_lists.append(counts)
+
+    # 25% local Pauli error on reset of qubit 1
+    counts = [3 * shots / 4, 0, shots / 4, 0]
+    counts_lists.append(counts)
+
+    # Convert to counts dict
+    return [list2dict(i, hex_counts) for i in counts_lists]