diff --git a/.pylintrc b/.pylintrc
index 8562476456..37e7904a5c 100644
--- a/.pylintrc
+++ b/.pylintrc
@@ -119,7 +119,7 @@ evaluation=10.0 - ((float(5 * error + warning + refactor + convention) / stateme
 # op = operation iterator
 # b = basis iterator
 good-names=a,b,i,j,k,d,n,m,ex,v,w,x,y,z,Run,_,logger,q,c,r,qr,cr,qc,nd,pi,op,b,ar,br,p,cp,dt,
-    __unittest,iSwapGate,mu
+    __unittest,iSwapGate,mu,t1,t2
 
 # Bad variable names which should always be refused, separated by a comma
 bad-names=foo,bar,toto,tutu,tata
diff --git a/qiskit_experiments/library/randomized_benchmarking/rb_utils.py b/qiskit_experiments/library/randomized_benchmarking/rb_utils.py
index db52b4a554..e60ec24e20 100644
--- a/qiskit_experiments/library/randomized_benchmarking/rb_utils.py
+++ b/qiskit_experiments/library/randomized_benchmarking/rb_utils.py
@@ -13,13 +13,16 @@
 """
 RB Helper functions
 """
+import functools
+import operator
 from typing import Tuple, Dict, Optional, List, Union, Sequence
 
 import numpy as np
 import uncertainties
+
+import qiskit.quantum_info as qi
 from qiskit import QiskitError, QuantumCircuit
 from qiskit.providers.backend import Backend
-
 from qiskit_experiments.database_service.device_component import Qubit
 from qiskit_experiments.framework import DbAnalysisResultV1, AnalysisResultData
 from qiskit_experiments.warnings import deprecated_function
@@ -138,6 +141,10 @@ def gates_per_clifford(
         return {key: value[0] / value[1] for (key, value) in result.items()}
 
     @staticmethod
+    @deprecated_function(
+        last_version="0.4",
+        msg="Please use coherence_limit_error function that can handle three or more qubits instead.",
+    )
     def coherence_limit(nQ=2, T1_list=None, T2_list=None, gatelen=0.1):
         """
         The error per gate (1-average_gate_fidelity) given by the T1,T2 limit.
@@ -199,6 +206,83 @@ def coherence_limit(nQ=2, T1_list=None, T2_list=None, gatelen=0.1):
 
         return coherence_limit_err
 
+    @staticmethod
+    def coherence_limit_error(
+        num_qubits: int, gate_length: float, t1s: Sequence, t2s: Optional[Sequence] = None
+    ):
+        r"""
+        The error per gate (1 - average_gate_fidelity) given by the T1,T2 limit
+        assuming qubit-wise gate-independent amplitude damping error
+        (i.e. thermal relaxation error with no excitation).
+
+        That means, suppose the gate length $t$, we are considering a quantum error channel
+        whose Choi matrix representation for a single qubit with $T_1$ and $T_2$ is give by
+
+        .. math::
+
+            \begin{bmatrix}
+                1 & 0 & 0 & e^{-\frac{t}{T_2}} \\
+                0 & 0 & 0 & 0 \\
+                0 & 0 & 1-e^{-\frac{t}{T_1}} & 0 \\
+                e^{-\frac{t}{T_2}} & 0 & 0 & e^{-\frac{t}{T_1}} \\
+            \end{bmatrix}
+
+        The coherence limit error computed by this function is
+        :math:`1 - F_{\text{avg}}(\mathcal{E}, U)` and the following equalities hold for the value.
+
+        .. math::
+
+             \begin{align}
+             1 - F_{\text{avg}}(\mathcal{E}, U)
+                    &= \frac{d}{d+1} \left(1 - F_{\text{pro}}(\mathcal{E}, U)\right) \\
+                    &= \frac{d}{d+1} \left(1 - \frac{Tr[S_U^\dagger S_{\mathcal{E}}]}{d^2}\right) \\
+                    &= \frac{d}{d+1} \left(1 - \frac{Tr[S_{\Lambda}]}{d^2}\right)
+             \end{align}
+
+        where :math:`F_{\text{avg}}(\mathcal{E}, U)` and :math:`F_{\text{pro}}(\mathcal{E}, U)` are
+        the average gate fidelity and the process fidelity of a quantum channel :math:`\mathcal{E}`
+        with a target unitary $U$ such that :math:`\mathcal{E}=\Lambda(U)`, respectively,
+        $d$ is the dimension of Hilbert space of the considering qubit system, and
+        :math:`S_{\Lambda}` is the Liouville Superoperator representation of a channel :math:`\Lambda`.
+
+        Args:
+            num_qubits: Number of qubits.
+            gate_length: Duration of the gate in seconds.
+            t1s: List of T1's from qubit 0 to num_qubits-1.
+            t2s: List of T2's (as measured, not Tphi). If not given, assume T2 = 2 * T1.
+                 Each T2 value is truncated down to 2 * T1 if T2 > 2 * T1.
+
+        Returns:
+            float: coherence limited error per gate.
+        Raises:
+            ValueError: if there are invalid inputs
+        """
+        t1s = np.array(t1s)
+        if t2s is None:
+            t2s = 2 * t1s
+        else:
+            t2s = np.array([min(t2, 2 * t1) for t1, t2 in zip(t1s, t2s)])
+
+        if len(t1s) != num_qubits or len(t2s) != num_qubits:
+            raise ValueError("Length of t1s/t2s must equal num_qubits")
+
+        def amplitude_damping_choi(t1, t2, time):
+            return qi.Choi(
+                np.array(
+                    [
+                        [1, 0, 0, np.exp(-time / t2)],
+                        [0, 0, 0, 0],
+                        [0, 0, 1 - np.exp(-time / t1), 0],
+                        [np.exp(-time / t2), 0, 0, np.exp(-time / t1)],
+                    ]
+                )
+            )
+
+        chois = [amplitude_damping_choi(t1, t2, gate_length) for t1, t2 in zip(t1s, t2s)]
+        traces = [np.real(np.trace(np.array(qi.SuperOp(choi)))) for choi in chois]
+        d = 2**num_qubits
+        return d / (d + 1) * (1 - functools.reduce(operator.mul, traces) / (d * d))
+
     @staticmethod
     @deprecated_function(
         last_version="0.4",
diff --git a/releasenotes/notes/upgrade-coherence-limit-9a1fccac402b4f5c.yaml b/releasenotes/notes/upgrade-coherence-limit-9a1fccac402b4f5c.yaml
new file mode 100644
index 0000000000..cd55b2ce6d
--- /dev/null
+++ b/releasenotes/notes/upgrade-coherence-limit-9a1fccac402b4f5c.yaml
@@ -0,0 +1,15 @@
+---
+upgrade:
+  - |
+    A function to compute the coherence limit error, :meth:`.RBUtils.coherence_limit_error`,
+    is added and now it can take any number of qubits. It replaces the deprecated
+    :meth:`.RBUtils.coherence_limit`, which can take only one or two qubits.
+    Note that the names and order of the arguments have also changed.
+
+    .. code-block:: python
+
+      # New function
+      RBUtils.coherence_limit_error(num_qubits=num_qubits, gate_length=gate_length, t1s=t1s, t2s=t2s)
+
+      # Deprecated method
+      # RBUtils.coherence_limit(nQ=num_qubits, T1_list=t1s, T2_list=t2s, gatelen=gate_length)
diff --git a/test/randomized_benchmarking/test_rb_utils.py b/test/randomized_benchmarking/test_rb_utils.py
index c9163e287d..0c2290226d 100644
--- a/test/randomized_benchmarking/test_rb_utils.py
+++ b/test/randomized_benchmarking/test_rb_utils.py
@@ -167,14 +167,52 @@ def test_coherence_limit(self):
         t2 = 100.0
         gate_2_qubits = 0.5
         gate_1_qubit = 0.1
-        twoq_coherence_err = rb.RBUtils.coherence_limit(2, [t1, t1], [t2, t2], gate_2_qubits)
 
-        oneq_coherence_err = rb.RBUtils.coherence_limit(1, [t1], [t2], gate_1_qubit)
+        with self.assertWarns(DeprecationWarning):
+            oneq_coherence_err = rb.RBUtils.coherence_limit(1, [t1], [t2], gate_1_qubit)
+            twoq_coherence_err = rb.RBUtils.coherence_limit(2, [t1, t1], [t2, t2], gate_2_qubits)
+            # random test to ensure ole and new coherence_limit yield the same value for 1q and 2q cases
+            import random
 
-        self.assertAlmostEqual(oneq_coherence_err, 0.00049975, 6, "Error: 1Q Coherence Limit")
+            random.seed(123)
+            for num_qubits in [1, 2]:
+                for _ in range(100):
+                    t1s = [random.randint(100, 200) for _ in range(num_qubits)]
+                    t2s = [random.randint(100, 200) for _ in range(num_qubits)]
+                    time = random.randint(1, 10)
+                    self.assertAlmostEqual(
+                        rb.RBUtils.coherence_limit(num_qubits, t1s, t2s, time),
+                        rb.RBUtils.coherence_limit_error(num_qubits, time, t1s, t2s),
+                    )
 
+        self.assertAlmostEqual(oneq_coherence_err, 0.00049975, 6, "Error: 1Q Coherence Limit")
         self.assertAlmostEqual(twoq_coherence_err, 0.00597, 5, "Error: 2Q Coherence Limit")
 
+    def test_coherence_limit_error(self):
+        """Test coherence_limit_error."""
+        t1 = 100.0
+        t2 = 150.0
+        coherence_err_1q = rb.RBUtils.coherence_limit_error(1, 0.1, [t1], [t2])
+        coherence_err_2q = rb.RBUtils.coherence_limit_error(2, 0.5, [t1] * 2, [t2] * 2)
+        coherence_err_9q = rb.RBUtils.coherence_limit_error(9, 0.9, [t1] * 9, [t2] * 9)
+        self.assertAlmostEqual(coherence_err_1q, 0.00038873, 6, "Error: 1Q Coherence Limit")
+        self.assertAlmostEqual(coherence_err_2q, 0.00465046, 6, "Error: 2Q Coherence Limit")
+        self.assertAlmostEqual(coherence_err_9q, 0.04601531, 6, "Error: 9Q Coherence Limit")
+
+        self.assertAlmostEqual(
+            rb.RBUtils.coherence_limit_error(num_qubits=1, gate_length=0.1, t1s=[10], t2s=[50]),
+            rb.RBUtils.coherence_limit_error(num_qubits=1, gate_length=0.1, t1s=[10], t2s=[20]),
+            "T2 value must be truncated down to 2 * T1",
+        )
+        self.assertAlmostEqual(
+            rb.RBUtils.coherence_limit_error(num_qubits=1, gate_length=0.1, t1s=[10]),
+            rb.RBUtils.coherence_limit_error(num_qubits=1, gate_length=0.1, t1s=[10], t2s=[20]),
+            "If T2 values are not given, assume 2 * T1",
+        )
+
+        with self.assertRaises(ValueError):
+            rb.RBUtils.coherence_limit_error(num_qubits=2, gate_length=0.1, t1s=[10])
+
     def test_clifford_1_qubit_generation(self):
         """Verify 1-qubit clifford indeed generates the correct group"""
         clifford_dicts = [