Repository: PennyLaneAI/pennylane
Pull request #7617 (2025-06-05 19:30:28)

[Decomposition] decomposition of `SelectPauliRot` as quantum function
**Context:**

**Description of the Change:**

Translates the decomposition of `SelectPauliRot` into a quantum function.

**Benefits:**

**Possible Drawbacks:**

**Related GitHub Issues:**

[sc-92660]

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7617?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
Attention: Patch coverage is `97.22222%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 99.68%. Comparing base [(`a4a14be`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/a4a14be9b21c6cb2bbf949b62189f99ca137a2c9?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) to head [(`10f907c`)](https://app.codecov.io/gh/PennyLaneAI/pennylane/commit/10f907c88932a89eb72c2b5606e29b2c0ad8e701?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

| [Files with missing lines](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7617?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Patch % | Lines |
|---|---|---|
| [...ennylane/templates/subroutines/select\_pauli\_rot.py](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7617?src=pr&el=tree&filepath=pennylane%2Ftemplates%2Fsubroutines%2Fselect_pauli_rot.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3RlbXBsYXRlcy9zdWJyb3V0aW5lcy9zZWxlY3RfcGF1bGlfcm90LnB5) | 95.45% | [1 Missing :warning: ](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7617?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master    #7617      +/-   ##
==========================================
- Coverage   99.68%   99.68%   -0.01%     
==========================================
  Files         531      531              
  Lines       51896    51912      +16     
==========================================
+ Hits        51732    51747      +15     
- Misses        164      165       +1     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/PennyLaneAI/pennylane/pull/7617?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

<details><summary> :rocket: New features to boost your workflow: </summary>

- :snowflake: [Test Analytics](https://docs.codecov.com/docs/test-analytics): Detect flaky tests, report on failures, and find test suite problems.
</details>

Patch:
diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md
index 03c689629..46ef92ad3 100644
--- a/doc/releases/changelog-dev.md
+++ b/doc/releases/changelog-dev.md
@@ -93,6 +93,7 @@
 * A new template called :class:`~.SelectPauliRot` that applies a sequence of uniformly controlled rotations to a target qubit 
   is now available. This operator appears frequently in unitary decomposition and block encoding techniques. 
   [(#7206)](https://github.com/PennyLaneAI/pennylane/pull/7206)
+  [(#7617)](https://github.com/PennyLaneAI/pennylane/pull/7617)
 
   ```python
   angles = np.array([1.0, 2.0, 3.0, 4.0])
@@ -254,6 +255,9 @@
   * :class:`~.PCPhase`
     [(#7591)](https://github.com/PennyLaneAI/pennylane/pull/7591)
 
+  * :class:`~.IntegerComparator`
+    [(#7636)](https://github.com/PennyLaneAI/pennylane/pull/7636)
+
 * A new decomposition rule that uses a single work wire for decomposing multi-controlled operators is added.
   [(#7383)](https://github.com/PennyLaneAI/pennylane/pull/7383)
 
@@ -361,6 +365,13 @@
 
 <h3>Improvements 🛠</h3>
 
+* Caching with finite shots now always warns about the lack of expected noise.
+  [(#7644)](https://github.com/PennyLaneAI/pennylane/pull/7644)
+
+* `cache` now defaults to `"auto"` with `qml.execute`, matching the behavior of `QNode` and reducing the 
+  performance cost of using `qml.execute` for standard executions.
+  [(#7644)](https://github.com/PennyLaneAI/pennylane/pull/7644)
+
 * `qml.grad` and `qml.jacobian` can now handle inputs with dynamic shapes being captured into plxpr.
   [(#7544)](https://github.com/PennyLaneAI/pennylane/pull/7544/)
 
@@ -373,6 +384,9 @@
 * The decomposition of `qml.PCPhase` is now significantly more efficient for more than 2 qubits.
   [(#7166)](https://github.com/PennyLaneAI/pennylane/pull/7166)
 
+* The decomposition of :class:`~.IntegerComparator` is now significantly more efficient.
+  [(#7636)](https://github.com/PennyLaneAI/pennylane/pull/7636)
+
 * :class:`~.QubitUnitary` now supports a decomposition that is compatible with an arbitrary number of qubits. 
   This represents a fundamental improvement over the previous implementation, which was limited to two-qubit systems.
   [(#7277)](https://github.com/PennyLaneAI/pennylane/pull/7277)
@@ -500,8 +514,17 @@
 
 * Updated documentation check to remove duplicate docstring references. [(#7453)](https://github.com/PennyLaneAI/pennylane/pull/7453)
 
+* Improved performance for `qml.clifford_t_decomposition` transform by introducing caching support and changed the
+  default basis set of `qml.ops.sk_decomposition` to `(H, S, T)`, resulting in shorter decomposition sequences.
+  [(#7454)](https://github.com/PennyLaneAI/pennylane/pull/7454)
+
 <h3>Labs: a place for unified and rapid prototyping of research software 🧪</h3>
 
+* The imports of dependencies introduced by ``labs`` functionalities have been modified such that
+  these dependencies only have to be installed for the functions that use them, not to use
+  ``labs`` functionalities in general. This decouples the various submodules, and even functions
+  within the same submodule, from each other.
+  [(#7xxx)](https://github.com/PennyLaneAI/pennylane/pull/7xxx)
 
 * A new module :mod:`pennylane.labs.intermediate_reps <pennylane.labs.intermediate_reps>`
   provides functionality to compute intermediate representations for particular circuits.
@@ -537,6 +560,10 @@
 
 <h3>Breaking changes 💔</h3>
 
+* The default value of `cache` is now `"auto"` with `qml.execute`. Like `QNode`, `"auto"` only turns on caching
+  when `max_diff > 1`.
+  [(#7644)](https://github.com/PennyLaneAI/pennylane/pull/7644)
+
 * A new decomposition for two-qubit unitaries was implemented in `two_qubit_decomposition`.
   It ensures the correctness of the decomposition in some edge cases but uses 3 CNOT gates
   even if 2 CNOTs would suffice theoretically.
@@ -682,6 +709,11 @@ Here's a list of deprecations made this release. For a more detailed breakdown o
 
 * Fixed a failing integration test for `qml.QDrift`  which multiplied the operators of the decomposition incorrectly to evolve the state.
   [(#7621)](https://github.com/PennyLaneAI/pennylane/pull/7621)
+
+* The decomposition test in `assert_valid` no longer checks the matrix of the decomposition if the operator
+  does not define a matrix representation.
+  [(#7655)](https://github.com/PennyLaneAI/pennylane/pull/7655)
+
 <h3>Documentation 📝</h3>
 
 * Updated the circuit drawing for `qml.Select` to include two commonly used symbols for 
@@ -826,6 +858,7 @@ Here's a list of deprecations made this release. For a more detailed breakdown o
 This release contains contributions from (in alphabetical order):
 
 Guillermo Alonso-Linaje,
+Utkarsh Azad,
 Astral Cai,
 Yushao Chen,
 Marcus Edwards,
diff --git a/pennylane/_version.py b/pennylane/_version.py
index f03aebbd0..fd7da99da 100644
--- a/pennylane/_version.py
+++ b/pennylane/_version.py
@@ -16,4 +16,4 @@
 Version number (major.minor.patch[-label])
 """
 
-__version__ = "0.42.0-dev49"
+__version__ = "0.42.0-dev50"
diff --git a/pennylane/capture/__init__.py b/pennylane/capture/__init__.py
index ae61f7be0..beb5d9c39 100644
--- a/pennylane/capture/__init__.py
+++ b/pennylane/capture/__init__.py
@@ -161,6 +161,7 @@ If needed, developers can also override the implementation method of the primiti
 from typing import Callable, Type
 
 from .switches import disable, enable, enabled, pause
+from .promote_consts import promote_consts
 from .capture_meta import CaptureMeta, ABCCaptureMeta
 from .flatfn import FlatFn
 from .make_plxpr import make_plxpr, run_autograph
diff --git a/pennylane/capture/base_interpreter.py b/pennylane/capture/base_interpreter.py
index d0c2307d9..47cfe3946 100644
--- a/pennylane/capture/base_interpreter.py
+++ b/pennylane/capture/base_interpreter.py
@@ -35,6 +35,7 @@ from .primitives import (
     qnode_prim,
     while_loop_prim,
 )
+from .promote_consts import promote_consts
 
 FlattenedHigherOrderPrimitives: dict["jax.extend.core.Primitive", Callable] = {}
 """
@@ -458,34 +459,27 @@ def _(self, *dyn_shape, dimension, dtype, shape, sharding):
 
 
 @PlxprInterpreter.register_primitive(adjoint_transform_prim)
-def handle_adjoint_transform(self, *invals, jaxpr, lazy, n_consts):
+def handle_adjoint_transform(self, *invals, jaxpr, lazy):
     """Interpret an adjoint transform primitive."""
-    consts = invals[:n_consts]
-    args = invals[n_consts:]
-    jaxpr = jaxpr_to_jaxpr(copy(self), jaxpr, consts, *args)
-
-    return adjoint_transform_prim.bind(
-        *jaxpr.consts, *args, jaxpr=jaxpr.jaxpr, lazy=lazy, n_consts=len(jaxpr.consts)
-    )
+    jaxpr = jaxpr_to_jaxpr(copy(self), jaxpr, [], *invals)
+    jaxpr, new_invals = promote_consts(jaxpr, invals)
+    return adjoint_transform_prim.bind(*new_invals, jaxpr=jaxpr, lazy=lazy)
 
 
 # pylint: disable=too-many-arguments
 @PlxprInterpreter.register_primitive(ctrl_transform_prim)
-def handle_ctrl_transform(self, *invals, n_control, jaxpr, control_values, work_wires, n_consts):
+def handle_ctrl_transform(self, *invals, n_control, jaxpr, control_values, work_wires):
     """Interpret a ctrl transform primitive."""
-    consts = invals[:n_consts]
-    args = invals[n_consts:-n_control]
-    jaxpr = jaxpr_to_jaxpr(copy(self), jaxpr, consts, *args)
+    args = invals[:-n_control]
+    jaxpr = jaxpr_to_jaxpr(copy(self), jaxpr, [], *args)
+    jaxpr, new_invals = promote_consts(jaxpr, invals)
 
     return ctrl_transform_prim.bind(
-        *jaxpr.consts,
-        *args,
-        *invals[-n_control:],
+        *new_invals,
         n_control=n_control,
-        jaxpr=jaxpr.jaxpr,
+        jaxpr=jaxpr,
         control_values=control_values,
         work_wires=work_wires,
-        n_consts=len(jaxpr.consts),
     )
 
 
@@ -587,22 +581,18 @@ def handle_while_loop(
 
 # pylint: disable=too-many-arguments
 @PlxprInterpreter.register_primitive(qnode_prim)
-def handle_qnode(self, *invals, shots, qnode, device, execution_config, qfunc_jaxpr, n_consts):
+def handle_qnode(self, *invals, shots, qnode, device, execution_config, qfunc_jaxpr):
     """Handle a qnode primitive."""
-    consts = invals[:n_consts]
-    args = invals[n_consts:]
-
-    new_qfunc_jaxpr = jaxpr_to_jaxpr(copy(self), qfunc_jaxpr, consts, *args)
 
+    qfunc_jaxpr = jaxpr_to_jaxpr(copy(self), qfunc_jaxpr, [], *invals)
+    jaxpr, new_invals = promote_consts(qfunc_jaxpr, invals)
     return qnode_prim.bind(
-        *new_qfunc_jaxpr.consts,
-        *args,
+        *new_invals,
         shots=shots,
         qnode=qnode,
         device=device,
         execution_config=execution_config,
-        qfunc_jaxpr=new_qfunc_jaxpr.jaxpr,
-        n_consts=len(new_qfunc_jaxpr.consts),
+        qfunc_jaxpr=jaxpr,
     )
 
 
diff --git a/pennylane/capture/promote_consts.py b/pennylane/capture/promote_consts.py
new file mode 100644
index 000000000..5443d94d3
--- /dev/null
+++ b/pennylane/capture/promote_consts.py
@@ -0,0 +1,49 @@
+# Copyright 2018-2025 Xanadu Quantum Technologies Inc.
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+# pylint: disable=protected-access
+"""
+This submodule contains the `promote_consts` util.
+"""
+
+from pennylane.typing import TensorLike
+
+
+def promote_consts(
+    jaxpr: "jax.extend.core.ClosedJaxpr", args: list[TensorLike] | tuple[TensorLike, ...]
+) -> tuple["jax.extend.core.Jaxpr", tuple[TensorLike]]:
+    """Convert a closed jaxpr into a jaxpr without constants and the constants prepended
+    to the arguments.
+
+    Args:
+        jaxpr (jax.extend.core.ClosedJaxpr): a closed jaxpr with bound constants
+        args (tuple[TensorLike]): the arguments we will add the constants to the start of
+
+    Returns:
+
+        (jax.extend.core.Jaxpr, tuple[TensorLike]): A Jaxpr and new arguments with the constants pre-pended.
+
+    >>> import jax
+    >>> def f(x):
+    ...     return jax.numpy.array([2.0]) + x
+    >>> args = (jax.numpy.array([0.5]), )
+    >>> jaxpr = jax.make_jaxpr(f)(*args)
+    >>> jaxpr
+    { lambda a:f32[1]; b:f32[1]. let c:f32[1] = add a b in (c,) }
+    >>> qml.capture.promote_consts(jaxpr, args)
+    ({ lambda ; a:f32[1] b:f32[1]. let c:f32[1] = add a b in (c,) },
+    (array([2.0], dtype=float32), Array([0.5], dtype=float32)))
+
+    """
+    new_args = tuple(jaxpr.consts) + tuple(args)
+    jaxpr = jaxpr.jaxpr.replace(constvars=(), invars=jaxpr.jaxpr.constvars + jaxpr.jaxpr.invars)
+    return jaxpr, new_args
diff --git a/pennylane/decomposition/collect_resource_ops.py b/pennylane/decomposition/collect_resource_ops.py
index d3a26dffb..19668dd5a 100644
--- a/pennylane/decomposition/collect_resource_ops.py
+++ b/pennylane/decomposition/collect_resource_ops.py
@@ -35,26 +35,23 @@ class CollectResourceOps(FlattenedInterpreter):
 
 
 @CollectResourceOps.register_primitive(adjoint_transform_prim)
-def _(self, *invals, jaxpr, lazy, n_consts):  # pylint: disable=unused-argument
+def _(self, *invals, jaxpr, lazy):  # pylint: disable=unused-argument
     """Collect all operations in the base plxpr and create adjoint resource ops with them."""
-    consts = invals[:n_consts]
-    args = invals[n_consts:]
     child = CollectResourceOps()
-    child.eval(jaxpr, consts, *args)
+    child.eval(jaxpr, [], *invals)
     for op in child.state["ops"]:
         self.state["ops"].add(adjoint_resource_rep(op.op_type, op.params))
     return []
 
 
 @CollectResourceOps.register_primitive(ctrl_transform_prim)
-def _(self, *invals, n_control, jaxpr, n_consts, **params):
+def _(self, *invals, n_control, jaxpr, **params):
     """Collect all operations in the target plxpr and create controlled resource ops with them."""
 
-    consts = invals[:n_consts]
-    args = invals[n_consts:-n_control]
+    args = invals[:-n_control]
     control = invals[-n_control:]
     child = CollectResourceOps()
-    child.eval(jaxpr, consts, *args)
+    child.eval(jaxpr, [], *args)
 
     # Extract the resource parameters of this control transform
     control_values = params.get("control_values")
diff --git a/pennylane/devices/default_qubit.py b/pennylane/devices/default_qubit.py
index 060b31574..30e8de7dd 100644
--- a/pennylane/devices/default_qubit.py
+++ b/pennylane/devices/default_qubit.py
@@ -1004,6 +1004,7 @@ class DefaultQubit(Device):
             key=key,
             execution_config=execution_config,
         )
+
         return interpreter.eval(jaxpr, consts, *args)
 
     def _backprop_jvp(self, jaxpr, args, tangents, execution_config=None):
diff --git a/pennylane/devices/qubit/dq_interpreter.py b/pennylane/devices/qubit/dq_interpreter.py
index 5755df63e..dc6d750f8 100644
--- a/pennylane/devices/qubit/dq_interpreter.py
+++ b/pennylane/devices/qubit/dq_interpreter.py
@@ -231,12 +231,9 @@ def _(self, *invals, reset, postselect):
 
 
 @DefaultQubitInterpreter.register_primitive(adjoint_transform_prim)
-def _(self, *invals, jaxpr, n_consts, lazy=True):
-    # TODO: requires jaxpr -> list of ops first
-    consts = invals[:n_consts]
-    args = invals[n_consts:]
+def _(self, *invals, jaxpr, lazy=True):
     recorder = CollectOpsandMeas()
-    recorder.eval(jaxpr, consts, *args)
+    recorder.eval(jaxpr, [], *invals)
 
     ops = recorder.state["ops"]
     with pause():
@@ -252,13 +249,12 @@ def _(self, *invals, jaxpr, n_consts, lazy=True):
 
 # pylint: disable=too-many-arguments
 @DefaultQubitInterpreter.register_primitive(ctrl_transform_prim)
-def _(self, *invals, n_control, jaxpr, control_values, work_wires, n_consts):
+def _(self, *invals, n_control, jaxpr, control_values, work_wires):
     # TODO: requires jaxpr -> list of ops first
-    consts = invals[:n_consts]
     control_wires = invals[-n_control:]
-    args = invals[n_consts:-n_control]
+    args = invals[:-n_control]
     recorder = CollectOpsandMeas()
-    recorder.eval(jaxpr, consts, *args)
+    recorder.eval(jaxpr, [], *args)
 
     ops = recorder.state["ops"]
     with pause():
diff --git a/pennylane/labs/dla/variational_kak.py b/pennylane/labs/dla/variational_kak.py
index c3af78086..8c3230bb4 100644
--- a/pennylane/labs/dla/variational_kak.py
+++ b/pennylane/labs/dla/variational_kak.py
@@ -17,7 +17,6 @@ import warnings
 from datetime import datetime
 from functools import partial
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 import pennylane as qml
@@ -25,16 +24,23 @@ from pennylane.liealg import adjvec_to_op, op_to_adjvec
 from pennylane.operation import Operator
 from pennylane.pauli import PauliSentence
 
-has_jax = True
 try:
     import jax
     import jax.numpy as jnp
     import optax
 
     jax.config.update("jax_enable_x64", True)
+    has_jax = True
 except ImportError:
     has_jax = False
 
+try:
+    import matplotlib.pyplot as plt
+
+    has_plt = True
+except ImportError:
+    has_plt = False
+
 
 def variational_kak_adj(H, g, dims, adj, verbose=False, opt_kwargs=None, pick_min=False):
     r"""
@@ -90,7 +96,7 @@ def variational_kak_adj(H, g, dims, adj, verbose=False, opt_kwargs=None, pick_mi
             Cartan decomposition :math:`\mathfrak{g} = \mathfrak{k} \oplus (\tilde{\mathfrak{m}} \oplus \mathfrak{a})`
         adj (np.ndarray): Adjoint representation of dimension ``(dim_g, dim_g, dim_g)``,
             with the implicit ordering ``(k, mtilde, a)``.
-        verbose (bool): Plot the optimization
+        verbose (bool): Plot the optimization. Requires matplotlib to be installed (``pip install matplotlib``)
         opt_kwargs (dict): Keyword arguments for the optimization like initial starting values
             for :math:`\theta` of dimension ``(dim_k,)``, given as ``theta0``.
             Also includes ``n_epochs``, ``lr``, ``b1``, ``b2``, ``verbose``, ``interrupt_tol``, see :func:`~run_opt`
@@ -218,6 +224,11 @@ def variational_kak_adj(H, g, dims, adj, verbose=False, opt_kwargs=None, pick_mi
         raise ImportError(
             "jax and optax are required for variational_kak_adj. You can install them with pip install jax jaxlib optax."
         )  # pragma: no cover
+    if verbose >= 1 and not has_plt:  # pragma: no cover
+        print(
+            "variational_kak_adj requires matplotlib to display a figure with the optimization "
+            "progress (for verbose>=1). You can install it with pip install matplotlib"
+        )
 
     if opt_kwargs is None:
         opt_kwargs = {}
diff --git a/pennylane/labs/intermediate_reps/rowcol.py b/pennylane/labs/intermediate_reps/rowcol.py
index 81d82685f..6b8e29dd1 100644
--- a/pennylane/labs/intermediate_reps/rowcol.py
+++ b/pennylane/labs/intermediate_reps/rowcol.py
@@ -15,11 +15,20 @@
 
 from typing import Iterable
 
-import galois
 import networkx as nx
 import numpy as np
 from networkx.algorithms.approximation import steiner_tree
 
+try:
+    import galois
+
+    # Create the Galois number field F_2, in which we can create arrays in _get_S.
+    F_2 = galois.GF(2)
+    has_galois = True
+
+except ImportError:
+    has_galois = False
+
 
 def postorder_traverse(tree: nx.Graph, source: int, source_parent: int = None):
     """Post-order traverse a tree graph, starting from (but excluding) the node ``source``.
@@ -201,10 +210,6 @@ def _update(P: np.ndarray, cnots: list[tuple[int]], control: int, target: int):
     return P, cnots
 
 
-# Create the Galois number field F_2, in which we can create arrays in _get_S.
-F_2 = galois.GF(2)
-
-
 def _get_S(P: np.ndarray, idx: int, node_set: Iterable[int], mode: str):
     # Find S (S') either by simply extracting a column or by solving a linear system for the row
     if mode == "column":
@@ -298,6 +303,12 @@ def rowcol(P: np.ndarray, connectivity: nx.Graph = None, verbose: bool = False)
         list[tuple[int]]: Wire pairs for CNOTs that implement the parity matrix (control first,
         target second).
 
+    .. note::
+
+        This function requires the package ``galois`` to be installed. It can be installed via
+        ``pip install galois``, for details see
+        `its documentation <https://mhostetter.github.io/galois/latest/>`__.
+
     **Example**
 
     Here we compute the example 1 from `Wu et al. <https://arxiv.org/abs/1910.14478>`__ in code.
@@ -632,6 +643,11 @@ def rowcol(P: np.ndarray, connectivity: nx.Graph = None, verbose: bool = False)
         4: ─╰●────╰●────╰X─╰●────╰●─────────────╰●──────────╰X────┤
 
     """
+    if not has_galois:  # pragma: no cover
+        raise ImportError(
+            "rowcol requires the package galois. You can install it with pip install galois."
+        )  # pragma: no cover
+
     P = P.copy()
     connectivity = connectivity.copy()
     n = len(P)
diff --git a/pennylane/labs/vibrational/christiansen_utils.py b/pennylane/labs/vibrational/christiansen_utils.py
index 653229b2b..724c58d98 100644
--- a/pennylane/labs/vibrational/christiansen_utils.py
+++ b/pennylane/labs/vibrational/christiansen_utils.py
@@ -16,13 +16,19 @@ import itertools
 from pathlib import Path
 from tempfile import TemporaryDirectory
 
-import h5py
 import numpy as np
 from scipy.special import factorial
 
 from pennylane import concurrency
 
-# pylint: disable = redefined-outer-name,
+# pylint: disable = redefined-outer-name
+
+try:
+    import h5py
+
+    has_h5py = True
+except ImportError:
+    has_h5py = False
 
 
 def _cform_onemode_kinetic(freqs, n_states, num_workers=1, backend="serial", path=None):
@@ -1045,10 +1051,20 @@ def christiansen_integrals(pes, n_states=16, cubic=False, num_workers=1, backend
         backend (string): the executor backend from the list of supported backends.
             Available options : "mp_pool", "cf_procpool", "cf_threadpool", "serial", "mpi4py_pool", "mpi4py_comm". Default value is set to "serial".
 
-
     Returns:
         TensorLike[float]: the integrals for the Christiansen Hamiltonian
+
+    .. note::
+
+        This function requires the ``h5py`` package to be installed.
+        It can be installed with ``pip install h5py``.
+
     """
+    if not has_h5py:  # pragma: no cover
+        raise ImportError(
+            "christiansen_integrals requires the h5py package. "
+            "You can install it with pip install h5py"
+        )  # pragma: no cover
     with TemporaryDirectory() as path:
         ham_cform_onebody = _cform_onemode(
             pes, n_states, num_workers=num_workers, backend=backend, path=path
@@ -1084,10 +1100,20 @@ def christiansen_integrals_dipole(pes, n_states=16, num_workers=1, backend="seri
         backend (string): the executor backend from the list of supported backends.
             Available options : "mp_pool", "cf_procpool", "cf_threadpool", "serial", "mpi4py_pool", "mpi4py_comm". Default value is set to "serial".
 
-
     Returns:
         TensorLike[float]: the integrals for the Christiansen dipole operator
+
+    .. note::
+
+        This function requires the ``h5py`` package to be installed.
+        It can be installed with ``pip install h5py``.
+
     """
+    if not has_h5py:  # pragma: no cover
+        raise ImportError(
+            "christiansen_integrals_dipole requires the h5py package. "
+            "You can install it with pip install h5py"
+        )  # pragma: no cover
     with TemporaryDirectory() as path:
         dipole_cform_onebody = _cform_onemode_dipole(
             pes, n_states, num_workers=num_workers, backend=backend, path=path
diff --git a/pennylane/labs/zxopt/basic_optimization.py b/pennylane/labs/zxopt/basic_optimization.py
index 660ddfddc..bcfc0c2ca 100644
--- a/pennylane/labs/zxopt/basic_optimization.py
+++ b/pennylane/labs/zxopt/basic_optimization.py
@@ -13,8 +13,6 @@
 # limitations under the License.
 """Optimization pass ``basic_optimization`` from pyzx using ZX calculus."""
 
-import pyzx as zx
-
 import pennylane as qml
 from pennylane.tape import QuantumScript, QuantumScriptBatch
 from pennylane.transforms import transform
@@ -22,6 +20,13 @@ from pennylane.typing import PostprocessingFn
 
 from .util import _tape2pyzx
 
+try:
+    import pyzx as zx
+
+    has_zx = True
+except ImportError:
+    has_zx = False
+
 
 def null_postprocessing(results):
     """A postprocesing function returned by a transform that only converts the batch of results
@@ -93,6 +98,12 @@ def basic_optimization(tape: QuantumScript) -> tuple[QuantumScriptBatch, Postpro
         3: ─╰X──H──T─╰●──H───────────┤
 
     """
+    if not has_zx:  # pragma: no cover
+        raise ImportError(
+            "basic_optimization requires the package pyzx. "
+            "You can install it with pip install pyzx"
+        )  # pragma: no cover
+
     pyzx_circ = _tape2pyzx(tape)
 
     pyzx_circ = pyzx_circ.to_basic_gates()
diff --git a/pennylane/labs/zxopt/full_optimize.py b/pennylane/labs/zxopt/full_optimize.py
index 46820e1cf..210eb28a4 100644
--- a/pennylane/labs/zxopt/full_optimize.py
+++ b/pennylane/labs/zxopt/full_optimize.py
@@ -14,8 +14,6 @@
 """Optimization pass ``full_optimize`` from pyzx using ZX calculus."""
 import warnings
 
-import pyzx as zx
-
 import pennylane as qml
 from pennylane.tape import QuantumScript, QuantumScriptBatch
 from pennylane.transforms import transform
@@ -23,6 +21,13 @@ from pennylane.typing import PostprocessingFn
 
 from .util import _tape2pyzx
 
+try:
+    import pyzx as zx
+
+    has_zx = True
+except ImportError:
+    has_zx = False
+
 
 def null_postprocessing(results):
     """A postprocesing function returned by a transform that only converts the batch of results
@@ -118,6 +123,11 @@ def full_optimize(
             new_circ = full_optimize(circ, clifford_t_args = {"epsilon": 0.1})
 
     """
+    if not has_zx:  # pragma: no cover
+        raise ImportError(
+            "full_optimize requires the package pyzx. " "You can install it with pip install pyzx"
+        )  # pragma: no cover
+
     try:
         pyzx_circ = _tape2pyzx(tape)
 
diff --git a/pennylane/labs/zxopt/full_reduce.py b/pennylane/labs/zxopt/full_reduce.py
index 182a7ea22..7391b1ca8 100644
--- a/pennylane/labs/zxopt/full_reduce.py
+++ b/pennylane/labs/zxopt/full_reduce.py
@@ -13,9 +13,6 @@
 # limitations under the License.
 """Optimization pass ``full_reduce`` from pyzx using ZX calculus."""
 
-import pyzx as zx
-from pyzx.graph.base import BaseGraph
-
 import pennylane as qml
 from pennylane.tape import QuantumScript, QuantumScriptBatch
 from pennylane.transforms import transform
@@ -23,6 +20,14 @@ from pennylane.typing import PostprocessingFn
 
 from .util import _tape2pyzx
 
+try:
+    import pyzx as zx
+    from pyzx.graph.base import BaseGraph
+
+    has_zx = True
+except ImportError:
+    has_zx = False
+
 
 def null_postprocessing(results):
     """A postprocesing function returned by a transform that only converts the batch of results
@@ -103,6 +108,10 @@ def full_reduce(tape: QuantumScript) -> tuple[QuantumScriptBatch, Postprocessing
 
     The original circuit has five :class:`~T` gates which are reduced to just one.
     """
+    if not has_zx:  # pragma: no cover
+        raise ImportError(
+            "full_reduce requires the package pyzx. " "You can install it with pip install pyzx"
+        )  # pragma: no cover
 
     pyzx_circ = _tape2pyzx(tape)
 
diff --git a/pennylane/labs/zxopt/todd.py b/pennylane/labs/zxopt/todd.py
index bf629d9b7..dd0cf1e4a 100644
--- a/pennylane/labs/zxopt/todd.py
+++ b/pennylane/labs/zxopt/todd.py
@@ -13,8 +13,6 @@
 # limitations under the License.
 """Third order duplicate and destroy (TODD) optimization method from pyzx, using ZX calculus."""
 
-import pyzx as zx
-
 import pennylane as qml
 from pennylane.tape import QuantumScript, QuantumScriptBatch
 from pennylane.transforms import transform
@@ -22,6 +20,13 @@ from pennylane.typing import PostprocessingFn
 
 from .util import _tape2pyzx
 
+try:
+    import pyzx as zx
+
+    has_zx = True
+except ImportError:
+    has_zx = False
+
 
 def null_postprocessing(results):
     """A postprocesing function returned by a transform that only converts the batch of results
@@ -100,6 +105,10 @@ def todd(
     The original five T gates are reduced to just one.
 
     """
+    if not has_zx:  # pragma: no cover
+        raise ImportError(
+            "todd requires the pyzx package. You can install it with pip install pyzx"
+        )  # pragma: no cover
 
     pyzx_circ = _tape2pyzx(tape)
 
diff --git a/pennylane/labs/zxopt/util.py b/pennylane/labs/zxopt/util.py
index ae9e35d57..831f7f185 100644
--- a/pennylane/labs/zxopt/util.py
+++ b/pennylane/labs/zxopt/util.py
@@ -15,7 +15,12 @@
 
 import re
 
-from pyzx import Circuit
+try:
+    from pyzx import Circuit
+
+    has_zx = True
+except ImportError:
+    has_zx = False
 
 
 def _remove_measurement_patterns(text):
@@ -48,6 +53,11 @@ def _tape2pyzx(tape):
         pyzx.Graph: graph instance of the pyzx circuit
 
     """
+    if not has_zx:  # pragma: no cover
+        raise ImportError(
+            "_tape2pyzx requires the package pyzx. " "You can install it with pip install pyzx"
+        )  # pragma: no cover
+
     g = Circuit.from_qasm(_remove_measurement_patterns(tape.to_openqasm()))
 
     return g
diff --git a/pennylane/ops/functions/assert_valid.py b/pennylane/ops/functions/assert_valid.py
index af11f9b35..455697cc4 100644
--- a/pennylane/ops/functions/assert_valid.py
+++ b/pennylane/ops/functions/assert_valid.py
@@ -167,12 +167,12 @@ def _test_decomposition_rule(op, rule: DecompositionRule, heuristic_resources=Fa
         tape.operations.insert(0, qml.Projector([0] * len(work_wires), wires=work_wires))
 
     # Tests that the decomposition produces the same matrix
-    op_matrix = qml.matrix(op, wire_order=all_wires)
-
-    decomp_matrix = qml.matrix(tape, wire_order=all_wires)
-    assert qml.math.allclose(
-        op_matrix, decomp_matrix
-    ), "decomposition must produce the same matrix as the operator."
+    if op.has_matrix:
+        op_matrix = op.matrix(wire_order=all_wires)
+        decomp_matrix = qml.matrix(tape, wire_order=all_wires)
+        assert qml.math.allclose(
+            op_matrix, decomp_matrix
+        ), "decomposition must produce the same matrix as the operator."
 
 
 def _check_matrix(op):
diff --git a/pennylane/ops/op_math/adjoint.py b/pennylane/ops/op_math/adjoint.py
index f0eead95b..1af97d1c7 100644
--- a/pennylane/ops/op_math/adjoint.py
+++ b/pennylane/ops/op_math/adjoint.py
@@ -198,13 +198,11 @@ def _get_adjoint_qfunc_prim():
     adjoint_prim.prim_type = "higher_order"
 
     @adjoint_prim.def_impl
-    def _(*args, jaxpr, lazy, n_consts):
+    def _(*args, jaxpr, lazy):
         from pennylane.tape.plxpr_conversion import CollectOpsandMeas
 
-        consts = args[:n_consts]
-        args = args[n_consts:]
         collector = CollectOpsandMeas()
-        collector.eval(jaxpr, consts, *args)
+        collector.eval(jaxpr, [], *args)
         for op in reversed(collector.state["ops"]):
             adjoint(op, lazy=lazy)
         return []
@@ -228,13 +226,13 @@ def _capture_adjoint_transform(qfunc: Callable, lazy=True) -> Callable:
         abstracted_axes, abstract_shapes = qml.capture.determine_abstracted_axes(args)
         jaxpr = jax.make_jaxpr(partial(qfunc, **kwargs), abstracted_axes=abstracted_axes)(*args)
         flat_args = jax.tree_util.tree_leaves(args)
+        jaxpr, new_args = qml.capture.promote_consts(
+            jaxpr, tuple(abstract_shapes) + tuple(flat_args)
+        )
         adjoint_prim.bind(
-            *jaxpr.consts,
-            *abstract_shapes,
-            *flat_args,
-            jaxpr=jaxpr.jaxpr,
+            *new_args,
+            jaxpr=jaxpr,
             lazy=lazy,
-            n_consts=len(jaxpr.consts),
         )
 
     return new_qfunc
diff --git a/pennylane/ops/op_math/controlled.py b/pennylane/ops/op_math/controlled.py
index 537ead1b2..468005178 100644
--- a/pennylane/ops/op_math/controlled.py
+++ b/pennylane/ops/op_math/controlled.py
@@ -288,16 +288,15 @@ def _get_ctrl_qfunc_prim():
     ctrl_prim.prim_type = "higher_order"
 
     @ctrl_prim.def_impl
-    def _(*args, n_control, jaxpr, control_values, work_wires, n_consts):
+    def _(*args, n_control, jaxpr, control_values, work_wires):
         from pennylane.tape.plxpr_conversion import CollectOpsandMeas
 
-        consts = args[:n_consts]
         control_wires = args[-n_control:]
-        args = args[n_consts:-n_control]
+        args = args[:-n_control]
 
         collector = CollectOpsandMeas()
         with qml.QueuingManager.stop_recording():
-            collector.eval(jaxpr, consts, *args)
+            collector.eval(jaxpr, [], *args)
 
         for op in collector.state["ops"]:
             ctrl(op, control_wires, control_values, work_wires)
@@ -324,17 +323,18 @@ def _capture_ctrl_transform(qfunc: Callable, control, control_values, work_wires
             *args
         )
         flat_args = jax.tree_util.tree_leaves(args)
+        jaxpr, new_args = qml.capture.promote_consts(
+            jaxpr, tuple(abstract_shapes) + tuple(flat_args)
+        )
+
         control_wires = qml.wires.Wires(control)  # make sure is iterable
         ctrl_prim.bind(
-            *jaxpr.consts,
-            *abstract_shapes,
-            *flat_args,
+            *new_args,
             *control_wires,
-            jaxpr=jaxpr.jaxpr,
+            jaxpr=jaxpr,
             n_control=len(control_wires),
             control_values=control_values,
             work_wires=work_wires,
-            n_consts=len(jaxpr.consts),
         )
 
     return new_qfunc
diff --git a/pennylane/ops/op_math/decompositions/solovay_kitaev.py b/pennylane/ops/op_math/decompositions/solovay_kitaev.py
index b537e3f64..8f2572e61 100644
--- a/pennylane/ops/op_math/decompositions/solovay_kitaev.py
+++ b/pennylane/ops/op_math/decompositions/solovay_kitaev.py
@@ -275,7 +275,7 @@ def _group_commutator_decompose(matrix, tol=1e-5):
     return w_hat, v_hat
 
 
-def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), basis_length=10):
+def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("H", "S", "T"), basis_length=10):
     r"""Approximate an arbitrary single-qubit gate in the Clifford+T basis using the `Solovay-Kitaev algorithm <https://arxiv.org/abs/quant-ph/0505030>`_.
 
     This method implements the Solovay-Kitaev decomposition algorithm that approximates any single-qubit
@@ -294,7 +294,7 @@ def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), ba
             a greater number of passes. Default is ``5``.
         basis_set (tuple[str]): Basis set to be used for the decomposition and building an approximate set internally.
             It accepts the following gate terms: ``('X', 'Y', 'Z', 'H', 'T', 'T*', 'S', 'S*')``, where ``*`` refers
-            to the gate adjoint. Default value is ``('T', 'T*', 'H')``.
+            to the gate adjoint. Default value is ``('H', 'S', 'T')``.
         basis_length (int): Maximum expansion length of Clifford+T sequences in the internally-built approximate set.
             Default is ``10``.
 
@@ -316,7 +316,7 @@ def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), ba
 
         op  = qml.RZ(np.pi/3, wires=0)
 
-        # Get the gate decomposition in ['T', 'T*', 'H']
+        # Get the gate decomposition in ['H', 'S', 'T']
         ops = qml.ops.sk_decomposition(op, epsilon=1e-3)
 
         # Get the approximate matrix from the ops
@@ -408,11 +408,18 @@ def sk_decomposition(op, epsilon, *, max_depth=5, basis_set=("T", "T*", "H"), ba
         )
 
     # Map the wires to that of the operation and queue
-    [map_tape], _ = qml.map_wires(new_tape, wire_map={0: op.wires[0]}, queue=True)
+    if queuing := QueuingManager.recording():
+        QueuingManager.remove(op)
+
+    if (op_wire := op.wires[0]) != 0:
+        [new_tape], _ = qml.map_wires(new_tape, wire_map={0: op_wire}, queue=True)
+    else:
+        if queuing:
+            _ = [qml.apply(op) for op in new_tape.operations]
 
     # Get phase information based on the decomposition effort
     phase = approx_set_gph[index] - gate_gph
     global_phase = qml.GlobalPhase(qml.math.array(phase, like=interface))
 
     # Return the gates from the mapped tape and global phase
-    return map_tape.operations + [global_phase]
+    return new_tape.operations + [global_phase]
diff --git a/pennylane/ops/qubit/arithmetic_ops.py b/pennylane/ops/qubit/arithmetic_ops.py
index e96e1e2ce..50e515ee4 100644
--- a/pennylane/ops/qubit/arithmetic_ops.py
+++ b/pennylane/ops/qubit/arithmetic_ops.py
@@ -15,6 +15,8 @@
 This submodule contains the discrete-variable quantum operations that perform
 arithmetic operations on their input states.
 """
+from collections import Counter
+
 # pylint: disable=arguments-differ
 from copy import copy
 from typing import Optional
@@ -22,10 +24,14 @@ from typing import Optional
 import numpy as np
 
 import pennylane as qml
-from pennylane.decomposition import add_decomps, register_resources
+from pennylane.decomposition import (
+    add_decomps,
+    register_condition,
+    register_resources,
+    resource_rep,
+)
 from pennylane.decomposition.symbolic_decomposition import pow_involutory, self_adjoint
 from pennylane.operation import FlatPytree, Operation
-from pennylane.ops import Identity
 from pennylane.typing import TensorLike
 from pennylane.wires import Wires, WiresLike
 
@@ -409,6 +415,8 @@ class IntegerComparator(Operation):
 
     grad_method = None
 
+    resource_keys = {"num_wires", "value", "geq", "num_work_wires"}
+
     def _flatten(self) -> FlatPytree:
         hp = self.hyperparameters
         metadata = (
@@ -456,6 +464,15 @@ class IntegerComparator(Operation):
 
         super().__init__(wires=total_wires)
 
+    @property
+    def resource_params(self) -> dict:
+        return {
+            "num_wires": len(self.wires),
+            "value": self.value,
+            "geq": self.geq,
+            "num_work_wires": len(self.hyperparameters["work_wires"]),
+        }
+
     def label(
         self,
         decimals: Optional[int] = None,
@@ -560,11 +577,12 @@ class IntegerComparator(Operation):
 
         **Example:**
 
-        >>> print(qml.IntegerComparator.compute_decomposition(4, wires=[0, 1, 2, 3]))
-        [MultiControlledX(wires=[0, 1, 2, 3], control_values=[1, 0, 0]),
-         MultiControlledX(wires=[0, 1, 2, 3], control_values=[1, 0, 1]),
-         MultiControlledX(wires=[0, 1, 2, 3], control_values=[1, 1, 0]),
-         MultiControlledX(wires=[0, 1, 2, 3], control_values=[1, 1, 1])]
+        >>> print(qml.draw(qml.IntegerComparator.compute_decomposition)(4, wires=[0, 1, 2, 3]))
+        0: ─╭●────╭●────╭●────┤
+        1: ─├●──X─├●────├●──X─┤
+        2: ─│─────├●──X─├●──X─┤
+        3: ─╰X────╰X────╰X────┤
+
         """
 
         if not isinstance(value, int):
@@ -573,35 +591,23 @@ class IntegerComparator(Operation):
         if wires is None:
             raise ValueError("Must specify the wires that the operation acts on.")
 
-        if len(wires) > 1:
-            control_wires = Wires(wires[:-1])
-            wires = Wires(wires[-1])
-        else:
+        if len(wires) < 2:
             raise ValueError(
                 f"IntegerComparator: wrong number of wires. {len(wires)} wire(s) given. Need at least 2."
             )
 
-        small_val = not geq and value == 0
-        large_val = geq and value > 2 ** len(control_wires) - 1
-        if small_val or large_val:
-            gates = [Identity(wires[0])]
+        work_wires = Wires([]) if work_wires is None else Wires(work_wires)
+        with qml.queuing.AnnotatedQueue() as q:
+            if geq:
+                _integer_comparator_ge_decomposition(wires, value, work_wires)
+            else:
+                _integer_comparator_lt_decomposition(wires, value, work_wires)
 
-        else:
-            values = range(value, 2 ** (len(control_wires))) if geq else range(value)
-            binary = "0" + str(len(control_wires)) + "b"
-            control_values_list = [format(n, binary) for n in values]
-            gates = []
-            for control_values in control_values_list:
-                control_values = [int(n) for n in control_values]
-                gates.append(
-                    qml.MultiControlledX(
-                        wires=control_wires + wires,
-                        control_values=control_values,
-                        work_wires=work_wires,
-                    )
-                )
+        if qml.queuing.QueuingManager.recording():
+            for op in q.queue:
+                qml.apply(op)
 
-        return gates
+        return q.queue
 
     @property
     def control_wires(self) -> Wires:
@@ -612,3 +618,284 @@ class IntegerComparator(Operation):
 
     def pow(self, z: int) -> list["IntegerComparator"]:
         return super().pow(z % 2)
+
+
+def _integer_comparator_lt_resource(num_wires, value, num_work_wires, **_):
+
+    if value == 0:
+        return {}
+
+    if value > 2 ** (num_wires - 1) - 1:
+        return {qml.X: 1}
+
+    num_controls = num_wires - 1
+    binary_str = format(value, f"0{num_controls}b")
+    last_significant = binary_str.rfind("1")
+    gate_counts = {resource_rep(qml.X): (last_significant + 1) * 2}
+
+    first_significant = binary_str.find("1")
+    gate_counts[
+        resource_rep(
+            qml.MultiControlledX,
+            num_control_wires=first_significant + 1,
+            num_work_wires=num_work_wires + num_wires - 2 - first_significant,
+            num_zero_control_values=0,
+            work_wire_type="dirty",
+        )
+    ] = 1
+
+    while (first_significant := binary_str.find("1", first_significant + 1)) != -1:
+        gate_counts[
+            resource_rep(
+                qml.MultiControlledX,
+                num_control_wires=first_significant + 1,
+                num_work_wires=num_work_wires + num_wires - 2 - first_significant,
+                num_zero_control_values=0,
+                work_wire_type="dirty",
+            )
+        ] = 1
+
+    return gate_counts
+
+
+@register_condition(lambda geq, **_: not geq)
+@register_resources(_integer_comparator_lt_resource)
+def _integer_comparator_lt_decomposition(wires, value, work_wires, **_):
+    """Decompose the ``IntegerComparator`` for when the flipping condition is ``n < value``.
+
+    This decomposition uses the minimum number of ``MultiControlledX`` gates. For a given value,
+    we first convert it to binary, and iteratively look for the significant bits. For example,
+    with 6 control wires, if the value is 22, which is 010110 in 6-bit binary, we observe
+    that all 6-bit numbers that start with 00 will satisfy the flipping condition, so we apply
+    a ``MultiControlledX`` with only the first two wires as controls, and 00 as the control values.
+    Then we look for the next significant bit, and observe that 22 starts with 0101. Therefore,
+    all 6-bit numbers that start with 0100 will also satisfy the flipping condition, so we apply
+    a ``MultiControlledX`` with the first four wires as controls, and 0100 as the control values.
+    This continues until we add a ``MultiControlledX`` for every significant bit in the value.
+
+    .. code-block:: pycon
+
+        0: ─╭○─╭○─╭○─┤
+        1: ─├○─├●─├●─┤
+        2: ─│──├○─├○─┤
+        3: ─│──├○─├●─┤
+        4: ─│──│──├○─┤
+        6: ─╰X─╰X─╰X─┤
+
+    If we decompose this circuit one level further, we get
+
+    .. code-block:: pycon
+
+        0: ──X─╭●──X──X─╭●──X──X─╭●──X─┤
+        1: ──X─├●──X────├●───────├●────┤
+        2: ────│───X────├●──X──X─├●──X─┤
+        3: ────│───X────├●──X────├●────┤
+        4: ────│────────│───X────├●──X─┤
+        6: ────╰X───────╰X───────╰X────┤
+
+    And we observe that the ``PauliX`` gates used to flip the control values can be merged:
+
+    .. code-block:: pycon
+
+        0: ──X─╭●────╭●────╭●──X─┤
+        1: ──X─├●──X─├●────├●────┤
+        2: ──X─│─────├●────├●──X─┤
+        3: ──X─│─────├●──X─├●────┤
+        4: ──X─│─────│─────├●──X─┤
+        6: ────╰X────╰X────╰X────┤
+
+    """
+
+    # If the value is zero, the flipping condition is never satisfied.
+    if value == 0:
+        return
+
+    num_controls = len(wires) - 1
+
+    # If the value is larger than the maximum value that can be represented by
+    # the number of control bits, the flipping condition is always satisfied, in
+    # which case we apply an X gate to the target wire and terminate.
+    if value > 2**num_controls - 1:
+        qml.X(wires[-1])
+        return
+
+    # Track which control bits have been flipped back
+    control_value_tracker = [0] * num_controls
+
+    # First apply X to all wires until the last significant bit to flip control values to 1.
+    binary_str = format(value, f"0{num_controls}b")
+    last_significant = binary_str.rfind("1")
+    for i in range(last_significant + 1):
+        qml.X(wires[i])
+
+    # The flipping condition is satisfied if all bits from the first bit to the
+    # first non-zero bit of the value are zeroes.
+    first_significant = binary_str.find("1")
+    qml.MultiControlledX(
+        wires=wires[: first_significant + 1] + wires[-1:],
+        work_wires=wires[first_significant + 1 : -1] + work_wires,
+    )
+    control_value_tracker[first_significant] = 1
+    qml.X(wires[first_significant])
+
+    # If the wire corresponding to the first significant bit of the value is 1, then we
+    # iteratively look for the next significant bit, and apply a flip conditioned on all
+    # bits from the last significant bit to this next significant bit being zeroes.
+    while (first_significant := binary_str.find("1", first_significant + 1)) != -1:
+        qml.MultiControlledX(
+            wires=wires[: first_significant + 1] + wires[-1:],
+            work_wires=wires[first_significant + 1 : -1] + work_wires,
+        )
+        control_value_tracker[first_significant] = 1
+        qml.X(wires[first_significant])
+
+    for i in range(last_significant + 1):
+        if control_value_tracker[i] == 0:
+            qml.X(wires[i])
+
+
+def _integer_comparator_ge_resource(num_wires, value, num_work_wires, **_):
+
+    # If the value is 0, the flipping condition is always satisfied.
+    if value == 0:
+        return {qml.X: 1}
+
+    num_controls = num_wires - 1
+
+    if value > 2**num_controls - 1:
+        return {}
+
+    binary_str = format(value, f"0{num_controls}b")
+    first_zero = binary_str.find("0")
+
+    if first_zero == -1:
+        return {
+            resource_rep(
+                qml.MultiControlledX,
+                num_control_wires=num_controls,
+                num_work_wires=num_work_wires,
+                num_zero_control_values=0,
+                work_wire_type="dirty",
+            ): 1
+        }
+
+    gate_set = Counter()
+
+    gate_set[
+        resource_rep(
+            qml.MultiControlledX,
+            num_control_wires=first_zero + 1,
+            num_work_wires=num_work_wires + num_wires - 2 - first_zero,
+            num_zero_control_values=0,
+            work_wire_type="dirty",
+        )
+    ] = 1
+    gate_set[resource_rep(qml.X)] = 2
+
+    while (first_zero := binary_str.find("0", first_zero + 1)) != -1:
+        gate_set[
+            resource_rep(
+                qml.MultiControlledX,
+                num_control_wires=first_zero + 1,
+                num_work_wires=num_work_wires + num_wires - 2 - first_zero,
+                num_zero_control_values=0,
+                work_wire_type="dirty",
+            )
+        ] = 1
+        gate_set[resource_rep(qml.X)] += 2
+
+    gate_set[
+        resource_rep(
+            qml.MultiControlledX,
+            num_control_wires=num_controls,
+            num_work_wires=num_work_wires,
+            num_zero_control_values=0,
+            work_wire_type="dirty",
+        )
+    ] += 1
+
+    return dict(gate_set)
+
+
+@register_condition(lambda geq, **_: geq)
+@register_resources(_integer_comparator_ge_resource)
+def _integer_comparator_ge_decomposition(wires, value, work_wires, **_):
+    """Decompose the ``IntegerComparator`` for when the flipping condition is ``n >= value``.
+
+    This decomposition rule mirrors the implementation for the ``n < value`` case.
+
+    """
+
+    # If the value is 0, the flipping condition is always satisfied.
+    if value == 0:
+        qml.X(wires[-1])
+        return
+
+    num_controls = len(wires) - 1
+
+    # If the value is larger than the maximum value that can be represented by
+    # the number of control bits, the flipping condition is never satisfied,
+    if value > 2**num_controls - 1:
+        return
+
+    # Track which control bits have been flipped
+    control_value_tracker = [0] * num_controls
+
+    binary_str = format(value, f"0{num_controls}b")
+    first_zero = binary_str.find("0")
+
+    if first_zero == -1:
+        # If the value happens to be the all-one state, then we apply a single MCX
+        qml.MultiControlledX(wires=wires, work_wires=work_wires)
+        return
+
+    qml.MultiControlledX(
+        wires=wires[: first_zero + 1] + wires[-1:],
+        work_wires=wires[first_zero + 1 : -1] + work_wires,
+    )
+    control_value_tracker[first_zero] = 1
+    qml.X(wires[first_zero])
+
+    while (first_zero := binary_str.find("0", first_zero + 1)) != -1:
+        qml.MultiControlledX(
+            wires=wires[: first_zero + 1] + wires[-1:],
+            work_wires=wires[first_zero + 1 : -1] + work_wires,
+        )
+        control_value_tracker[first_zero] = 1
+        qml.X(wires[first_zero])
+
+    # The last MCX corresponds to the equal case.
+    qml.MultiControlledX(wires=wires, work_wires=work_wires)
+
+    for i in range(num_controls):
+        if control_value_tracker[i]:
+            qml.X(wires[i])
+
+
+def _integer_comparator_flip_geq_resource(num_wires, value, num_work_wires, geq, **_):
+    """Resource estimation for flipping the geq condition."""
+    return {
+        qml.X: 1,
+        resource_rep(
+            qml.IntegerComparator,
+            num_wires=num_wires,
+            value=value,
+            geq=not geq,
+            num_work_wires=num_work_wires,
+        ): 1,
+    }
+
+
+@register_resources(_integer_comparator_flip_geq_resource)
+def _integer_comparator_flip_geq(value, geq, wires, work_wires, **_):
+    """Decompose the IntegerComparator by flipping geq to lt or vice versa."""
+    qml.X(wires[-1])
+    IntegerComparator(value, wires, geq=not geq, work_wires=work_wires)
+
+
+add_decomps(
+    IntegerComparator,
+    _integer_comparator_lt_decomposition,
+    _integer_comparator_ge_decomposition,
+    _integer_comparator_flip_geq,
+)
diff --git a/pennylane/tape/plxpr_conversion.py b/pennylane/tape/plxpr_conversion.py
index 590eeaa4c..f9cf7eb3d 100644
--- a/pennylane/tape/plxpr_conversion.py
+++ b/pennylane/tape/plxpr_conversion.py
@@ -97,13 +97,11 @@ class CollectOpsandMeas(FlattenedInterpreter):
 
 
 @CollectOpsandMeas.register_primitive(adjoint_transform_prim)
-def _(self, *invals, jaxpr, lazy, n_consts):
+def _(self, *invals, jaxpr, lazy):
     """Handle an adjoint transform primitive by collecting the operations in the jaxpr, and
     then applying their adjoint in reverse order."""
-    consts = invals[:n_consts]
-    args = invals[n_consts:]
     child = CollectOpsandMeas()
-    child.eval(jaxpr, consts, *args)
+    child.eval(jaxpr, [], *invals)
     assert child.state
 
     for op in reversed(child.state["ops"]):
@@ -113,16 +111,15 @@ def _(self, *invals, jaxpr, lazy, n_consts):
 
 
 @CollectOpsandMeas.register_primitive(ctrl_transform_prim)
-def _(self, *invals, n_control, jaxpr, n_consts, **params):
+def _(self, *invals, n_control, jaxpr, **params):
     """Handle a control transform primitive by collecting the operations in the jaxpr,
     and then applying their controlled versions.
     """
-    consts = invals[:n_consts]
-    args = invals[n_consts:-n_control]
+    args = invals[:-n_control]
     control = invals[-n_control:]
 
     child = CollectOpsandMeas()
-    child.eval(jaxpr, consts, *args)
+    child.eval(jaxpr, [], *args)
     assert child.state
 
     for op in child.state["ops"]:
@@ -190,13 +187,10 @@ def _(self, *invals, jaxpr, n_consts, **params):
 
 @CollectOpsandMeas.register_primitive(qnode_prim)
 def _(
-    self, *invals, shots, qnode, device, execution_config, qfunc_jaxpr, n_consts
+    self, *invals, shots, qnode, device, execution_config, qfunc_jaxpr
 ):  # pylint: disable=too-many-arguments
-    consts = invals[:n_consts]
-    args = invals[n_consts:]
-
     child = CollectOpsandMeas()
-    out = child.eval(qfunc_jaxpr, consts, *args)
+    out = child.eval(qfunc_jaxpr, [], *invals)
     assert child.state
     self.state["ops"].extend(child.state["ops"])
     self.state["measurements"].extend(child.state["measurements"])
diff --git a/pennylane/templates/state_preparations/mottonen.py b/pennylane/templates/state_preparations/mottonen.py
index 9571632cb..68b8fe9af 100644
--- a/pennylane/templates/state_preparations/mottonen.py
+++ b/pennylane/templates/state_preparations/mottonen.py
@@ -141,12 +141,9 @@ def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
         control_wires (array[int]): wires that act as control
         target_wire (int): wire that acts as target
 
-    Returns:
-          list[.Operator]: sequence of operators defined by this function
     """
-    op_list = []
-    gray_code_rank = len(control_wires)
 
+    gray_code_rank = len(control_wires)
     theta = compute_theta(alpha, num_qubits=gray_code_rank)
 
     if gray_code_rank == 0:
@@ -155,14 +152,14 @@ def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
             or qml.math.requires_grad(theta)
             or qml.math.all(theta[..., 0] != 0.0)
         ):
-            op_list.append(gate(theta[..., 0], wires=[target_wire]))
-        return op_list
+            gate(theta[..., 0], wires=[target_wire])
+        return
 
     code = gray_code(gray_code_rank)
     control_indices = np.log2(code ^ np.roll(code, -1)).astype(int)
 
-    # For abstract or differentiated theta we will never skip a rotation. Likewise if there is at
-    # least one non-zero angle (per batch if batched) for all rotations.
+    # For abstract or differentiated theta we will never skip a rotation. Likewise, if there
+    # is at least one non-zero angle (per batch if batched) for all rotations.
     skip_none = qml.math.is_abstract(theta) or qml.math.requires_grad(theta)
     if not skip_none:
         nonzero = (
@@ -172,9 +169,32 @@ def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
     for i, control_index in enumerate(control_indices):
         # If we do not _never_ skip, we might skip _some_ rotation
         if skip_none or qml.math.all(theta[..., i] != 0.0):
-            op_list.append(gate(theta[..., i], wires=[target_wire]))
-        op_list.append(qml.CNOT(wires=[control_wires[control_index], target_wire]))
-    return op_list
+            gate(theta[..., i], wires=[target_wire])
+        qml.CNOT(wires=[control_wires[control_index], target_wire])
+
+
+def _uniform_rotation_dagger_ops(gate, alpha, control_wires, target_wire):
+    r"""Returns a list of operators that applies a uniformly-controlled rotation to the target qubit.
+
+    Args:
+        gate (.Operation): gate to be applied, needs to have exactly one parameter
+        alpha (tensor_like): angles to decompose the uniformly-controlled rotation into multi-controlled rotations
+        control_wires (array[int]): wires that act as control
+        target_wire (int): wire that acts as target
+
+    Returns:
+          list[.Operator]: sequence of operators defined by this function
+
+    """
+
+    with qml.queuing.AnnotatedQueue() as q:
+        _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire)
+
+    if qml.queuing.QueuingManager.recording():
+        for op in q.queue:
+            qml.apply(op)
+
+    return q.queue
 
 
 def _get_alpha_z(omega, n, k):
@@ -396,7 +416,7 @@ class MottonenStatePreparation(Operation):
             alpha_y_k = _get_alpha_y(a, len(wires_reverse), k)
             control = wires_reverse[k:]
             target = wires_reverse[k - 1]
-            op_list.extend(_apply_uniform_rotation_dagger(qml.RY, alpha_y_k, control, target))
+            op_list.extend(_uniform_rotation_dagger_ops(qml.RY, alpha_y_k, control, target))
 
         # If necessary, apply inverse z rotation cascade to prepare correct phases of amplitudes
         if (
@@ -409,9 +429,7 @@ class MottonenStatePreparation(Operation):
                 control = wires_reverse[k:]
                 target = wires_reverse[k - 1]
                 if len(alpha_z_k) > 0:
-                    op_list.extend(
-                        _apply_uniform_rotation_dagger(qml.RZ, alpha_z_k, control, target)
-                    )
+                    op_list.extend(_uniform_rotation_dagger_ops(qml.RZ, alpha_z_k, control, target))
 
             global_phase = qml.math.sum(-1 * qml.math.angle(state_vector) / len(state_vector))
             op_list.extend([qml.GlobalPhase(global_phase, wires=wires)])
diff --git a/pennylane/templates/subroutines/select_pauli_rot.py b/pennylane/templates/subroutines/select_pauli_rot.py
index 0c8a7f785..8b46e45ba 100644
--- a/pennylane/templates/subroutines/select_pauli_rot.py
+++ b/pennylane/templates/subroutines/select_pauli_rot.py
@@ -16,6 +16,7 @@ Contains the SelectPauliRot template.
 """
 
 import pennylane as qml
+from pennylane.decomposition import add_decomps, adjoint_resource_rep, register_resources
 from pennylane.operation import Operation
 from pennylane.templates.state_preparations.mottonen import _apply_uniform_rotation_dagger
 
@@ -167,20 +168,43 @@ class SelectPauliRot(Operation):
         control_wires = qml.wires.Wires(control_wires)
         target_wire = qml.wires.Wires(target_wire)
 
-        op_list = []
+        with qml.queuing.AnnotatedQueue() as q:
+            decompose_select_pauli_rot(angles, control_wires + target_wire, rot_axis)
 
-        if rot_axis == "X":
-            op_list.append(qml.Hadamard(target_wire))
-        elif rot_axis == "Y":
-            op_list.extend([qml.adjoint(qml.S(target_wire)), qml.Hadamard(target_wire)])
+        if qml.queuing.QueuingManager.recording():
+            for op in q.queue:
+                qml.apply(op)
 
-        op_list.extend(
-            _apply_uniform_rotation_dagger(qml.RZ, angles, control_wires[::-1], target_wire[0])
-        )
+        return q.queue
+
+
+def _select_pauli_rot_resource(num_wires, rot_axis):
+    return {
+        qml.RZ: 2 ** (num_wires - 1),
+        qml.CNOT: 2 ** (num_wires - 1) if num_wires > 1 else 0,
+        qml.H: 0 if rot_axis == "Z" else 2,
+        qml.S: 1 if rot_axis == "Y" else 0,
+        adjoint_resource_rep(qml.S, {}): 1 if rot_axis == "Y" else 0,
+    }
+
+
+@register_resources(_select_pauli_rot_resource)
+def decompose_select_pauli_rot(angles, wires, rot_axis, **__):
+    r"""Decomposes the SelectPauliRot"""
+
+    if rot_axis == "X":
+        qml.Hadamard(wires[-1])
+    elif rot_axis == "Y":
+        qml.adjoint(qml.S(wires[-1]))
+        qml.Hadamard(wires[-1])
+
+    _apply_uniform_rotation_dagger(qml.RZ, angles, wires[-2::-1], wires[-1])
+
+    if rot_axis == "X":
+        qml.Hadamard(wires[-1])
+    elif rot_axis == "Y":
+        qml.Hadamard(wires[-1])
+        qml.S(wires[-1])
 
-        if rot_axis == "X":
-            op_list.append(qml.Hadamard(target_wire))
-        elif rot_axis == "Y":
-            op_list.extend([qml.Hadamard(target_wire), qml.S(target_wire)])
 
-        return op_list
+add_decomps(SelectPauliRot, decompose_select_pauli_rot)
diff --git a/pennylane/transforms/decompose.py b/pennylane/transforms/decompose.py
index 8e0f9f986..7757276d2 100644
--- a/pennylane/transforms/decompose.py
+++ b/pennylane/transforms/decompose.py
@@ -331,9 +331,8 @@ def _get_plxpr_decompose():  # pylint: disable=missing-docstring, too-many-state
 
     # pylint: disable=too-many-arguments
     @DecomposeInterpreter.register_primitive(ctrl_transform_prim)
-    def _(self, *invals, n_control, jaxpr, control_values, work_wires, n_consts):
-        consts = invals[:n_consts]
-        args = invals[n_consts:-n_control]
+    def _(self, *invals, n_control, jaxpr, control_values, work_wires):
+        args = invals[:-n_control]
         control_wires = invals[-n_control:]
 
         unroller = ControlTransformInterpreter(
@@ -341,7 +340,7 @@ def _get_plxpr_decompose():  # pylint: disable=missing-docstring, too-many-state
         )
 
         def wrapper(*inner_args):
-            return unroller.eval(jaxpr, consts, *inner_args)
+            return unroller.eval(jaxpr, [], *inner_args)
 
         jaxpr = jax.make_jaxpr(wrapper)(*args)
         return self.eval(jaxpr.jaxpr, jaxpr.consts, *args)
diff --git a/pennylane/transforms/decompositions/clifford_t_transform.py b/pennylane/transforms/decompositions/clifford_t_transform.py
index b2d7d4adc..5b6ba4cfc 100644
--- a/pennylane/transforms/decompositions/clifford_t_transform.py
+++ b/pennylane/transforms/decompositions/clifford_t_transform.py
@@ -15,6 +15,7 @@
 
 import math
 import warnings
+from functools import lru_cache, partial
 from itertools import product
 
 import pennylane as qml
@@ -64,6 +65,10 @@ _CLIFFORD_T_GATES = tuple(_CLIFFORD_T_ONE_GATES + _CLIFFORD_T_TWO_GATES) + (qml.
 # Gates to be skipped during decomposition
 _SKIP_OP_TYPES = (qml.Barrier, qml.Snapshot, qml.WireCut)
 
+# Stores the cache of a specified size for the decomposition function
+# that is used to decompose the RZ gates in the Clifford+T basis.
+_CLIFFORD_T_CACHE = None
+
 
 def _check_clifford_op(op, use_decomposition=False):
     r"""Checks if an operator is Clifford or not.
@@ -330,12 +335,71 @@ def _merge_param_gates(operations, merge_ops=None):
     return merged_ops, number_ops
 
 
-# pylint: disable=too-many-branches
+# NOTE: The cache size is set to 2000, which is a reasonable size for the mapping
+# decomposition in circuit with up to 500 wires. This is based on the fact that the
+# decompositions being mapped would contain {H, S, T} (with or without {Sadj, Tadj}).
+@lru_cache(maxsize=2000)
+def _map_wires(op, wire):
+    """Maps the operator to the provided wire."""
+    return op.map_wires({0: wire})
+
+
+class _CachedCallable:
+    """A class to cache the decomposition of the operators.
+
+    Args:
+        method (str): The method to be used for decomposition.
+        epsilon (float): The maximum permissible operator norm error for the decomposition.
+        cache_size (int): The size of the cache built for the decomposition function based on the angle.
+        **method_kwargs: Keyword argument to pass options for the ``method`` used for decompositions.
+    """
+
+    def __init__(self, method, epsilon, cache_size, **method_kwargs):
+        match method:
+            case "sk":
+                self.decompose_fn = lru_cache(maxsize=cache_size)(
+                    partial(sk_decomposition, epsilon=epsilon, **method_kwargs)
+                )
+            case _:
+                raise NotImplementedError(
+                    f"Currently we only support Solovay-Kitaev ('sk') decomposition, got {method}"
+                )
+
+        self.method = method
+        self.epsilon = epsilon
+        self.cache_size = cache_size
+        self.method_kwargs = method_kwargs
+        self.query = lru_cache(maxsize=cache_size)(self.cached_decompose)
+
+    def compatible(self, method, epsilon, cache_size, **method_kwargs):
+        """Check equality based on `method`, `epsilon` and `method_kwargs`."""
+        return (
+            self.method == method
+            and self.epsilon <= epsilon
+            and self.cache_size <= cache_size
+            and self.method_kwargs == method_kwargs
+        )
+
+    def cached_decompose(self, op):
+        """Decomposes the angle into a sequence of gates."""
+        f_adj = op.data[0] < 2 * math.pi  # 4 * pi / 2
+        cached_op = op if f_adj else op.adjoint()
+
+        seq = self.decompose_fn(cached_op)
+        if f_adj:
+            return seq
+
+        adj = [qml.adjoint(s, lazy=False) for s in reversed(seq)]
+        return adj[1:] + adj[:1]
+
+
+# pylint: disable= too-many-nested-blocks, too-many-branches, too-many-statements, unnecessary-lambda-assignment
 @transform
 def clifford_t_decomposition(
     tape: QuantumScript,
     epsilon=1e-4,
     method="sk",
+    cache_size=1000,
     **method_kwargs,
 ) -> tuple[QuantumScriptBatch, PostprocessingFn]:
     r"""Decomposes a circuit into the Clifford+T basis.
@@ -355,6 +419,7 @@ def clifford_t_decomposition(
         tape (QNode or QuantumTape or Callable): The quantum circuit to be decomposed.
         epsilon (float): The maximum permissible operator norm error of the complete circuit decomposition. Defaults to ``0.0001``.
         method (str): Method to be used for Clifford+T decomposition. Default value is ``"sk"`` for Solovay-Kitaev.
+        cache_size (int): The size of the cache built for the decomposition function based on the angle. Defaults to ``1000``.
         **method_kwargs: Keyword argument to pass options for the ``method`` used for decompositions.
 
     Returns:
@@ -450,26 +515,32 @@ def clifford_t_decomposition(
         # Compute the per-gate epsilon value
         epsilon /= number_ops or 1
 
-        # Every decomposition implementation should have the following shape:
+        # _CACHED_DECOMPOSE is a global variable that caches the decomposition function,
+        # where the implementation of each function should have the following signature:
         # def decompose_fn(op: Operator, epsilon: float, **method_kwargs) -> List[Operator]
         # note: the last operator in the decomposition must be a GlobalPhase
 
-        # Build the approximation set for Solovay-Kitaev decomposition
-        if method == "sk":
-            decompose_fn = sk_decomposition
-
-        else:
-            raise NotImplementedError(
-                f"Currently we only support Solovay-Kitaev ('sk') decomposition, got {method}"
-            )
+        # Build the decomposition cache based on the method
+        global _CLIFFORD_T_CACHE  # pylint: disable=global-statement
+        if _CLIFFORD_T_CACHE is None or not _CLIFFORD_T_CACHE.compatible(
+            method, epsilon, cache_size, **method_kwargs
+        ):
+            _CLIFFORD_T_CACHE = _CachedCallable(method, epsilon, cache_size, **method_kwargs)
 
         decomp_ops = []
         phase = new_operations.pop().data[0]
         for op in new_operations:
             if isinstance(op, qml.RZ):
-                clifford_ops = decompose_fn(op, epsilon, **method_kwargs)
-                phase += qml.math.convert_like(clifford_ops.pop().data[0], phase)
-                decomp_ops.extend(clifford_ops)
+                # If simplifies to Identity, skip it
+                if not (op_param := op.simplify().data):
+                    continue
+                # Decompose the RZ operation with a default wire
+                clifford_ops = _CLIFFORD_T_CACHE.query(qml.RZ(op_param[0], [0]))
+                # Extract the global phase from the last operation
+                phase += qml.math.convert_like(clifford_ops[-1].data[0], phase)
+                # Map the operations to the original wires
+                op_wire = op.wires[0]
+                decomp_ops.extend([_map_wires(cl_op, op_wire) for cl_op in clifford_ops[:-1]])
             else:
                 decomp_ops.append(op)
 
@@ -478,7 +549,9 @@ def clifford_t_decomposition(
             decomp_ops.append(qml.GlobalPhase(phase))
 
     # Construct a new tape with the expanded set of operations
+    # and then clear `decomp_ops` list to free up the memory
     new_tape = compiled_tape.copy(operations=decomp_ops)
+    decomp_ops.clear()
 
     # Perform a final attempt of simplification before return
     [new_tape], _ = cancel_inverses(new_tape)
diff --git a/pennylane/transforms/defer_measurements.py b/pennylane/transforms/defer_measurements.py
index 3f1a4eb2e..086e29c54 100644
--- a/pennylane/transforms/defer_measurements.py
+++ b/pennylane/transforms/defer_measurements.py
@@ -419,6 +419,7 @@ def _get_plxpr_defer_measurements():
                 for branch, value in condition.items():
                     # When reduce_postselected is True, some branches can be ()
                     cur_consts = invals[consts_slices[i]]
+                    jaxpr = jaxpr.replace(constvars=(), invars=jaxpr.constvars + jaxpr.invars)
                     qml.cond(value, ctrl_transform_prim.bind)(
                         *cur_consts,
                         *args,
@@ -427,7 +428,6 @@ def _get_plxpr_defer_measurements():
                         n_control=len(control_wires),
                         control_values=branch,
                         work_wires=None,
-                        n_consts=len(cur_consts),
                     )
 
         return [None] * len(jaxpr_branches[0].outvars)
diff --git a/pennylane/workflow/_cache_transform.py b/pennylane/workflow/_cache_transform.py
index 03dc85640..70262b5cf 100644
--- a/pennylane/workflow/_cache_transform.py
+++ b/pennylane/workflow/_cache_transform.py
@@ -43,7 +43,7 @@ def _cache_transform(tape: QuantumScript, cache: MutableMapping):
     def cache_hit_postprocessing(_results: ResultBatch) -> Result:
         result = cache[tape.hash]
         if result is not None:
-            if tape.shots and getattr(cache, "_persistent_cache", True):
+            if tape.shots:
                 warnings.warn(_CACHED_EXECUTION_WITH_FINITE_SHOTS_WARNINGS, UserWarning)
             return result
 
diff --git a/pennylane/workflow/_capture_qnode.py b/pennylane/workflow/_capture_qnode.py
index 6cab0575e..85e55ec74 100644
--- a/pennylane/workflow/_capture_qnode.py
+++ b/pennylane/workflow/_capture_qnode.py
@@ -116,7 +116,7 @@ import jax
 from jax.interpreters import ad, batching, mlir
 
 import pennylane as qml
-from pennylane.capture import FlatFn, QmlPrimitive
+from pennylane.capture import FlatFn, QmlPrimitive, promote_consts
 from pennylane.exceptions import CaptureError
 from pennylane.logging import debug_logger
 from pennylane.typing import TensorLike
@@ -192,15 +192,12 @@ qnode_prim.prim_type = "higher_order"
 # pylint: disable=too-many-arguments
 @debug_logger
 @qnode_prim.def_impl
-def _(*args, qnode, shots, device, execution_config, qfunc_jaxpr, n_consts, batch_dims=None):
+def _(*args, qnode, shots, device, execution_config, qfunc_jaxpr, batch_dims=None):
     if shots != device.shots:
         raise NotImplementedError(
             "Overriding shots is not yet supported with the program capture execution."
         )
 
-    consts = args[:n_consts]
-    non_const_args = args[n_consts:]
-
     device_program = device.preprocess_transforms(execution_config)
     if batch_dims is not None:
         temp_all_args = []
@@ -211,56 +208,48 @@ def _(*args, qnode, shots, device, execution_config, qfunc_jaxpr, n_consts, batc
                 temp_all_args.append(a[tuple(slices)])
             else:
                 temp_all_args.append(a)
-        temp_consts = temp_all_args[:n_consts]
-        temp_args = temp_all_args[n_consts:]
     else:
-        temp_consts = consts
-        temp_args = non_const_args
+        temp_all_args = args
 
     # Expand user transforms applied to the qfunc
+    accumulated_consts = ()
     if getattr(qfunc_jaxpr.eqns[0].primitive, "prim_type", "") == "transform":
         transformed_func = qml.capture.expand_plxpr_transforms(
-            partial(qml.capture.eval_jaxpr, qfunc_jaxpr, temp_consts)
+            partial(qml.capture.eval_jaxpr, qfunc_jaxpr, [])
         )
 
-        qfunc_jaxpr = jax.make_jaxpr(transformed_func)(*temp_args)
-        temp_consts = qfunc_jaxpr.consts
-        qfunc_jaxpr = qfunc_jaxpr.jaxpr
+        qfunc_jaxpr = jax.make_jaxpr(transformed_func)(*temp_all_args)
+        qfunc_jaxpr, accumulated_consts = promote_consts(qfunc_jaxpr, accumulated_consts)
 
     # Expand user transforms applied to the qnode
-    qfunc_jaxpr = qnode.transform_program(qfunc_jaxpr, temp_consts, *temp_args)
-    temp_consts = qfunc_jaxpr.consts
-    qfunc_jaxpr = qfunc_jaxpr.jaxpr
+    qfunc_jaxpr = qnode.transform_program(qfunc_jaxpr, [], *accumulated_consts, *temp_all_args)
+    qfunc_jaxpr, accumulated_consts = promote_consts(qfunc_jaxpr, accumulated_consts)
 
     # Apply device preprocessing transforms
     graph_enabled = qml.decomposition.enabled_graph()
     try:
         qml.decomposition.disable_graph()
-        qfunc_jaxpr = device_program(qfunc_jaxpr, temp_consts, *temp_args)
+        qfunc_jaxpr = device_program(qfunc_jaxpr, [], *temp_all_args)
+        qfunc_jaxpr, accumulated_consts = promote_consts(qfunc_jaxpr, accumulated_consts)
     finally:
         if graph_enabled:
             qml.decomposition.enable_graph()
-    consts = qfunc_jaxpr.consts
-    qfunc_jaxpr = qfunc_jaxpr.jaxpr
 
-    partial_eval = partial(
-        device.eval_jaxpr, qfunc_jaxpr, consts, execution_config=execution_config
-    )
+    partial_eval = partial(device.eval_jaxpr, qfunc_jaxpr, [], execution_config=execution_config)
     if batch_dims is None:
-        return partial_eval(*non_const_args)
-    return jax.vmap(partial_eval, batch_dims[n_consts:])(*non_const_args)
+        return partial_eval(*accumulated_consts, *args)
+
+    return jax.vmap(partial_eval, batch_dims)(*accumulated_consts, *args)
 
 
 # pylint: disable=unused-argument
 @debug_logger
 @qnode_prim.def_abstract_eval
-def _(*args, qnode, shots, device, execution_config, qfunc_jaxpr, n_consts, batch_dims=None):
+def _(*args, qnode, shots, device, execution_config, qfunc_jaxpr, batch_dims=None):
 
     mps = qfunc_jaxpr.outvars
 
-    batch_shape = (
-        _get_batch_shape(args[n_consts:], batch_dims[n_consts:]) if batch_dims is not None else ()
-    )
+    batch_shape = _get_batch_shape(args, batch_dims) if batch_dims is not None else ()
 
     return _get_shapes_for(
         *mps, shots=shots, num_device_wires=len(device.wires), batch_shape=batch_shape
@@ -277,7 +266,6 @@ def _qnode_batching_rule(
     device,
     execution_config,
     qfunc_jaxpr,
-    n_consts,
 ):
     """
     Batching rule for the ``qnode`` primitive.
@@ -286,25 +274,10 @@ def _qnode_batching_rule(
     """
 
     for idx, (arg, batch_dim) in enumerate(zip(batched_args, batch_dims)):
-
-        if _is_scalar_tensor(arg):
-            continue
-
-        # Regardless of their shape, jax.vmap automatically inserts `None` as the batch dimension for constants.
-        # However, if the constant is not a standard JAX type, the batch dimension is not inserted at all.
-        # How to handle this case is still an open question. For now, we raise a warning and give the user full flexibility.
-        if idx < n_consts:
-            warn(
-                f"Constant argument at index {idx} is not scalar. "
-                "This may lead to unintended behavior or wrong results if the argument is provided "
-                "using parameter broadcasting to a quantum operation that supports batching.",
-                UserWarning,
-            )
-
         # To resolve this ambiguity, we might add more properties to the AbstractOperator
         # class to indicate which operators support batching and check them here.
         # As above, at this stage we raise a warning and give the user full flexibility.
-        elif arg.size > 1 and batch_dim is None:
+        if not _is_scalar_tensor(arg) and arg.size > 1 and batch_dim is None:
             warn(
                 f"Argument at index {idx} has size > 1 but its batch dimension is None. "
                 "This may lead to unintended behavior or wrong results if the argument is provided "
@@ -319,7 +292,6 @@ def _qnode_batching_rule(
         device=device,
         execution_config=execution_config,
         qfunc_jaxpr=qfunc_jaxpr,
-        n_consts=n_consts,
         batch_dims=batch_dims,
     )
 
@@ -476,8 +448,10 @@ def _extract_qfunc_jaxpr(qnode, abstracted_axes, *args, **kwargs):
             "flow functions like for_loop, while_loop, etc."
         ) from exc
 
+    qfunc_jaxpr, consts = promote_consts(qfunc_jaxpr, ())
+
     assert flat_fn.out_tree is not None, "out_tree should be set by call to flat_fn"
-    return qfunc_jaxpr, flat_fn.out_tree
+    return qfunc_jaxpr, flat_fn.out_tree, consts
 
 
 def capture_qnode(qnode: "qml.QNode", *args, **kwargs) -> "qml.typing.Result":
@@ -569,7 +543,7 @@ def capture_qnode(qnode: "qml.QNode", *args, **kwargs) -> "qml.typing.Result":
     cache_key = _get_jaxpr_cache_key(flat_dynamic_args, flat_static_args, kwargs, abstracted_axes)
 
     if cached_value := qnode.capture_cache.get(cache_key, None):
-        qfunc_jaxpr, config, out_tree = cached_value
+        qfunc_jaxpr, config, out_tree, consts = cached_value
     else:
         config = construct_execution_config(
             qnode, resolve=False
@@ -581,20 +555,21 @@ def capture_qnode(qnode: "qml.QNode", *args, **kwargs) -> "qml.typing.Result":
             # as the original dynamic arguments
             abstracted_axes = jax.tree_util.tree_unflatten(dynamic_args_struct, abstracted_axes)
 
-        qfunc_jaxpr, out_tree = _extract_qfunc_jaxpr(qnode, abstracted_axes, *args, **kwargs)
+        qfunc_jaxpr, out_tree, consts = _extract_qfunc_jaxpr(
+            qnode, abstracted_axes, *args, **kwargs
+        )
 
-        qnode.capture_cache[cache_key] = (qfunc_jaxpr, config, out_tree)
+        qnode.capture_cache[cache_key] = (qfunc_jaxpr, config, out_tree, consts)
 
     res = qnode_prim.bind(
-        *qfunc_jaxpr.consts,
         *abstract_shapes,
+        *consts,
         *flat_dynamic_args,
         shots=shots,
         qnode=qnode,
         device=qnode.device,
         execution_config=config,
-        qfunc_jaxpr=qfunc_jaxpr.jaxpr,
-        n_consts=len(qfunc_jaxpr.consts),
+        qfunc_jaxpr=qfunc_jaxpr,
     )
 
     return jax.tree_util.tree_unflatten(out_tree, res)
diff --git a/pennylane/workflow/_setup_transform_program.py b/pennylane/workflow/_setup_transform_program.py
index fd1ba0c3c..325e5e0be 100644
--- a/pennylane/workflow/_setup_transform_program.py
+++ b/pennylane/workflow/_setup_transform_program.py
@@ -90,9 +90,13 @@ def _setup_transform_program(
     _prune_dynamic_transform(outer_transform_program, inner_transform_program)
 
     # If caching is desired but an explicit cache is not provided, use an ``LRUCache``.
+    if cache == "auto":
+        if resolved_execution_config.derivative_order == 1:
+            cache = None
+        else:
+            cache = True
     if cache is True:
         cache = LRUCache(maxsize=cachesize)
-        setattr(cache, "_persistent_cache", False)
 
     # Ensure that ``cache`` is not a Boolean to simplify downstream code.
     elif cache is False:
diff --git a/pennylane/workflow/execution.py b/pennylane/workflow/execution.py
index d42450d56..97aa4dbc0 100644
--- a/pennylane/workflow/execution.py
+++ b/pennylane/workflow/execution.py
@@ -47,7 +47,7 @@ def execute(
     *,
     transform_program: TransformProgram = None,
     grad_on_execution: Literal[True, False, "best"] = "best",
-    cache: Union[None, bool, dict, Cache] = True,
+    cache: Union[None, bool, dict, Cache, Literal["auto"]] = "auto",
     cachesize: int = 10000,
     max_diff: int = 1,
     device_vjp: Union[bool, None] = False,
@@ -75,8 +75,13 @@ def execute(
             if the device is queried for the gradient; gradient transform
             functions available in ``qml.gradients`` are only supported on the backward
             pass. The 'best' option chooses automatically between the two options and is default.
-        cache (None, bool, dict, Cache): Whether to cache evaluations. This can result in
-            a significant reduction in quantum evaluations during gradient computations.
+        cache="auto" (str or bool or dict or Cache): Whether to cache evalulations.
+            ``"auto"`` indicates to cache only when ``max_diff > 1``. This can result in
+            a reduction in quantum evaluations during higher order gradient computations.
+            If ``True``, a cache with corresponding ``cachesize`` is created for each batch
+            execution. If ``False``, no caching is used. You may also pass your own cache
+            to be used; this can be any object that implements the special methods
+            ``__getitem__()``, ``__setitem__()``, and ``__delitem__()``, such as a dictionary.
         cachesize (int): the size of the cache.
         max_diff (int): If ``diff_method`` is a gradient transform, this option specifies
             the maximum number of derivatives to support. Increasing this value allows
diff --git a/pennylane/workflow/qnode.py b/pennylane/workflow/qnode.py
index 7563bb9dc..47bbb6979 100644
--- a/pennylane/workflow/qnode.py
+++ b/pennylane/workflow/qnode.py
@@ -613,7 +613,6 @@ class QNode:
             _validate_jax_version()
         self.diff_method = diff_method
         _validate_diff_method(self.device, self.diff_method)
-        cache = (max_diff > 1) if cache == "auto" else cache
 
         self.capture_cache = LRUCache(maxsize=1000)
         if isinstance(static_argnums, int):


Test patch:
diff --git a/tests/capture/test_base_interpreter.py b/tests/capture/test_base_interpreter.py
index 3a39e58b0..d429aad90 100644
--- a/tests/capture/test_base_interpreter.py
+++ b/tests/capture/test_base_interpreter.py
@@ -402,7 +402,7 @@ class TestHigherOrderPrimitiveRegistrations:
 
         jaxpr2 = jax.make_jaxpr(ConstAdder()(f))(0.5)
         assert jaxpr2.consts == [scalar]
-        assert len(jaxpr2.eqns[0].params["jaxpr"].constvars) == 1
+        assert len(jaxpr2.eqns[0].params["jaxpr"].constvars) == 0
 
     def test_ctrl_transform(self):
         """Test the higher order ctrl transform."""
@@ -439,12 +439,14 @@ class TestHigherOrderPrimitiveRegistrations:
             qml.ctrl(g, control)(x)
 
         jaxpr = jax.make_jaxpr(f)(0.5, 1)
+
         assert len(jaxpr.consts) == 0
         assert len(jaxpr.eqns[0].params["jaxpr"].constvars) == 0
+        assert len(jaxpr.eqns[0].params["jaxpr"].invars) == 1
 
         jaxpr2 = jax.make_jaxpr(ConstAdder()(f))(0.5, 1)
         assert jaxpr2.consts == [scalar]
-        assert len(jaxpr2.eqns[0].params["jaxpr"].constvars) == 1
+        assert len(jaxpr2.eqns[0].params["jaxpr"].constvars) == 0
 
     def test_cond(self):
         """Test the cond higher order primitive."""
@@ -697,7 +699,7 @@ class TestHigherOrderPrimitiveRegistrations:
 
         jaxpr2 = jax.make_jaxpr(ConstAdder()(f))()
         assert jaxpr2.consts == [scalar]
-        assert len(jaxpr2.eqns[0].params["qfunc_jaxpr"].constvars) == 1
+        assert len(jaxpr2.eqns[0].params["qfunc_jaxpr"].constvars) == 0
 
     @pytest.mark.parametrize("grad_f", (qml.grad, qml.jacobian))
     def test_grad_and_jac(self, grad_f):
diff --git a/tests/capture/test_nested_plxpr.py b/tests/capture/test_nested_plxpr.py
index de2e975c9..392b1885a 100644
--- a/tests/capture/test_nested_plxpr.py
+++ b/tests/capture/test_nested_plxpr.py
@@ -146,7 +146,6 @@ class TestAdjointQfunc:
         jaxpr = jax.make_jaxpr(workflow)(0.5)
 
         assert jaxpr.eqns[0].primitive == adjoint_transform_prim
-        assert jaxpr.eqns[0].params["n_consts"] == 1
         assert len(jaxpr.eqns[0].invars) == 2  # one const, one arg
 
         with qml.queuing.AnnotatedQueue() as q:
@@ -170,9 +169,8 @@ class TestAdjointQfunc:
         assert len(plxpr.eqns) == 1
         grad_eqn = plxpr.eqns[0]
         assert grad_eqn.primitive == grad_prim
-        assert set(grad_eqn.params.keys()) == {"argnum", "n_consts", "jaxpr", "method", "h"}
+        assert set(grad_eqn.params.keys()) == {"argnum", "jaxpr", "n_consts", "method", "h"}
         assert grad_eqn.params["argnum"] == [0]
-        assert grad_eqn.params["n_consts"] == 0
         assert grad_eqn.params["method"] is None
         assert grad_eqn.params["h"] is None
         assert len(grad_eqn.params["jaxpr"].eqns) == 1
@@ -296,7 +294,6 @@ class TestCtrlQfunc:
         assert plxpr.eqns[0].params["control_values"] == [True]
         assert plxpr.eqns[0].params["n_control"] == 1
         assert plxpr.eqns[0].params["work_wires"] is None
-        assert plxpr.eqns[0].params["n_consts"] == 0
 
     def test_dynamic_control_wires(self):
         """Test that control wires can be dynamic."""
@@ -366,7 +363,6 @@ class TestCtrlQfunc:
         plxpr = jax.make_jaxpr(f)(-0.5, 1, 2)
 
         # First equation of plxpr is the multiplication of x by 2
-        assert plxpr.eqns[1].params["n_consts"] == 1  # w1 is a const for the outer `ctrl`
         assert (
             plxpr.eqns[1].invars[0] is plxpr.jaxpr.invars[1]
         )  # first input is first control wire, const
@@ -409,7 +405,6 @@ class TestCtrlQfunc:
 
         eqn = jaxpr.eqns[0]
         assert eqn.params["control_values"] == [True]
-        assert eqn.params["n_consts"] == 0
         assert eqn.params["n_control"] == 1
         assert eqn.params["work_wires"] is None
 
diff --git a/tests/capture/test_promote_consts.py b/tests/capture/test_promote_consts.py
new file mode 100644
index 000000000..9ba364967
--- /dev/null
+++ b/tests/capture/test_promote_consts.py
@@ -0,0 +1,44 @@
+# Copyright 2018-2025 Xanadu Quantum Technologies Inc.
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+# pylint: disable=protected-access
+"""Unit tests for qml.capture.promote_consts."""
+
+import pytest
+
+from pennylane.capture import promote_consts
+
+jax = pytest.importorskip("jax")
+jnp = pytest.importorskip("jax.numpy")
+
+
+@pytest.mark.jax
+def test_promote_consts():
+    """Test that promote consts works as expected."""
+
+    const = jnp.array([2.0])
+
+    def f(x):
+        return const + x
+
+    args = (jnp.array([0.5]),)
+
+    jaxpr = jax.make_jaxpr(f)(*args)
+
+    new_jaxpr, new_args = promote_consts(jaxpr, args)
+
+    assert isinstance(new_jaxpr, jax.extend.core.Jaxpr)
+    assert new_jaxpr.constvars == []
+
+    assert len(new_args) == 2
+    assert jnp.allclose(new_args[0], const)
+    assert jnp.allclose(new_args[1], args[0])
diff --git a/tests/capture/transforms/test_capture_defer_measurements.py b/tests/capture/transforms/test_capture_defer_measurements.py
index 7ca44c736..c7d8d01f9 100644
--- a/tests/capture/transforms/test_capture_defer_measurements.py
+++ b/tests/capture/transforms/test_capture_defer_measurements.py
@@ -675,7 +675,7 @@ class TestDeferMeasurementsHigherOrderPrimitives:
 
         inner_jaxpr = jaxpr.eqns[0].params["qfunc_jaxpr"]
         collector = CollectOpsandMeas()
-        collector.eval(inner_jaxpr, jaxpr.consts, x)
+        collector.eval(inner_jaxpr, [], *jaxpr.consts, x)
 
         ops = collector.state["ops"]
         expected_ops = [
@@ -727,7 +727,10 @@ class TestDeferMeasurementsHigherOrderPrimitives:
         qfunc_jaxpr = inner_jaxpr.eqns[0].params["qfunc_jaxpr"]
 
         collector = CollectOpsandMeas()
-        collector.eval(qfunc_jaxpr, jaxpr.consts, x)
+        if postselect is None:
+            collector.eval(qfunc_jaxpr, [], x)
+        else:
+            collector.eval(qfunc_jaxpr, [], qml.math.array([postselect]), x)
 
         ops = collector.state["ops"]
         expected_ops = [qml.RX(x, 0), qml.CNOT([0, 3]), qml.CRX(x, [3, 0])]
diff --git a/tests/capture/transforms/test_capture_merge_amplitude_embedding.py b/tests/capture/transforms/test_capture_merge_amplitude_embedding.py
index 691d9363e..2f8ae87ef 100644
--- a/tests/capture/transforms/test_capture_merge_amplitude_embedding.py
+++ b/tests/capture/transforms/test_capture_merge_amplitude_embedding.py
@@ -622,9 +622,10 @@ class TestHigherOrderPrimitiveIntegration:
         assert jaxpr.eqns[2].primitive == qml.RY._primitive
 
         inner_jaxpr = jaxpr.eqns[1].params["jaxpr"]
-        assert inner_jaxpr.eqns[0].primitive == qml.AmplitudeEmbedding._primitive
-        assert qml.math.allclose(inner_jaxpr.eqns[0].params["n_wires"], 2)
-        assert inner_jaxpr.eqns[1].primitive == qml.X._primitive
+
+        assert inner_jaxpr.eqns[-2].primitive == qml.AmplitudeEmbedding._primitive
+        assert qml.math.allclose(inner_jaxpr.eqns[-2].params["n_wires"], 2)
+        assert inner_jaxpr.eqns[-1].primitive == qml.X._primitive
 
     @pytest.mark.parametrize("lazy", [True, False])
     def test_adjoint_transform_prim(self, lazy):
@@ -645,10 +646,11 @@ class TestHigherOrderPrimitiveIntegration:
         assert jaxpr.eqns[0].params["lazy"] == lazy
 
         inner_jaxpr = jaxpr.eqns[0].params["jaxpr"]
-        assert len(inner_jaxpr.eqns) == 2
-        assert inner_jaxpr.eqns[0].primitive == qml.AmplitudeEmbedding._primitive
-        assert qml.math.allclose(inner_jaxpr.eqns[0].params["n_wires"], 2)
-        assert inner_jaxpr.eqns[1].primitive == qml.X._primitive
+
+        assert len(inner_jaxpr.eqns) == 9
+        assert inner_jaxpr.eqns[-2].primitive == qml.AmplitudeEmbedding._primitive
+        assert qml.math.allclose(inner_jaxpr.eqns[-2].params["n_wires"], 2)
+        assert inner_jaxpr.eqns[-1].primitive == qml.X._primitive
 
     def test_cond_prim_only_true_branch(self):
         """Test that the transform works correctly when applied with cond_prim with just a true branch."""
@@ -814,12 +816,12 @@ class TestHigherOrderPrimitiveIntegration:
 
         assert jaxpr.eqns[0].primitive == qnode_prim
         qfunc_jaxpr = jaxpr.eqns[0].params["qfunc_jaxpr"]
-        assert len(qfunc_jaxpr.eqns) == 4
-        assert qfunc_jaxpr.eqns[0].primitive == qml.AmplitudeEmbedding._primitive
-        assert qml.math.allclose(qfunc_jaxpr.eqns[0].params["n_wires"], 2)
-        assert qfunc_jaxpr.eqns[1].primitive == qml.Hadamard._primitive
-        assert qfunc_jaxpr.eqns[2].primitive == qml.PauliZ._primitive
-        assert qfunc_jaxpr.eqns[3].primitive == qml.measurements.ExpectationMP._obs_primitive
+        assert len(qfunc_jaxpr.eqns) == 11
+        assert qfunc_jaxpr.eqns[-4].primitive == qml.AmplitudeEmbedding._primitive
+        assert qml.math.allclose(qfunc_jaxpr.eqns[-4].params["n_wires"], 2)
+        assert qfunc_jaxpr.eqns[-3].primitive == qml.Hadamard._primitive
+        assert qfunc_jaxpr.eqns[-2].primitive == qml.PauliZ._primitive
+        assert qfunc_jaxpr.eqns[-1].primitive == qml.measurements.ExpectationMP._obs_primitive
 
     def test_grad_prim(self):
         """Test that the transform works correctly when applied with grad_prim."""
diff --git a/tests/capture/workflow/test_capture_qnode.py b/tests/capture/workflow/test_capture_qnode.py
index d86d24b26..7641a3d89 100644
--- a/tests/capture/workflow/test_capture_qnode.py
+++ b/tests/capture/workflow/test_capture_qnode.py
@@ -318,7 +318,6 @@ def test_qnode_closure_variables():
 
     jaxpr = jax.make_jaxpr(circuit)(1)
     assert len(jaxpr.eqns[0].invars) == 2  # one closure variable, one arg
-    assert jaxpr.eqns[0].params["n_consts"] == 1
 
     out = jax.core.eval_jaxpr(jaxpr.jaxpr, [jnp.array(0.5)], 0)
     assert qml.math.allclose(out, jnp.cos(0.5))
@@ -630,6 +629,89 @@ class TestDevicePreprocessing:
         assert all(sample == 0 for sample in results) or all(sample == 1 for sample in results)
 
 
+def _interpret1(jaxpr, consts, *args):
+    args = [jnp.array([2]) * a for a in args]
+    return qml.capture.PlxprInterpreter().eval(jaxpr, consts, *args)
+
+
+# pylint: disable=unused-argument
+def _dummy_transform1(jaxpr, consts, targs, tkwargs, *args):
+    return jax.make_jaxpr(partial(_interpret1, jaxpr, consts))(*args)
+
+
+@partial(qml.transform, plxpr_transform=_dummy_transform1)
+def transform1(tape):
+    raise NotImplementedError
+
+
+def _interpret2(jaxpr, consts, *args):
+    args = [jnp.array([3.0, 3.0]) * a for a in args]
+    return qml.capture.PlxprInterpreter().eval(jaxpr, consts, *args)
+
+
+# pylint: disable=unused-argument
+def _dummy_transform2(jaxpr, consts, targs, tkwargs, *args):
+    return jax.make_jaxpr(partial(_interpret2, jaxpr, consts))(*args)
+
+
+@partial(qml.transform, plxpr_transform=_dummy_transform2)
+def transform2(tape):
+    raise NotImplementedError
+
+
+class JaxprDevice(qml.devices.Device):
+
+    def __init__(self, wires, shots=None):
+        super().__init__(wires=wires, shots=shots)
+        self.jaxpr_history = []
+        self.consts_history = []
+        self.args_history = []
+
+    def preprocess_transforms(self, execution_config=None):
+        program = qml.transforms.core.TransformProgram()
+        program.add_transform(transform1)
+        program.add_transform(transform2)
+        return program
+
+    def execute(self, circuits, execution_config=None):
+        raise NotImplementedError
+
+    def eval_jaxpr(self, jaxpr, consts, *args, execution_config=None):
+        self.jaxpr_history.append(jaxpr)
+        self.consts_history.append(consts)
+        self.args_history.append(args)
+        return qml.device("default.qubit", wires=self.wires).eval_jaxpr(
+            jaxpr, consts, *args, execution_config=execution_config
+        )
+
+
+def test_device_preprocessing_consts():
+    """Test that device preprocessing can add consts to the jaxpr multiple times."""
+
+    dev = JaxprDevice(wires=2)
+
+    @qml.qnode(dev)
+    def c(x):
+        qml.RX(x, 0)
+        return qml.expval(qml.Z(0))
+
+    _ = c(jnp.array(0.5))
+
+    jaxpr = dev.jaxpr_history[-1]
+    assert len(jaxpr.constvars) == 0
+    assert len(jaxpr.invars) == 3
+
+    assert jaxpr.invars[0].aval.shape == (2,)  # [3.0, 3.0]
+    assert jaxpr.invars[1].aval.shape == (1,)  # [2]
+
+    assert len(dev.consts_history[-1]) == 0
+    assert len(dev.args_history[-1]) == 3
+
+    assert jnp.allclose(dev.args_history[-1][0], jnp.array([3.0, 3.0]))
+    assert jnp.allclose(dev.args_history[-1][1], jnp.array([2.0]))
+    assert jnp.allclose(dev.args_history[-1][2], jnp.array(0.5))
+
+
 class TestDifferentiation:
 
     def test_error_backprop_unsupported(self):
@@ -901,7 +983,7 @@ class TestQNodeVmapIntegration:
             qml.RX(const, wires=0)
             return qml.expval(qml.PauliZ(0))
 
-        with pytest.warns(UserWarning, match="Constant argument at index 0 is not scalar. "):
+        with pytest.warns(UserWarning, match="Argument at index 0 has size > 1 "):
             jax.make_jaxpr(jax.vmap(circuit))(jnp.array([0.1, 0.2]))
 
     def test_vmap_overriding_shots(self):
diff --git a/tests/ops/qubit/test_arithmetic_ops.py b/tests/ops/qubit/test_arithmetic_ops.py
index 8b701b0fa..58b3de4ae 100644
--- a/tests/ops/qubit/test_arithmetic_ops.py
+++ b/tests/ops/qubit/test_arithmetic_ops.py
@@ -22,6 +22,7 @@ import numpy as np
 import pytest
 
 import pennylane as qml
+from pennylane.ops.functions.assert_valid import _test_decomposition_rule
 from pennylane.wires import Wires
 
 label_data = [
@@ -410,38 +411,50 @@ class TestIntegerComparator:
         with pytest.raises(ValueError, match=expected_error_message):
             qml.IntegerComparator.compute_decomposition(value, wires=wires, geq=geq)
 
-    def test_decomposition(self):
-        """Test operator's ``compute_decomposition()`` method."""
-        with qml.queuing.AnnotatedQueue() as q1:
-            qml.IntegerComparator.compute_decomposition(2, wires=[0, 1, 2])
-        tape1 = qml.tape.QuantumScript.from_queue(q1)
-        assert all(isinstance(op, qml.MultiControlledX) for op in tape1.operations)
+    @pytest.mark.parametrize("geq", [True, False])
+    @pytest.mark.parametrize("value", list(range(9)))
+    @pytest.mark.parametrize("initial_state", list(itertools.product([0, 1], repeat=3)))
+    def test_decomposition(self, value, geq, initial_state):
+        """Tests the decomposition of the IntegerComparator operator."""
 
-        with qml.queuing.AnnotatedQueue() as q2:
-            qml.IntegerComparator.compute_decomposition(2, wires=[0, 1, 2], geq=False)
-        tape2 = qml.tape.QuantumScript.from_queue(q2)
-        assert all(isinstance(op, qml.MultiControlledX) for op in tape2.operations)
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.Projector(initial_state, wires=[0, 1, 2])
+            qml.IntegerComparator.compute_decomposition(value, wires=[0, 1, 2, 3], geq=geq)
+        tape = qml.tape.QuantumScript.from_queue(q)
 
-        num_wires = 3
-        num_workers = 1
-        dev = qml.device("default.qubit", wires=num_wires + num_workers)
+        with qml.queuing.AnnotatedQueue() as q:
+            qml.Projector(initial_state, wires=[0, 1, 2])
+            val = int("".join(str(i) for i in initial_state), 2)
+            if (geq and val >= value) or (not geq and val < value):
+                qml.PauliX(3)
+        equivalent_tape = qml.tape.QuantumScript.from_queue(q)
+
+        matrix = qml.matrix(tape, wire_order=[0, 1, 2, 3])
+        equivalent_matrix = qml.matrix(equivalent_tape, wire_order=[0, 1, 2, 3])
+        assert qml.math.allclose(matrix, equivalent_matrix)
+
+    @pytest.mark.parametrize("geq", [True, False])
+    @pytest.mark.parametrize("value", list(range(9)))
+    def test_decomposition_qfuncs(self, value, geq):
+        """Tests the qfunc decompositions of the IntegerComparator operator."""
+        op = qml.IntegerComparator(value, wires=[0, 1, 2, 3], geq=geq)
+        for rule in qml.list_decomps(qml.IntegerComparator):
+            _test_decomposition_rule(op, rule)
+
+    @pytest.mark.integration
+    @pytest.mark.usefixtures("enable_graph_decomposition")
+    def test_decomposition_graph_integration(self):
+        """Tests that the decomposition is correct under the new system."""
+
+        tape = qml.tape.QuantumScript([qml.IntegerComparator(42, wires=[0, 1, 2, 3, 4, 5, 6])])
+        [decomp], _ = qml.transforms.decompose(
+            tape,
+            gate_set={qml.X, qml.CNOT, qml.Hadamard, qml.T, "Adjoint(T)", qml.RX, qml.RY, qml.RZ},
+        )
 
-        @qml.qnode(dev)
-        def f(bitstring, tape, geq):
-            qml.BasisState(bitstring, wires=range(num_wires + num_workers))
-            qml.adjoint(qml.IntegerComparator(2, wires=(0, 1, 2), geq=geq), lazy=False)
-            for op in tape.operations:
-                op.queue()
-            return qml.probs(wires=range(num_wires + num_workers))
-
-        for tape, geq in zip([tape1, tape2], [True, False]):
-            u = np.array(
-                [
-                    f(np.array(b), tape, geq)
-                    for b in itertools.product(range(2), repeat=num_wires + num_workers)
-                ]
-            ).T
-            assert np.allclose(u, np.eye(2 ** (num_wires + num_workers)))
+        mat = qml.matrix(decomp, wire_order=[0, 1, 2, 3, 4, 5, 6])
+        expected_mat = qml.matrix(tape, wire_order=[0, 1, 2, 3, 4, 5, 6])
+        assert qml.math.allclose(mat, expected_mat)
 
     def test_decomposition_extraneous_value(self):
         """Test operator's ``compute_decomposition()`` method when ``value`` is such that
diff --git a/tests/templates/test_subroutines/test_select_pauli_rot.py b/tests/templates/test_subroutines/test_select_pauli_rot.py
index 1fcb510df..281f1cdbe 100644
--- a/tests/templates/test_subroutines/test_select_pauli_rot.py
+++ b/tests/templates/test_subroutines/test_select_pauli_rot.py
@@ -21,6 +21,7 @@ import pytest
 
 import pennylane as qml
 from pennylane import numpy as np
+from pennylane.ops.functions.assert_valid import _test_decomposition_rule
 
 
 def get_tape(angles, wires):
@@ -53,7 +54,7 @@ class TestSelectPauliRot:
             rot_axis="X",
         )
 
-        qml.ops.functions.assert_valid(op)
+        qml.ops.functions.assert_valid(op, heuristic_resources=True)
 
     @pytest.mark.parametrize(
         ("angles", "rot_axis", "target_wire", "msg_match"),
@@ -163,6 +164,16 @@ class TestSelectPauliRot:
             for gate in dec[1::2]:
                 assert gate.name == "CNOT"
 
+    @pytest.mark.parametrize("n", [1, 2, 3, 4])
+    @pytest.mark.parametrize("axis", "XYZ")
+    def test_decomposition_new(self, n, axis):
+        """Tests that the decomposition of SelectPauliRot is correct with the new system."""
+
+        angles = np.random.random((2**n,))
+        op = qml.SelectPauliRot(angles, control_wires=range(n), target_wire=n, rot_axis=axis)
+        for rule in qml.list_decomps("SelectPauliRot"):
+            _test_decomposition_rule(op, rule)
+
     @pytest.mark.jax
     def test_interface_jax(self):
         """Test that SelectPauliRot works with jax"""
diff --git a/tests/test_qnode.py b/tests/test_qnode.py
index ceb261d00..d801fc771 100644
--- a/tests/test_qnode.py
+++ b/tests/test_qnode.py
@@ -229,7 +229,7 @@ class TestInitialization:
         def f():
             return qml.state()
 
-        assert f.execute_kwargs["cache"] is False
+        assert f.execute_kwargs["cache"] == "auto"
 
     def test_cache_initialization_maxdiff_2(self):
         """Test that when max_diff = 2, the cache initialization to True."""
@@ -238,7 +238,7 @@ class TestInitialization:
         def f():
             return qml.state()
 
-        assert f.execute_kwargs["cache"] is True
+        assert f.execute_kwargs["cache"] == "auto"
 
 
 # pylint: disable=too-many-public-methods
@@ -1394,6 +1394,20 @@ class TestShots:
         with pytest.warns(UserWarning, match="Cached execution with finite shots detected"):
             qml.execute([tape], dev, None, cache=cache)
 
+        cache2 = {}
+        with pytest.warns(UserWarning, match="Cached execution with finite shots detected"):
+            qml.execute([tape, tape], dev, cache=cache2)
+
+    def test_no_caching_by_default(self):
+        """Test that caching is turned off by default."""
+
+        dev = qml.device("default.qubit", wires=1)
+
+        tape = qml.tape.QuantumScript([qml.H(0)], [qml.sample(wires=0)], shots=1000)
+
+        res1, res2 = qml.execute([tape, tape], dev)
+        assert not np.allclose(res1, res2)
+
     def test_no_warning_infinite_shots(self):
         """Tests that no warning is raised when caching is used with infinite shots."""
         dev = qml.device("default.qubit", wires=1)
diff --git a/tests/transforms/test_cliffordt_transform.py b/tests/transforms/test_cliffordt_transform.py
index 0fe38aa4c..d75830921 100644
--- a/tests/transforms/test_cliffordt_transform.py
+++ b/tests/transforms/test_cliffordt_transform.py
@@ -21,6 +21,8 @@ import pytest
 import pennylane as qml
 from pennylane.transforms.decompositions.clifford_t_transform import (
     _CLIFFORD_T_GATES,
+    _CachedCallable,
+    _map_wires,
     _merge_param_gates,
     _one_qubit_decompose,
     _rot_decompose,
@@ -200,7 +202,9 @@ class TestCliffordCompile:
         dev = qml.device("default.qubit")
 
         qnode_basic = qml.QNode(circuit, dev)
-        qnode_transformed = clifford_t_decomposition(qnode_basic, epsilon=epsilon, max_depth=10)
+        qnode_transformed = clifford_t_decomposition(
+            qnode_basic, epsilon=epsilon, max_depth=10, basis_set=("T", "T*", "H")
+        )
         mat_exact = qml.matrix(qnode_basic, wire_order=[0, 1])()
         mat_approx = qml.matrix(qnode_transformed, wire_order=[0, 1])()
         phase = qml.math.divide(
@@ -458,7 +462,7 @@ class TestCliffordCompile:
     def test_clifford_decompose_interfaces(self):
         """Test that unwrap converts lists to lists and interface variables to numpy."""
 
-        dev = qml.device("default.qubit", wires=3)
+        dev = qml.device("default.qubit")
 
         def circuit(x):
             qml.RZ(x[0], wires=[0])
@@ -510,3 +514,71 @@ class TestCliffordCompile:
         # Compare results
         assert all(qml.math.allclose(res1, res2, atol=1e-2) for res1, res2 in zip(*funres))
         assert all(qml.math.allclose(res1, res2, atol=1e-2) for res1, res2 in zip(*igrads))
+
+
+def circuit_7(num_repeat, rand_angles):
+    """Circuit 7 with a repeated operations"""
+    for angle in rand_angles:
+        for idx in range(num_repeat):
+            qml.RZ(angle, idx)
+        for idx in range(num_repeat):
+            qml.CNOT([idx, (idx + 1) % num_repeat])
+    return qml.expval(qml.Z(0))
+
+
+class TestCliffordCached:
+    """Unit tests for clifford caching function."""
+
+    # pylint: disable=protected-access, import-outside-toplevel, reimported
+    def test_clifford_cached(self):
+        """Test that the cached version of the circuit is equivalent to the original one."""
+
+        import pennylane.transforms.decompositions.clifford_t_transform as clt2
+
+        clt2._CLIFFORD_T_CACHE = None
+
+        num_angles = 1
+        rand_angles = qml.math.random.random.rand(num_angles)
+        rand_angles = qml.math.concatenate((rand_angles, -rand_angles))
+
+        num_repeat = 2
+        old_tape = qml.tape.make_qscript(circuit_7)(num_repeat, rand_angles)
+        _ = clifford_t_decomposition(old_tape, epsilon=10)
+
+        assert isinstance(clt2._CLIFFORD_T_CACHE, _CachedCallable)
+        cache_info = clt2._CLIFFORD_T_CACHE.query.cache_info()
+        assert cache_info.misses == 2 * num_angles
+        assert cache_info.hits == 2 * num_angles * (num_repeat - 1)
+
+        num_repeat = 2
+        old_tape = qml.tape.make_qscript(circuit_7)(num_repeat, rand_angles)
+        _ = clifford_t_decomposition(old_tape, epsilon=10)
+
+        assert isinstance(clt2._CLIFFORD_T_CACHE, _CachedCallable)
+        cache_info = clt2._CLIFFORD_T_CACHE.query.cache_info()
+        assert cache_info.misses == 2 * num_angles
+        assert cache_info.hits == 2 * num_angles * (2 * num_repeat - 1)
+
+        num_repeat = 2
+        old_tape = qml.tape.make_qscript(circuit_7)(num_repeat, rand_angles)
+        _ = clifford_t_decomposition(old_tape, cache_size=100)
+
+        assert isinstance(clt2._CLIFFORD_T_CACHE, _CachedCallable)
+        cache_info = clt2._CLIFFORD_T_CACHE.query.cache_info()
+        assert cache_info.misses == 2 * num_angles
+        assert cache_info.hits == 2 * num_angles * (num_repeat - 1)
+        assert cache_info.maxsize == 100
+
+    def test_wire_mapping(self):
+        """Test that wire mapping is being cached correctly."""
+        _map_wires.cache_clear()  # Clear the cache before testing
+
+        for wire in range(5):
+            assert _map_wires(qml.X(0), wire) == qml.X(wire)
+        assert _map_wires.cache_info().hits == 0
+        assert _map_wires.cache_info().misses == 5
+
+        for wire in range(10):
+            assert _map_wires(qml.X(0), wire) == qml.X(wire)
+        assert _map_wires.cache_info().hits == 5
+        assert _map_wires.cache_info().misses == 10
diff --git a/tests/workflow/interfaces/execute/test_jax_jit.py b/tests/workflow/interfaces/execute/test_jax_jit.py
index bb3285a87..c8aa45fc3 100644
--- a/tests/workflow/interfaces/execute/test_jax_jit.py
+++ b/tests/workflow/interfaces/execute/test_jax_jit.py
@@ -195,6 +195,7 @@ class TestCaching:
                 [tape],
                 dev,
                 diff_method=param_shift,
+                cache=True,
                 cachesize=cachesize,
             )[0]
 
diff --git a/tests/workflow/interfaces/test_transform_program_integration.py b/tests/workflow/interfaces/test_transform_program_integration.py
index c346ea868..bf0dc9b5c 100644
--- a/tests/workflow/interfaces/test_transform_program_integration.py
+++ b/tests/workflow/interfaces/test_transform_program_integration.py
@@ -156,7 +156,7 @@ class TestTransformProgram:
         )
         with dev.tracker:
             results = qml.execute(
-                (tape1, tape2, tape3), dev, interface=interface, transform_program=prog
+                (tape1, tape2, tape3), dev, interface=interface, transform_program=prog, cache=True
             )
 
         assert qml.math.allclose(results[0], np.cos(1.2))
diff --git a/tests/workflow/test_cache_transform.py b/tests/workflow/test_cache_transform.py
index a270387f1..41087a40e 100644
--- a/tests/workflow/test_cache_transform.py
+++ b/tests/workflow/test_cache_transform.py
@@ -101,13 +101,13 @@ def test_batch_of_identical_tapes(cache):
     assert batch_fns((result,)) == (result, result, result)
 
 
-def test_finite_shots_with_persistent_cache_warning(cache):
+def test_finite_shots_with_persistent_cache_warning():
     """Tests that a UserWarning is emitted if a cache hit occurs for a tape with
-    finite shots that uses a persistent cache.
+    finite shots that uses a cache.
     """
     tape = QuantumScript([], [qml.expval(qml.Z(0))], shots=1)
 
-    batch_tapes, batch_fns = _cache_transform([tape, tape], cache=cache)
+    batch_tapes, batch_fns = _cache_transform([tape, tape], cache={})
     assert len(batch_tapes) == 1
 
     with pytest.warns(UserWarning, match=r"Cached execution with finite shots detected!"):