diff --git a/demonstrations_v2/tutorial_liealgebra/demo.py b/demonstrations_v2/tutorial_liealgebra/demo.py
index 49566c0140..af9a664483 100644
--- a/demonstrations_v2/tutorial_liealgebra/demo.py
+++ b/demonstrations_v2/tutorial_liealgebra/demo.py
@@ -98,22 +98,22 @@
 First, let us do a linear combination of :math:`\{iX, iY, iZ\}` with some real values and check unitarity after putting them in the exponent.
 """
 import numpy as np
-import pennylane as qml
+import pennylane as qp
 from pennylane import X, Y, Z
 
 su2 = [1j * X(0), 1j * Y(0), 1j * Z(0)]
 
 coeffs = [1., 2., 3.]                           # some real coefficients
-exponent = qml.dot(coeffs, su2)                 # linear combination of operators
-U = qml.math.expm(exponent.matrix())            # compute matrix exponent of lin. comb.
+exponent = qp.dot(coeffs, su2)                 # linear combination of operators
+U = qp.math.expm(exponent.matrix())            # compute matrix exponent of lin. comb.
 print(np.allclose(U.conj().T @ U, np.eye(2)))   # check that result is unitary UU* = 1
 
 ##############################################################################
 # If we throw complex values in the mix, the resulting matrix is not unitary anymore.
 
 coeffs = [1., 2.+ 1j, 3.]                       # some complex coefficients
-exponent = qml.dot(coeffs, su2)
-U = qml.math.expm(exponent.matrix())
+exponent = qp.dot(coeffs, su2)
+U = qp.math.expm(exponent.matrix())
 print(np.allclose(U.conj().T @ U, np.eye(2)))   # result is not unitary anymore
 
 ##############################################################################
@@ -158,7 +158,7 @@
 #
 # Let us do a quick example and compute the Lie closure of :math:`\{iX, iY\}` (more examples later).
 
-print(qml.commutator(1j * X(0), 1j * Y(0)))
+print(qp.commutator(1j * X(0), 1j * Y(0)))
 
 ##############################################################################
 # We know that the commutator between :math:`iX` and :math:`iY` yields a new operator :math:`\propto iZ.`
@@ -168,15 +168,15 @@
 list_ops = [1j * X(0), 1j * Y(0), 1j * Z(0)]
 for op1 in list_ops:
     for op2 in list_ops:
-        print(qml.commutator(op1, op2))
+        print(qp.commutator(op1, op2))
 
 ##############################################################################
 # Since no new operators have been created we know the lie closure is complete and our dynamical Lie algebra
 # is :math:`\langle\{iX, iY\}\rangle_\text{Lie} = \{iX, iY, iZ\}( = \mathfrak{su}(2)).` 
 # 
-# PennyLane provides some dedicated functionality for Lie algebras. We can compute the Lie closure of the generators using ``qml.lie_closure``.
+# PennyLane provides some dedicated functionality for Lie algebras. We can compute the Lie closure of the generators using ``qp.lie_closure``.
 
-dla = qml.lie_closure([X(0), Y(0)])
+dla = qp.lie_closure([X(0), Y(0)])
 dla
 
 ##############################################################################
@@ -187,14 +187,14 @@
 # these :math:`X` and :math:`Y` rotations :math:`e^{-i \phi X}` and :math:`e^{-i \phi Y}.` 
 # For example, let us write a Pauli-Z rotation at non-trivial angle :math:`0.5` as a product of them.
 
-U_target = qml.matrix(qml.RZ(-0.5, 0))
-decomp = qml.ops.one_qubit_decomposition(U_target, 0, rotations="XYX")
+U_target = qp.matrix(qp.RZ(-0.5, 0))
+decomp = qp.ops.one_qubit_decomposition(U_target, 0, rotations="XYX")
 print(decomp)
 
 ##############################################################################
 # We can check that this is indeed a valid decomposition by computing the trace distance to the target.
 
-U = qml.prod(*decomp).matrix()
+U = qp.prod(*decomp).matrix()
 print(1 - np.real(np.trace(U_target @ U))/2)
 
 ##############################################################################
@@ -215,10 +215,10 @@
 generators = [1j * (X(0) @ X(1)), 1j * Z(0), 1j * Z(1)]
 
 # collection of linearly independent basis vectors, automatically discards linearly dependent ones
-dla = qml.pauli.PauliVSpace(generators, dtype=complex)
+dla = qp.pauli.PauliVSpace(generators, dtype=complex)
 for i, op1 in enumerate(generators):
     for op2 in generators[i+1:]:
-        res = qml.commutator(op1, op2)/2
+        res = qp.commutator(op1, op2)/2
         res = res.simplify() # ensures all products of scalars are executed
         print(f"[{op1}, {op2}] = {res}")
 
@@ -234,7 +234,7 @@
 
 for i, op1 in enumerate(dla.basis.copy()):
     for op2 in dla.basis.copy()[i+1:]:
-        res = qml.commutator(op1, op2)/2
+        res = qp.commutator(op1, op2)/2
         res = res.simplify()
         print(f"[{op1}, {op2}] = {res}")
 
@@ -257,7 +257,7 @@
 # We have constructed the DLA by hand to showcase the process. We can use the PennyLane function :func:`~lie_closure` for convenience.
 # In that case, we omit the explicit use of the imgaginary factor.
 
-dla2 = qml.lie_closure([X(0) @ X(1), Z(0), Z(1)])
+dla2 = qp.lie_closure([X(0) @ X(1), Z(0), Z(1)])
 for op in dla2:
     print(op)
 
@@ -277,8 +277,8 @@ def IsingGenerators(n, bc="open"):
     return gens
 
 for n in range(2, 5):
-    open_ = qml.lie_closure(IsingGenerators(n, "open"))
-    periodic_ = qml.lie_closure(IsingGenerators(n, "periodic"))
+    open_ = qp.lie_closure(IsingGenerators(n, "open"))
+    periodic_ = qp.lie_closure(IsingGenerators(n, "periodic"))
     print(f"Ising for n = {n}")
     print(f"open: {len(open_)} = {n*(2*n-1)} = 2n * (2n - 1)/2")
     print(f"open: {len(periodic_)} = {2*n*(2*n-1)} = 2 * 2n * (2n - 1)/2")
@@ -328,15 +328,15 @@ def IsingGenerators(n, bc="open"):
 # Let us briefly verify this for a small example for ``n = 3`` qubits that readily generalizes to arbitrary sizes.
 
 n = 3
-H = qml.sum(*(P(i) @ P(i+1) for i in range(n-1) for P in [X, Y, Z]))
+H = qp.sum(*(P(i) @ P(i+1) for i in range(n-1) for P in [X, Y, Z]))
 
-SX = qml.sum(*(X(i) for i in range(n)))
-SY = qml.sum(*(Y(i) for i in range(n)))
-SZ = qml.sum(*(Z(i) for i in range(n)))
+SX = qp.sum(*(X(i) for i in range(n)))
+SY = qp.sum(*(Y(i) for i in range(n)))
+SZ = qp.sum(*(Z(i) for i in range(n)))
 
-print(qml.commutator(H, SX))
-print(qml.commutator(H, SY))
-print(qml.commutator(H, SZ))
+print(qp.commutator(H, SX))
+print(qp.commutator(H, SY))
+print(qp.commutator(H, SZ))
 
 ##############################################################################
 #
@@ -363,9 +363,9 @@ def IsingGenerators(n, bc="open"):
 # This is easily verified by looking at the commutation relation between these operators that match :math:`[\hat{O}_i, \hat{O}_j] = 2i \varepsilon_{ij\ell} \hat{O}_\ell,` the defining
 # property of :math:`\mathfrak{su}(2).`
 
-print(qml.commutator(SX, SY) == (2j*SZ).simplify())
-print(qml.commutator(SZ, SX) == (2j*SY).simplify())
-print(qml.commutator(SY, SZ) == (2j*SX).simplify())
+print(qp.commutator(SX, SY) == (2j*SZ).simplify())
+print(qp.commutator(SZ, SX) == (2j*SY).simplify())
+print(qp.commutator(SY, SZ) == (2j*SX).simplify())
 
 ##############################################################################
 #
diff --git a/demonstrations_v2/tutorial_liealgebra/metadata.json b/demonstrations_v2/tutorial_liealgebra/metadata.json
index 34d266c35b..8f1a5ee484 100644
--- a/demonstrations_v2/tutorial_liealgebra/metadata.json
+++ b/demonstrations_v2/tutorial_liealgebra/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-02-27T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing",
         "Getting Started"
diff --git a/demonstrations_v2/tutorial_liesim/demo.py b/demonstrations_v2/tutorial_liesim/demo.py
index 1a04dbac6c..91ad534511 100644
--- a/demonstrations_v2/tutorial_liesim/demo.py
+++ b/demonstrations_v2/tutorial_liesim/demo.py
@@ -156,7 +156,7 @@
 
 """
 
-import pennylane as qml
+import pennylane as qp
 from pennylane import X, Z, I
 import numpy as np
 
@@ -186,7 +186,7 @@
 # work with PauliSentence instances for efficiency
 generators = [op.pauli_rep for op in generators]
 
-dla = qml.lie_closure(generators, pauli=True)
+dla = qp.lie_closure(generators, pauli=True)
 dim_g = len(dla)
 
 ##############################################################################
@@ -207,7 +207,7 @@
     # initial state |0x0| = (I + Z)/2, note that trace function
     # below already normalizes by the dimension,
     # so we can ommit the explicit factor /2
-    rho_in = qml.prod(*(I(i) + Z(i) for i in h_i.wires))
+    rho_in = qp.prod(*(I(i) + Z(i) for i in h_i.wires))
     rho_in = rho_in.pauli_rep
 
     e_in[i] = (h_i @ rho_in).trace()
@@ -239,7 +239,7 @@
 # the forward pass of the expectation value computation. For demonstration purposes,
 # we choose a random subset of ``depth=10`` generators for gates from the DLA.
 
-adjoint_repr = qml.structure_constants(dla)
+adjoint_repr = qp.structure_constants(dla)
 
 depth = 10
 gate_choice = np.random.choice(dim_g, size=depth)
@@ -264,13 +264,13 @@ def forward(theta):
 ##############################################################################
 # As a sanity check, we compare the computation with the full state vector equivalent circuit.
 
-H = 0.5 * qml.sum(*[op.operation() for op in generators])
+H = 0.5 * qp.sum(*[op.operation() for op in generators])
 
-@qml.qnode(qml.device("default.qubit"), interface="jax")
+@qp.qnode(qp.device("default.qubit"), interface="jax")
 def qnode(theta):
     for i, mu in enumerate(gate_choice):
-        qml.exp(-1j * theta[i] * dla[mu].operation())
-    return qml.expval(H)
+        qp.exp(-1j * theta[i] * dla[mu].operation())
+    return qp.expval(H)
 
 statevec_forward, statevec_backward = qnode(theta), jax.grad(qnode)(theta)
 statevec_forward, statevec_backward
@@ -282,8 +282,8 @@ def qnode(theta):
 # expectation value vector.
 
 print(
-    qml.math.allclose(statevec_forward, gsim_forward), 
-    qml.math.allclose(statevec_backward, gsim_backward),
+    qp.math.allclose(statevec_forward, gsim_forward), 
+    qp.math.allclose(statevec_backward, gsim_backward),
 )
 
 ##############################################################################
diff --git a/demonstrations_v2/tutorial_liesim/metadata.json b/demonstrations_v2/tutorial_liesim/metadata.json
index 7b2680231c..dc43b78ca2 100644
--- a/demonstrations_v2/tutorial_liesim/metadata.json
+++ b/demonstrations_v2/tutorial_liesim/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-06-07T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing",
         "Getting Started"
diff --git a/demonstrations_v2/tutorial_liesim_extension/demo.py b/demonstrations_v2/tutorial_liesim_extension/demo.py
index 177f0cd807..56c1a40446 100644
--- a/demonstrations_v2/tutorial_liesim_extension/demo.py
+++ b/demonstrations_v2/tutorial_liesim_extension/demo.py
@@ -84,7 +84,7 @@
 
 """
 
-import pennylane as qml
+import pennylane as qp
 import numpy as np
 
 from pennylane import X, Y, Z, I
@@ -99,7 +99,7 @@ def TFIM(n):
     generators += [Z(i) for i in range(n)]
     generators = [op.pauli_rep for op in generators]
 
-    dla = qml.lie_closure(generators, pauli=True)
+    dla = qp.lie_closure(generators, pauli=True)
     dim_g = len(dla)
     return generators, dla, dim_g
 
@@ -109,7 +109,7 @@ def TFIM(n):
 # In regular :math:`\mathfrak{g}`-sim, the unitary evolution :math:`\mathcal{U}` of the expectation vector
 # is simply generated by the adjoint representation :math:`U.`
 
-adjoint_repr = qml.structure_constants(dla)
+adjoint_repr = qp.structure_constants(dla)
 gate = adjoint_repr[-1]
 theta = 0.5
 
@@ -135,7 +135,7 @@ def TFIM(n):
 
 p = dla[-5] @ dla[-1]
 p = next(iter(p)) # strip any scalar coefficients
-dla_vspace = qml.pauli.PauliVSpace(dla, dtype=complex)
+dla_vspace = qp.pauli.PauliVSpace(dla, dtype=complex)
 dla_vspace.is_independent(p.pauli_rep)
 
 
@@ -248,7 +248,7 @@ def exppw(theta, ps):
 E_in = [E0_in, E1_in]
 
 for i, hi in enumerate(dla):
-    rho_in = qml.prod(*(I(i) + Z(i) for i in hi.wires))
+    rho_in = qp.prod(*(I(i) + Z(i) for i in hi.wires))
     rho_in = rho_in.pauli_rep
 
     E_in[0][i] = (hi @ rho_in).trace()
@@ -256,8 +256,8 @@ def exppw(theta, ps):
 for i, hi in enumerate(dla):
     for j, hj in enumerate(dla):
         prod = hi @ hj
-        if prod.wires != qml.wires.Wires([]):
-            rho_in = qml.prod(*(I(i) + Z(i) for i in prod.wires))
+        if prod.wires != qp.wires.Wires([]):
+            rho_in = qp.prod(*(I(i) + Z(i) for i in prod.wires))
         else:
             rho_in = PauliWord({}) # identity
         rho_in = rho_in.pauli_rep
@@ -301,12 +301,12 @@ def exppw(theta, ps):
 ##############################################################################
 # As a sanity check, let us compare this to the same circuit but using our default state vector simulator in PennyLane.
 
-@qml.qnode(qml.device("default.qubit"))
+@qp.qnode(qp.device("default.qubit"))
 def true():
-    qml.exp(-1j * theta * dla[-1].operation())
-    qml.exp(-1j * 0.5 * p.operation())
-    qml.exp(-1j * 0.5 * dla[-2].operation())
-    return [qml.expval(op.operation()) for op in dla]
+    qp.exp(-1j * theta * dla[-1].operation())
+    qp.exp(-1j * 0.5 * p.operation())
+    qp.exp(-1j * 0.5 * dla[-2].operation())
+    return [qp.expval(op.operation()) for op in dla]
 
 true_res = np.array(true())
 
@@ -362,7 +362,7 @@ def true():
 # We now set up the moment vector spaces starting from the DLA and keep adding linearly independent product operators.
 
 def Moment_step(ops, dla):
-    MomentX = qml.pauli.PauliVSpace(ops.copy())
+    MomentX = qp.pauli.PauliVSpace(ops.copy())
     for i, op1 in enumerate(dla):
         for op2 in ops[i+1:]:
             prod = op1 @ op2
@@ -399,7 +399,7 @@ def Moment_step(ops, dla):
 e_in = np.zeros((dim[comp_moment]), dtype=float)
 
 for i, hi in enumerate(Moment[comp_moment]):
-    rho_in = qml.prod(*(I(i) + Z(i) for i in hi.wires))
+    rho_in = qp.prod(*(I(i) + Z(i) for i in hi.wires))
     rho_in = rho_in.pauli_rep
 
     e_in[i] = (hi @ rho_in).trace()
@@ -408,7 +408,7 @@ def Moment_step(ops, dla):
 #
 # Next, we compute the (pseudo-)adjoint representation of :math:`\mathcal{M}^1.`
 
-adjoint_repr = qml.structure_constants(Moment[comp_moment])
+adjoint_repr = qp.structure_constants(Moment[comp_moment])
 
 ##############################################################################
 #
@@ -454,7 +454,7 @@ def Moment_step(ops, dla):
         Moment.append(Moment_step(Moment[-1], dla))
         dim.append(len(Moment[-1]))
 
-    tempdla = qml.lie_closure(dla + [Moment[1][-1]], pauli=True)
+    tempdla = qp.lie_closure(dla + [Moment[1][-1]], pauli=True)
 
     dims_dla.append(dim_g)
     dims_moment.append(dim)
diff --git a/demonstrations_v2/tutorial_liesim_extension/metadata.json b/demonstrations_v2/tutorial_liesim_extension/metadata.json
index f20d933038..ce84d69bc5 100644
--- a/demonstrations_v2/tutorial_liesim_extension/metadata.json
+++ b/demonstrations_v2/tutorial_liesim_extension/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-06-18T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing",
         "Quantum Machine Learning"
diff --git a/demonstrations_v2/tutorial_local_cost_functions/demo.py b/demonstrations_v2/tutorial_local_cost_functions/demo.py
index 96f1400645..ea6a56eef1 100644
--- a/demonstrations_v2/tutorial_local_cost_functions/demo.py
+++ b/demonstrations_v2/tutorial_local_cost_functions/demo.py
@@ -56,7 +56,7 @@
 We first need to import the following modules.
 """
 
-import pennylane as qml
+import pennylane as qp
 from pennylane import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.ticker import LinearLocator, FormatStrFormatter
@@ -76,7 +76,7 @@
 # how many qubits we train on will effect our results.
 
 wires = 6
-dev = qml.device("lightning.qubit", wires=wires)
+dev = qp.device("lightning.qubit", wires=wires)
 
 
 ######################################################################
@@ -111,19 +111,19 @@
 
 def global_cost_simple(rotations):
     for i in range(wires):
-        qml.RX(rotations[0][i], wires=i)
-        qml.RY(rotations[1][i], wires=i)
-    return qml.probs(wires=range(wires))
+        qp.RX(rotations[0][i], wires=i)
+        qp.RY(rotations[1][i], wires=i)
+    return qp.probs(wires=range(wires))
 
 def local_cost_simple(rotations):
     for i in range(wires):
-        qml.RX(rotations[0][i], wires=i)
-        qml.RY(rotations[1][i], wires=i)
-    return [qml.probs(wires=i) for i in range(wires)]
+        qp.RX(rotations[0][i], wires=i)
+        qp.RY(rotations[1][i], wires=i)
+    return [qp.probs(wires=i) for i in range(wires)]
 
-global_circuit = qml.set_shots(qml.QNode(global_cost_simple, dev, interface="autograd"), shots = 10000)
+global_circuit = qp.set_shots(qp.QNode(global_cost_simple, dev, interface="autograd"), shots = 10000)
 
-local_circuit = qml.set_shots(qml.QNode(local_cost_simple, dev, interface="autograd"), shots = 10000)
+local_circuit = qp.set_shots(qp.QNode(local_cost_simple, dev, interface="autograd"), shots = 10000)
 
 def cost_local(rotations):
     return 1 - np.sum([i for (i, _) in local_circuit(rotations)]) / wires
@@ -150,12 +150,12 @@ def cost_global(rotations):
 print("Global Cost: {: .7f}".format(cost_global(rotations)))
 print("Local Cost: {: .7f}".format(cost_local(rotations)))
 
-qml.drawer.use_style('black_white')
-fig1, ax1 = qml.draw_mpl(global_circuit, decimals=2)(rotations)
+qp.drawer.use_style('black_white')
+fig1, ax1 = qp.draw_mpl(global_circuit, decimals=2)(rotations)
 fig1.suptitle("Global Cost", fontsize='xx-large')
 plt.show()
 
-fig2, ax2 = qml.draw_mpl(local_circuit, decimals=2)(rotations)
+fig2, ax2 = qp.draw_mpl(local_circuit, decimals=2)(rotations)
 fig2.suptitle("Local Cost", fontsize='xx-large')
 plt.show()
 
@@ -237,23 +237,23 @@ def plot_surface(surface):
 
 def global_cost_simple(rotations):
     for i in range(wires):
-        qml.RX(rotations[0][i], wires=i)
-        qml.RY(rotations[1][i], wires=i)
+        qp.RX(rotations[0][i], wires=i)
+        qp.RY(rotations[1][i], wires=i)
     for i in range(wires - 1):
-        qml.CNOT([i, i + 1])
-    return qml.probs(wires=range(wires))
+        qp.CNOT([i, i + 1])
+    return qp.probs(wires=range(wires))
 
 def local_cost_simple(rotations):
     for i in range(wires):
-        qml.RX(rotations[0][i], wires=i)
-        qml.RY(rotations[1][i], wires=i)
+        qp.RX(rotations[0][i], wires=i)
+        qp.RY(rotations[1][i], wires=i)
     for i in range(wires - 1):
-        qml.CNOT([i, i + 1])
-    return qml.probs(wires=[0])
+        qp.CNOT([i, i + 1])
+    return qp.probs(wires=[0])
 
-global_circuit = qml.set_shots(qml.QNode(global_cost_simple, dev, interface="autograd"), shots = 10000)
+global_circuit = qp.set_shots(qp.QNode(global_cost_simple, dev, interface="autograd"), shots = 10000)
 
-local_circuit = qml.set_shots(qml.QNode(local_cost_simple, dev, interface="autograd"), shots = 10000)
+local_circuit = qp.set_shots(qp.QNode(local_cost_simple, dev, interface="autograd"), shots = 10000)
 
 def cost_local(rotations):
     return 1 - local_circuit(rotations)[0]
@@ -285,7 +285,7 @@ def cost_global(rotations):
 
 
 rotations = np.array([[3.] * len(range(wires)), [0.] * len(range(wires))], requires_grad=True)
-opt = qml.GradientDescentOptimizer(stepsize=0.2)
+opt = qp.GradientDescentOptimizer(stepsize=0.2)
 steps = 100
 params_global = rotations
 for i in range(steps):
@@ -296,7 +296,7 @@ def cost_global(rotations):
         print("Cost after step {:5d}: {: .7f}".format(i + 1, cost_global(params_global)))
     if cost_global(params_global) < 0.1:
         break
-fig, ax = qml.draw_mpl(global_circuit, decimals=2)(params_global)
+fig, ax = qp.draw_mpl(global_circuit, decimals=2)(params_global)
 plt.show()
 
 
@@ -307,7 +307,7 @@ def cost_global(rotations):
 #
 
 rotations = np.array([[3. for i in range(wires)], [0. for i in range(wires)]], requires_grad=True)
-opt = qml.GradientDescentOptimizer(stepsize=0.2)
+opt = qp.GradientDescentOptimizer(stepsize=0.2)
 steps = 100
 params_local = rotations
 for i in range(steps):
@@ -319,7 +319,7 @@ def cost_global(rotations):
     if cost_local(params_local) < 0.05:
         break
 
-fig, ax = qml.draw_mpl(local_circuit, decimals=2)(params_local)
+fig, ax = qp.draw_mpl(local_circuit, decimals=2)(params_local)
 plt.show()
 
 
@@ -343,8 +343,8 @@ def cost_global(rotations):
 # us the exact representation.
 #
 
-_dev = qml.device("lightning.qubit", wires=wires)
-global_circuit = qml.QNode(global_cost_simple, _dev, interface="autograd")
+_dev = qp.device("lightning.qubit", wires=wires)
+global_circuit = qp.QNode(global_cost_simple, _dev, interface="autograd")
 print(
     "Current cost: "
     + str(cost_global(params_local))
@@ -369,24 +369,24 @@ def cost_global(rotations):
 
 def tunable_cost_simple(rotations):
     for i in range(wires):
-        qml.RX(rotations[0][i], wires=i)
-        qml.RY(rotations[1][i], wires=i)
+        qp.RX(rotations[0][i], wires=i)
+        qp.RY(rotations[1][i], wires=i)
     for i in range(wires - 1):
-        qml.CNOT([i, i + 1])
-    return qml.probs(range(locality))
+        qp.CNOT([i, i + 1])
+    return qp.probs(range(locality))
 
 def cost_tunable(rotations):
     return 1 - tunable_circuit(rotations)[0]
 
-tunable_circuit = qml.set_shots(qml.QNode(tunable_cost_simple, dev, interface="autograd"), shots = 10000)
+tunable_circuit = qp.set_shots(qp.QNode(tunable_cost_simple, dev, interface="autograd"), shots = 10000)
 locality = 2
 params_tunable = params_local
-fig, ax = qml.draw_mpl(tunable_circuit, decimals=2)(params_tunable)
+fig, ax = qp.draw_mpl(tunable_circuit, decimals=2)(params_tunable)
 plt.show()
 print(cost_tunable(params_tunable))
 
 locality = 2
-opt = qml.GradientDescentOptimizer(stepsize=0.1)
+opt = qp.GradientDescentOptimizer(stepsize=0.1)
 steps = 600
 for i in range(steps):
     # update the circuit parameters
@@ -406,7 +406,7 @@ def cost_tunable(rotations):
         continue
     elif runCost < 0.1 and locality >= wires:
         break
-fig, ax = qml.draw_mpl(tunable_circuit, decimals=2)(params_tunable)
+fig, ax = qp.draw_mpl(tunable_circuit, decimals=2)(params_tunable)
 plt.show()
 
 
@@ -436,12 +436,12 @@ def cost_tunable(rotations):
 samples = 10
 plateau = 0
 trained = 0
-opt = qml.GradientDescentOptimizer(stepsize=0.2)
+opt = qp.GradientDescentOptimizer(stepsize=0.2)
 steps = 400
 wires = 8
 
-dev = qml.device("lightning.qubit", wires=wires)
-global_circuit = qml.set_shots(qml.QNode(global_cost_simple, dev, interface="autograd"), shots = 10000)
+dev = qp.device("lightning.qubit", wires=wires)
+global_circuit = qp.set_shots(qp.QNode(global_cost_simple, dev, interface="autograd"), shots = 10000)
 
 for runs in range(samples):
     print("--- New run! ---")
@@ -469,12 +469,12 @@ def cost_tunable(rotations):
 samples = 10
 plateau = 0
 trained = 0
-opt = qml.GradientDescentOptimizer(stepsize=0.2)
+opt = qp.GradientDescentOptimizer(stepsize=0.2)
 steps = 400
 wires = 8
 
-dev = qml.device("lightning.qubit", wires=wires)
-tunable_circuit = qml.set_shots(qml.QNode(tunable_cost_simple, dev, interface="autograd"), shots = 10000)
+dev = qp.device("lightning.qubit", wires=wires)
+tunable_circuit = qp.set_shots(qp.QNode(tunable_cost_simple, dev, interface="autograd"), shots = 10000)
 
 for runs in range(samples):
     locality = 1
diff --git a/demonstrations_v2/tutorial_local_cost_functions/metadata.json b/demonstrations_v2/tutorial_local_cost_functions/metadata.json
index ad7db05171..858306c14b 100644
--- a/demonstrations_v2/tutorial_local_cost_functions/metadata.json
+++ b/demonstrations_v2/tutorial_local_cost_functions/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2020-09-09T00:00:00+00:00",
-    "dateOfLastModification": "2025-10-15T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Optimization"
     ],
diff --git a/demonstrations_v2/tutorial_loom_catalyst/demo.py b/demonstrations_v2/tutorial_loom_catalyst/demo.py
index af1a96fbc2..cd5e8bbce3 100644
--- a/demonstrations_v2/tutorial_loom_catalyst/demo.py
+++ b/demonstrations_v2/tutorial_loom_catalyst/demo.py
@@ -187,7 +187,7 @@
 
 from catalyst import qjit, cond, measure, debug
 from jax import random, numpy as jnp
-import pennylane as qml
+import pennylane as qp
 
 distance = 3
 
@@ -201,7 +201,7 @@
 # the Pennylane backend.
 # 
 
-dev = qml.device("qrack.simulator", wires=n_qubits)
+dev = qp.device("qrack.simulator", wires=n_qubits)
 
 ######################################################################
 # Then, we proceed to generate the circuit. Note that we start by introducing an error by hand, and
@@ -209,7 +209,7 @@
 # 
 
 @qjit()
-@qml.qnode(dev)
+@qp.qnode(dev)
 def circuit(seed : int):
 
     # Based on the seed, apply an X gate to a random data qubit
@@ -220,14 +220,14 @@ def circuit(seed : int):
     @cond(random_qubit != -1)
     def apply_noise():
         debug.print("Applying noise to qubit: {}", random_qubit)
-        qml.X(random_qubit)
+        qp.X(random_qubit)
     
     apply_noise()
 
     # Syndrome extraction routine: entangle data and auxiliary, and measure auxiliary
     for i in range(distance-1):
-        qml.CNOT(wires=[data_qubits[i], aux_qubits[i]])
-        qml.CNOT(wires=[data_qubits[i+1], aux_qubits[i]])
+        qp.CNOT(wires=[data_qubits[i], aux_qubits[i]])
+        qp.CNOT(wires=[data_qubits[i+1], aux_qubits[i]])
     syndrome = [measure(aux_qubit) for aux_qubit in aux_qubits]
 
 
@@ -235,15 +235,15 @@ def apply_noise():
     @cond(jnp.logical_and(syndrome[0] == 0, syndrome[1] == 1))
     def fix_data_qubits():
         debug.print("Applying correction on data qubit 2")
-        qml.X(data_qubits[2])
+        qp.X(data_qubits[2])
     @fix_data_qubits.else_if( jnp.logical_and(syndrome[0] == 1, syndrome[1] == 0))
     def fix_data_qubits():
         debug.print("Applying correction on data qubit 0")
-        qml.X(data_qubits[0])
+        qp.X(data_qubits[0])
     @fix_data_qubits.else_if(jnp.logical_and(syndrome[0] == 1, syndrome[1] == 1))
     def fix_data_qubits():
         debug.print("Applying correction on data qubit 1")
-        qml.X(data_qubits[1])
+        qp.X(data_qubits[1])
     
     # Apply the the fix
     fix_data_qubits()
@@ -348,7 +348,7 @@ def fix_data_qubits():
 # 
 
 # Define device
-dev = qml.device(
+dev = qp.device(
     "qrack.simulator",
     wires=len(reg_rc),
     shots=1,
@@ -376,7 +376,7 @@ def fix_data_qubits():
 # 
 
 @qjit()
-@qml.qnode(dev)
+@qp.qnode(dev)
 def repetition_code_circuit():
 
     ## Run the circuit
@@ -462,11 +462,11 @@ def repetition_code_circuit():
 # Indeed, the probability of getting a logical error without decoding is much higher. Interestingly
 # enough, though, even with the QEC scheme, the probability of error did not go to zero! Wonder why?
 # 
-# When we defined the ``qml.device``, we set the noise level at at 0.05:
+# When we defined the ``qp.device``, we set the noise level at at 0.05:
 # 
 
 # Define device
-dev = qml.device(
+dev = qp.device(
     "qrack.simulator",
     wires=len(reg_rc),
     shots=1,
diff --git a/demonstrations_v2/tutorial_loom_catalyst/metadata.json b/demonstrations_v2/tutorial_loom_catalyst/metadata.json
index 2a41af3d3c..f567cfc4cc 100644
--- a/demonstrations_v2/tutorial_loom_catalyst/metadata.json
+++ b/demonstrations_v2/tutorial_loom_catalyst/metadata.json
@@ -11,7 +11,7 @@
     "executable_stable": false,
     "executable_latest": false,
     "dateOfPublication": "2025-08-29T14:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing",
         "Quantum Hardware"
diff --git a/demonstrations_v2/tutorial_magic_state_distillation/demo.py b/demonstrations_v2/tutorial_magic_state_distillation/demo.py
index 0138e99d8c..bb48736361 100644
--- a/demonstrations_v2/tutorial_magic_state_distillation/demo.py
+++ b/demonstrations_v2/tutorial_magic_state_distillation/demo.py
@@ -61,7 +61,7 @@
 and a phase shift:
 """
 
-import pennylane as qml
+import pennylane as qp
 
 import jax; jax.config.update('jax_platform_name', 'cpu')
 from jax import numpy as jnp
@@ -70,8 +70,8 @@
 b = 0.5 * jnp.arccos(1 / jnp.sqrt(3))
 
 def generate_T0(wire):
-    qml.RY(2 * b, wires=wire)
-    qml.PhaseShift(jnp.pi / 4, wires=wire)
+    qp.RY(2 * b, wires=wire)
+    qp.PhaseShift(jnp.pi / 4, wires=wire)
 
 ######################################################################
 # Now, let’s create up a *faulty* :math:`| T_0\rangle` state-generating function. We can do this by
@@ -84,11 +84,11 @@ def generate_T0(wire):
 def faulty_generate_T0(wire, key, r):
     key, subkey = random.split(key)
     disturbance = random.uniform(subkey, minval=-r, maxval=r)
-    qml.RY(2 * (b + disturbance), wires=wire)
+    qp.RY(2 * (b + disturbance), wires=wire)
 
     key, subkey = random.split(key)
     disturbance = random.uniform(subkey, minval=-r, maxval=r)
-    qml.PhaseShift(jnp.pi / 4 + disturbance, wires=wire)
+    qp.PhaseShift(jnp.pi / 4 + disturbance, wires=wire)
     return key
 
 ######################################################################
@@ -128,30 +128,30 @@ def error_correction_decoder(wires):
     """
     w0, w1, w2, w3, w4 = wires
 
-    qml.CNOT(wires=[w1, w0])
-    qml.CZ(wires=[w1, w0])
-    qml.CZ(wires=[w1, w2])
-    qml.CZ(wires=[w1, w4])
+    qp.CNOT(wires=[w1, w0])
+    qp.CZ(wires=[w1, w0])
+    qp.CZ(wires=[w1, w2])
+    qp.CZ(wires=[w1, w4])
 
-    qml.CNOT(wires=[w2, w0])
-    qml.CZ(wires=[w2, w3])
-    qml.CZ(wires=[w2, w4])
+    qp.CNOT(wires=[w2, w0])
+    qp.CZ(wires=[w2, w3])
+    qp.CZ(wires=[w2, w4])
 
-    qml.CNOT(wires=[w3, w0])
+    qp.CNOT(wires=[w3, w0])
 
-    qml.CNOT(wires=[w4, w0])
-    qml.CZ(wires=[w4, w0])
+    qp.CNOT(wires=[w4, w0])
+    qp.CZ(wires=[w4, w0])
 
-    qml.PauliZ(wires=w0)
-    qml.PauliZ(wires=w1)
-    qml.PauliZ(wires=w4)
+    qp.PauliZ(wires=w0)
+    qp.PauliZ(wires=w1)
+    qp.PauliZ(wires=w4)
 
-    qml.Hadamard(wires=w1)
-    qml.Hadamard(wires=w2)
-    qml.Hadamard(wires=w3)
-    qml.Hadamard(wires=w4)
+    qp.Hadamard(wires=w1)
+    qp.Hadamard(wires=w2)
+    qp.Hadamard(wires=w3)
+    qp.Hadamard(wires=w4)
 
-print(qml.draw(error_correction_decoder)(range(5)))
+print(qp.draw(error_correction_decoder)(range(5)))
 
 ######################################################################
 # We’ll also define some helper functions to perform a conditional reset of a qubit into the
@@ -168,7 +168,7 @@ def measure_and_reset(wire):
     m = measure(wire)
 
     if m:
-        qml.PauliX(wires=wire)
+        qp.PauliX(wires=wire)
 
     return m
 
@@ -191,10 +191,10 @@ def measure_and_reset(wire):
 # operations.
 #
 
-dev = qml.device("lightning.qubit", wires=5)
+dev = qp.device("lightning.qubit", wires=5)
 
-@qml.qjit(autograph=True)
-@qml.qnode(dev)
+@qp.qjit(autograph=True)
+@qp.qnode(dev)
 def state_distillation(random_key, r):
     key = random_key
     syndrome = True
@@ -228,10 +228,10 @@ def state_distillation(random_key, r):
 
     # We can convert the T1 state back to T0
     # by applying a Hadamard and Pauli-Y rotation
-    qml.Hadamard(wires=0)
-    qml.PauliY(wires=0)
+    qp.Hadamard(wires=0)
+    qp.PauliY(wires=0)
 
-    return qml.state()
+    return qp.state()
 
 ######################################################################
 # Benchmark
@@ -246,19 +246,19 @@ def state_distillation(random_key, r):
 import matplotlib.pyplot as plt
 import os
 
-dev_5_qubits = qml.device("default.qubit", wires=5)
+dev_5_qubits = qp.device("default.qubit", wires=5)
 
 @jax.jit
-@qml.qnode(dev_5_qubits, interface="jax")
+@qp.qnode(dev_5_qubits, interface="jax")
 def T0_state():
     generate_T0(0)
-    return qml.state()
+    return qp.state()
 
 @jax.jit
-@qml.qnode(dev_5_qubits, interface="jax")
+@qp.qnode(dev_5_qubits, interface="jax")
 def faulty_T0_state(random_key, r):
     faulty_generate_T0(0, random_key, r)
-    return qml.state()
+    return qp.state()
 
 exact_T0_state = T0_state()
 perturbations = jnp.linspace(0, 1, 20)
diff --git a/demonstrations_v2/tutorial_magic_state_distillation/metadata.json b/demonstrations_v2/tutorial_magic_state_distillation/metadata.json
index 28ec7547df..d85ead16fd 100644
--- a/demonstrations_v2/tutorial_magic_state_distillation/metadata.json
+++ b/demonstrations_v2/tutorial_magic_state_distillation/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-04-26T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Algorithms",
         "Quantum Computing"
diff --git a/demonstrations_v2/tutorial_magic_states/demo.py b/demonstrations_v2/tutorial_magic_states/demo.py
index 53af85b59e..a11a89f9d3 100644
--- a/demonstrations_v2/tutorial_magic_states/demo.py
+++ b/demonstrations_v2/tutorial_magic_states/demo.py
@@ -47,36 +47,36 @@
 
 """
 
-import pennylane as qml
+import pennylane as qp
 from pennylane import numpy as np
 
 
 def prepare_magic_state():
     """Prepares the |H> magic state on wire 1."""
-    qml.Hadamard(wires=1)
-    qml.T(wires=1)
+    qp.Hadamard(wires=1)
+    qp.T(wires=1)
 
 
-dev = qml.device("default.qubit")
+dev = qp.device("default.qubit")
 
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def magic_state_injection_circuit(target_state_params):
 
     # Prepare the initial target state (e.g., Ry rotation)
-    qml.RY(target_state_params, wires=0)
+    qp.RY(target_state_params, wires=0)
 
     # Prepare the Magic State on Qubit 1
     prepare_magic_state()
 
     # Apply the Clifford operations for injection
-    qml.CNOT(wires=[0, 1])
+    qp.CNOT(wires=[0, 1])
 
     # The outcome of m_1 dictates the final correction
-    m_1 = qml.measure(1)
-    qml.cond(m_1 == 1, qml.S)(wires=0)
+    m_1 = qp.measure(1)
+    qp.cond(m_1 == 1, qp.S)(wires=0)
 
-    return qml.density_matrix(wires=[0])
+    return qp.density_matrix(wires=[0])
 
 
 print(magic_state_injection_circuit(np.pi / 3))
diff --git a/demonstrations_v2/tutorial_magic_states/metadata.json b/demonstrations_v2/tutorial_magic_states/metadata.json
index 5009bb2fd9..d1dcddc80b 100644
--- a/demonstrations_v2/tutorial_magic_states/metadata.json
+++ b/demonstrations_v2/tutorial_magic_states/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2026-01-23T10:00:00+00:00",
-    "dateOfLastModification": "2026-01-23T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing"
     ],
diff --git a/demonstrations_v2/tutorial_mapping/demo.py b/demonstrations_v2/tutorial_mapping/demo.py
index 7e7795f331..cf642d3bae 100644
--- a/demonstrations_v2/tutorial_mapping/demo.py
+++ b/demonstrations_v2/tutorial_mapping/demo.py
@@ -84,7 +84,7 @@
 then map the operator using :func:`~.pennylane.fermi.jordan_wigner`.
 """
 
-import pennylane as qml
+import pennylane as qp
 from pennylane.fermi import from_string, jordan_wigner
 
 qubits = 10
@@ -109,8 +109,8 @@
 
 orbitals = 10
 electrons = 5
-state_number = qml.qchem.hf_state(electrons, orbitals)
-state_parity = qml.qchem.hf_state(electrons, orbitals, basis="parity")
+state_number = qp.qchem.hf_state(electrons, orbitals)
+state_parity = qp.qchem.hf_state(electrons, orbitals, basis="parity")
 
 print("State in occupation number basis:\n", state_number)
 print("State in parity basis:\n", state_parity)
@@ -143,7 +143,7 @@
 # Parity mapping using :func:`~.pennylane.fermi.parity_transform` in PennyLane.
 
 qubits = 10
-pauli_pr = qml.parity_transform(fermi_op, qubits, ps=True)
+pauli_pr = qp.parity_transform(fermi_op, qubits, ps=True)
 pauli_pr
 
 ##############################################################################
@@ -155,11 +155,11 @@
 # `qubit tapering <https://pennylane.ai/qml/demos/tutorial_qubit_tapering>`__ demo.
 # Let's look at an example.
 
-generators = [qml.prod(*[qml.Z(i) for i in range(qubits-1)]), qml.Z(qubits-1)]
-paulixops = qml.paulix_ops(generators, qubits)
+generators = [qp.prod(*[qp.Z(i) for i in range(qubits-1)]), qp.Z(qubits-1)]
+paulixops = qp.paulix_ops(generators, qubits)
 paulix_sector = [1, 1]
-taper_op = qml.taper(pauli_pr, generators, paulixops, paulix_sector)
-qml.simplify(taper_op)
+taper_op = qp.taper(pauli_pr, generators, paulixops, paulix_sector)
+qp.simplify(taper_op)
 
 ###############################################################################
 # Note that the tapered operator doesn't include qubit :math:`8` and :math:`9.`
@@ -175,7 +175,7 @@
 # Let's use the :func:`~.pennylane.fermi.bravyi_kitaev` function to map our :math:`a_{5}^{\dagger}`
 # operator.
 
-pauli_bk = qml.bravyi_kitaev(fermi_op, qubits, ps=True)
+pauli_bk = qp.bravyi_kitaev(fermi_op, qubits, ps=True)
 pauli_bk
 
 ##############################################################################
@@ -225,7 +225,7 @@
 
 electrons = 2
 qubits = len(h_fermi.wires)
-h_pauli = qml.bravyi_kitaev(h_fermi, qubits, tol=1e-16)
+h_pauli = qp.bravyi_kitaev(h_fermi, qubits, tol=1e-16)
 
 ##############################################################################
 # Initial state
@@ -282,11 +282,11 @@
 
 singles_pauli = []
 for op in singles_fermi:
-    singles_pauli.append(qml.bravyi_kitaev(op, qubits, ps=True))
+    singles_pauli.append(qp.bravyi_kitaev(op, qubits, ps=True))
 
 doubles_pauli = []
 for op in doubles_fermi:
-    doubles_pauli.append(qml.bravyi_kitaev(op, qubits, ps=True))
+    doubles_pauli.append(qp.bravyi_kitaev(op, qubits, ps=True))
 
 ##############################################################################
 # Note that we need to exponentiate these operators to be able to use them in the circuit
@@ -294,19 +294,19 @@
 
 params = jnp.array([0.22347661, 0.0, 0.0])
 
-dev = qml.device("default.qubit", wires=qubits)
+dev = qp.device("default.qubit", wires=qubits)
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def circuit(params):
-    qml.BasisState(hf_state, wires=range(qubits))
+    qp.BasisState(hf_state, wires=range(qubits))
 
     for i, excitation in enumerate(doubles_pauli):
-        qml.exp((excitation * params[i] / 2).operation()), range(qubits)
+        qp.exp((excitation * params[i] / 2).operation()), range(qubits)
 
     for j, excitation in enumerate(singles_pauli):
-        qml.exp((excitation * params[i + j + 1] / 2).operation()), range(qubits)
+        qp.exp((excitation * params[i + j + 1] / 2).operation()), range(qubits)
 
-    return qml.expval(h_pauli)
+    return qp.expval(h_pauli)
 
 print('Energy =', circuit(params))
 
diff --git a/demonstrations_v2/tutorial_mapping/metadata.json b/demonstrations_v2/tutorial_mapping/metadata.json
index fbb0c1de4d..d6a1428ef3 100644
--- a/demonstrations_v2/tutorial_mapping/metadata.json
+++ b/demonstrations_v2/tutorial_mapping/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-05-06T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Algorithms",
         "Quantum Computing",
diff --git a/demonstrations_v2/tutorial_mbqc/demo.py b/demonstrations_v2/tutorial_mbqc/demo.py
index 06c3e71d2c..0a5aaba3e1 100644
--- a/demonstrations_v2/tutorial_mbqc/demo.py
+++ b/demonstrations_v2/tutorial_mbqc/demo.py
@@ -89,25 +89,25 @@
 # and define a circuit to prepare the cluster state.
 #
 
-import pennylane as qml
+import pennylane as qp
 
 qubits = [str(node) for node in G.nodes]
-dev = qml.device("lightning.qubit", wires=qubits)
+dev = qp.device("lightning.qubit", wires=qubits)
 
 
-@qml.qnode(dev, interface="autograd")
+@qp.qnode(dev, interface="autograd")
 def cluster_state():
     for node in qubits:
-        qml.Hadamard(wires=[node])
+        qp.Hadamard(wires=[node])
 
     for edge in G.edges:
         i, j = edge
-        qml.CZ(wires=[str(i), str(j)])
+        qp.CZ(wires=[str(i), str(j)])
 
-    return qml.state()
+    return qp.state()
 
 
-print(qml.draw(cluster_state)())
+print(qp.draw(cluster_state)())
 
 ##############################################################################
 #
@@ -151,30 +151,30 @@ def cluster_state():
 #
 # Let's implement one-qubit teleportation in PennyLane.
 
-import pennylane as qml
+import pennylane as qp
 import pennylane.numpy as np
 
-dev = qml.device("lightning.qubit", wires=2)
+dev = qp.device("lightning.qubit", wires=2)
 
 
-@qml.qnode(dev, interface="autograd")
+@qp.qnode(dev, interface="autograd")
 def one_bit_teleportation(input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=0)
+    qp.StatePrep(input_state, wires=0)
 
     # Prepare the cluster state
-    qml.Hadamard(wires=1)
-    qml.CZ(wires=[0, 1])
+    qp.Hadamard(wires=1)
+    qp.CZ(wires=[0, 1])
 
     # Measure the first qubit in the Pauli-X basis
     # and apply an X-gate conditioned on the outcome
-    qml.Hadamard(wires=0)
-    m = qml.measure(wires=[0])
-    qml.cond(m == 1, qml.PauliX)(wires=1)
-    qml.Hadamard(wires=1)
+    qp.Hadamard(wires=0)
+    m = qp.measure(wires=[0])
+    qp.cond(m == 1, qp.PauliX)(wires=1)
+    qp.Hadamard(wires=1)
 
     # Return the density matrix of the output state
-    return qml.density_matrix(wires=[1])
+    return qp.density_matrix(wires=[1])
 
 
 ##############################################################################
@@ -297,19 +297,19 @@ def generate_random_state(n=1):
 # To start off, we define the :math:`R_z(\theta)` gate using two qubits with the gate-based approach
 # so we can later compare our MBQC approach to it.
 
-dev = qml.device("lightning.qubit", wires=1)
+dev = qp.device("lightning.qubit", wires=1)
 
 
-@qml.qnode(dev, interface="autograd")
+@qp.qnode(dev, interface="autograd")
 def RZ(theta, input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=0)
+    qp.StatePrep(input_state, wires=0)
 
     # Perform the Rz rotation
-    qml.RZ(theta, wires=0)
+    qp.RZ(theta, wires=0)
 
     # Return the density matrix of the output state
-    return qml.density_matrix(wires=[0])
+    return qp.density_matrix(wires=[0])
 
 
 ##############################################################################
@@ -318,28 +318,28 @@ def RZ(theta, input_state):
 # in the MBQC formalism.
 #
 
-mbqc_dev = qml.device("lightning.qubit", wires=2)
+mbqc_dev = qp.device("lightning.qubit", wires=2)
 
 
-@qml.qnode(mbqc_dev, interface="autograd")
+@qp.qnode(mbqc_dev, interface="autograd")
 def RZ_MBQC(theta, input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=0)
+    qp.StatePrep(input_state, wires=0)
 
     # Prepare the cluster state
-    qml.Hadamard(wires=1)
-    qml.CZ(wires=[0, 1])
+    qp.Hadamard(wires=1)
+    qp.CZ(wires=[0, 1])
 
     # Measure the first qubit and correct the state
-    qml.RZ(theta, wires=0)
-    qml.Hadamard(wires=0)
-    m = qml.measure(wires=[0])
+    qp.RZ(theta, wires=0)
+    qp.Hadamard(wires=0)
+    m = qp.measure(wires=[0])
 
-    qml.cond(m == 1, qml.PauliX)(wires=1)
-    qml.Hadamard(wires=1)
+    qp.cond(m == 1, qp.PauliX)(wires=1)
+    qp.Hadamard(wires=1)
 
     # Return the density matrix of the output state
-    return qml.density_matrix(wires=[1])
+    return qp.density_matrix(wires=[1])
 
 
 ##############################################################################
@@ -359,49 +359,49 @@ def RZ_MBQC(theta, input_state):
 # For the :math:`R_x(\theta)` gate we take a similar approach.
 #
 
-dev = qml.device("lightning.qubit", wires=1)
+dev = qp.device("lightning.qubit", wires=1)
 
 
-@qml.qnode(dev, interface="autograd")
+@qp.qnode(dev, interface="autograd")
 def RX(theta, input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=0)
+    qp.StatePrep(input_state, wires=0)
 
     # Perform the Rz rotation
-    qml.RX(theta, wires=0)
+    qp.RX(theta, wires=0)
 
     # Return the density matrix of the output state
-    return qml.density_matrix(wires=[0])
+    return qp.density_matrix(wires=[0])
 
 
-mbqc_dev = qml.device("lightning.qubit", wires=3)
+mbqc_dev = qp.device("lightning.qubit", wires=3)
 
 
-@qml.qnode(mbqc_dev, interface="autograd")
+@qp.qnode(mbqc_dev, interface="autograd")
 def RX_MBQC(theta, input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=0)
+    qp.StatePrep(input_state, wires=0)
 
     # Prepare the cluster state
-    qml.Hadamard(wires=1)
-    qml.Hadamard(wires=2)
-    qml.CZ(wires=[0, 1])
-    qml.CZ(wires=[1, 2])
+    qp.Hadamard(wires=1)
+    qp.Hadamard(wires=2)
+    qp.CZ(wires=[0, 1])
+    qp.CZ(wires=[1, 2])
 
     # Measure the qubits and perform corrections
-    qml.Hadamard(wires=0)
-    m1 = qml.measure(wires=[0])
+    qp.Hadamard(wires=0)
+    m1 = qp.measure(wires=[0])
 
-    qml.RZ(theta, wires=1)
-    qml.cond(m1 == 1, qml.RX)(-2 * theta, wires=2)
-    qml.Hadamard(wires=1)
-    m2 = qml.measure(wires=[1])
+    qp.RZ(theta, wires=1)
+    qp.cond(m1 == 1, qp.RX)(-2 * theta, wires=2)
+    qp.Hadamard(wires=1)
+    m2 = qp.measure(wires=[1])
 
-    qml.cond(m2 == 1, qml.PauliX)(wires=2)
-    qml.cond(m1 == 1, qml.PauliZ)(wires=2)
+    qp.cond(m2 == 1, qp.PauliX)(wires=2)
+    qp.cond(m1 == 1, qp.PauliZ)(wires=2)
 
     # Return the density matrix of the output state
-    return qml.density_matrix(wires=[2])
+    return qp.density_matrix(wires=[2])
 
 
 ##############################################################################
@@ -439,46 +439,46 @@ def RX_MBQC(theta, input_state):
 #
 # Let's see how one can do this in PennyLane.
 
-dev = qml.device("lightning.qubit", wires=2)
+dev = qp.device("lightning.qubit", wires=2)
 
 
-@qml.qnode(dev, interface="autograd")
+@qp.qnode(dev, interface="autograd")
 def CNOT(input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=[0, 1])
-    qml.CNOT(wires=[0, 1])
+    qp.StatePrep(input_state, wires=[0, 1])
+    qp.CNOT(wires=[0, 1])
 
-    return qml.density_matrix(wires=[0, 1])
+    return qp.density_matrix(wires=[0, 1])
 
 
-mbqc_dev = qml.device("lightning.qubit", wires=4)
+mbqc_dev = qp.device("lightning.qubit", wires=4)
 
 
-@qml.qnode(mbqc_dev, interface="autograd")
+@qp.qnode(mbqc_dev, interface="autograd")
 def CNOT_MBQC(input_state):
     # Prepare the input state
-    qml.StatePrep(input_state, wires=[0, 1])
+    qp.StatePrep(input_state, wires=[0, 1])
 
     # Prepare the cluster state
-    qml.Hadamard(wires=2)
-    qml.Hadamard(wires=3)
-    qml.CZ(wires=[2, 0])
-    qml.CZ(wires=[2, 1])
-    qml.CZ(wires=[2, 3])
+    qp.Hadamard(wires=2)
+    qp.Hadamard(wires=3)
+    qp.CZ(wires=[2, 0])
+    qp.CZ(wires=[2, 1])
+    qp.CZ(wires=[2, 3])
 
     # Measure the qubits in the appropriate bases
-    qml.Hadamard(wires=1)
-    m1 = qml.measure(wires=[1])
-    qml.Hadamard(wires=2)
-    m2 = qml.measure(wires=[2])
+    qp.Hadamard(wires=1)
+    m1 = qp.measure(wires=[1])
+    qp.Hadamard(wires=2)
+    m2 = qp.measure(wires=[2])
 
     # Correct the state
-    qml.cond(m1 == 1, qml.PauliZ)(wires=0)
-    qml.cond(m2 == 1, qml.PauliX)(wires=3)
-    qml.cond(m1 == 1, qml.PauliZ)(wires=3)
+    qp.cond(m1 == 1, qp.PauliZ)(wires=0)
+    qp.cond(m2 == 1, qp.PauliX)(wires=3)
+    qp.cond(m1 == 1, qp.PauliZ)(wires=3)
 
     # Return the density matrix of the output state
-    return qml.density_matrix(wires=[0, 3])
+    return qp.density_matrix(wires=[0, 3])
 
 
 ##############################################################################
diff --git a/demonstrations_v2/tutorial_mbqc/metadata.json b/demonstrations_v2/tutorial_mbqc/metadata.json
index ff3cfdb5b3..f0d2fdc653 100644
--- a/demonstrations_v2/tutorial_mbqc/metadata.json
+++ b/demonstrations_v2/tutorial_mbqc/metadata.json
@@ -11,7 +11,7 @@
     "executable_stable": false,
     "executable_latest": false,
     "dateOfPublication": "2022-12-05T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Hardware",
         "Quantum Computing"
diff --git a/demonstrations_v2/tutorial_mcm_introduction/demo.py b/demonstrations_v2/tutorial_mcm_introduction/demo.py
index 2728396fd7..39cd544803 100644
--- a/demonstrations_v2/tutorial_mcm_introduction/demo.py
+++ b/demonstrations_v2/tutorial_mcm_introduction/demo.py
@@ -122,14 +122,14 @@
 # them in detail.
 #
 
-import pennylane as qml
+import pennylane as qp
 
-dev = qml.device("default.qubit")
+dev = qp.device("default.qubit")
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def before():
-    qml.Hadamard(0)  # Create |+> state
-    return qml.expval(qml.X(0)), qml.expval(qml.Z(0))
+    qp.Hadamard(0)  # Create |+> state
+    return qp.expval(qp.X(0)), qp.expval(qp.Z(0))
 
 b = before()
 print(f"Expectation values before any measurement: {b[0]:.1f}, {b[1]:.1f}")
@@ -164,11 +164,11 @@ def before():
 # variable to its outcome.
 #
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def after():
-    qml.Hadamard(0)  # Create |+> state
-    qml.measure(0)  # Measure without recording the outcome
-    return qml.expval(qml.X(0)), qml.expval(qml.Z(0))
+    qp.Hadamard(0)  # Create |+> state
+    qp.measure(0)  # Measure without recording the outcome
+    return qp.expval(qp.X(0)), qp.expval(qp.Z(0))
 
 a = after()
 print(f"Expectation value after the measurement:  {a[0]:.1f}, {a[1]:.1f}")
@@ -188,14 +188,14 @@ def after():
 #
 # that is, the qubit is in a new, pure state. In PennyLane, we can postselect on the case
 # where we measured a :math:`0` using the ``postselect`` keyword argument of
-# ``qml.measure``:
+# ``qp.measure``:
 #
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def after():
-    qml.Hadamard(0)  # Create |+> state
-    qml.measure(0, postselect=0)  # Measure and only accept 0 as outcome
-    return qml.expval(qml.X(0)), qml.expval(qml.Z(0))
+    qp.Hadamard(0)  # Create |+> state
+    qp.measure(0, postselect=0)  # Measure and only accept 0 as outcome
+    return qp.expval(qp.X(0)), qp.expval(qp.Z(0))
 
 a = after()
 print(f"Expectation value after the postselected measurement:  {a[0]:.1f}, {a[1]:.1f}")
@@ -224,14 +224,14 @@ def after():
 # We code this circuit up similar to the one above, using an additional ``CNOT`` gate
 # to create the Bell state. We also include optional hyperparameters such as
 # ``postselect`` as keyword arguments to our quantum function and pass them on to
-# ``qml.measure``. Note that we can't complete the quantum function yet, because we still
+# ``qp.measure``. Note that we can't complete the quantum function yet, because we still
 # need to discuss what to return from it!
 #
 
 def bell_pair_preparation(**kwargs):
-    qml.Hadamard(0)
-    qml.CNOT([0, 1])  # Create a Bell pair
-    qml.measure(0, **kwargs)  # Measure first qubit, using keyword arguments
+    qp.Hadamard(0)
+    qp.CNOT([0, 1])  # Create a Bell pair
+    qp.measure(0, **kwargs)  # Measure first qubit, using keyword arguments
 
 ######################################################################
 # Without recording the outcome, i.e., ``postselect=None``, we obtain the state
@@ -251,10 +251,10 @@ def bell_pair_preparation(**kwargs):
 # And those will be the return types to complete our quantum function:
 #
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def bell_pair(postselect):
     bell_pair_preparation(postselect=postselect)
-    return qml.purity([0, 1]), qml.vn_entropy(0)
+    return qp.purity([0, 1]), qp.vn_entropy(0)
 
 ######################################################################
 # So let's compare the purities and von Neumann entropies of the Bell state
@@ -297,10 +297,10 @@ def bell_pair(postselect):
 # can simply pass the keyword argument ``reset`` to activate the qubit reset:
 #
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def bell_pair_with_reset(reset):
     bell_pair_preparation(reset=reset)
-    return qml.expval(qml.Z(0)), qml.expval(qml.Z(1)), qml.expval(qml.Z(0) @ qml.Z(1))
+    return qp.expval(qp.Z(0)), qp.expval(qp.Z(1)), qp.expval(qp.Z(0) @ qp.Z(1))
 
 no_reset = bell_pair_with_reset(reset=False)
 reset = bell_pair_with_reset(reset=True)
@@ -360,10 +360,10 @@ def bell_pair_with_reset(reset):
 magic_state = np.array([1, np.exp(1j * np.pi / 4)]) / np.sqrt(2)
 
 def t_gadget(wire, aux_wire):
-    qml.StatePrep(magic_state, aux_wire)
-    qml.CNOT([wire, aux_wire])
-    mcm = qml.measure(aux_wire, reset=True)  # Resetting disentangles aux qubit
-    qml.cond(mcm, qml.S)(wire)  # Apply qml.S(wire) if mcm was 1
+    qp.StatePrep(magic_state, aux_wire)
+    qp.CNOT([wire, aux_wire])
+    mcm = qp.measure(aux_wire, reset=True)  # Resetting disentangles aux qubit
+    qp.cond(mcm, qp.S)(wire)  # Apply qp.S(wire) if mcm was 1
 
 ######################################################################
 # We will not derive why this works (see, e.g., [#zhou]_ instead), but
@@ -380,16 +380,16 @@ def t_gadget(wire, aux_wire):
 # - return the expectation value :math:`\langle X_0\rangle.`
 #
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def test_t_gadget(init_state):
-    qml.Hadamard(0)  # Create |+> state
+    qp.Hadamard(0)  # Create |+> state
     if init_state == "-":
-        qml.Z(0)  # Flip to |-> state
+        qp.Z(0)  # Flip to |-> state
 
     t_gadget(0, 1)  # Apply T-gadget
-    qml.adjoint(qml.T)(0)  # Apply T^† to undo the gadget
+    qp.adjoint(qp.T)(0)  # Apply T^† to undo the gadget
 
-    return qml.expval(qml.X(0))
+    return qp.expval(qp.X(0))
 
 print(f"<X₀> with initial state |+>: {test_t_gadget('+'):4.1f}")
 print(f"<X₀> with initial state |->: {test_t_gadget('-'):4.1f}")
diff --git a/demonstrations_v2/tutorial_mcm_introduction/metadata.json b/demonstrations_v2/tutorial_mcm_introduction/metadata.json
index 3baae949fe..8717812448 100644
--- a/demonstrations_v2/tutorial_mcm_introduction/metadata.json
+++ b/demonstrations_v2/tutorial_mcm_introduction/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-05-10T00:00:00+00:00",
-    "dateOfLastModification": "2026-01-15T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Getting Started",
         "Quantum Computing"
diff --git a/demonstrations_v2/tutorial_measurement_optimize/demo.py b/demonstrations_v2/tutorial_measurement_optimize/demo.py
index 83d218cc4f..a1ccf3bf71 100644
--- a/demonstrations_v2/tutorial_measurement_optimize/demo.py
+++ b/demonstrations_v2/tutorial_measurement_optimize/demo.py
@@ -110,11 +110,11 @@
 import warnings
 import jax
 from jax import numpy as jnp
-import pennylane as qml
+import pennylane as qp
 
 jax.config.update("jax_enable_x64", True)
 
-dataset = qml.data.load("qchem", molname="H2", bondlength=0.7)[0]
+dataset = qp.data.load("qchem", molname="H2", bondlength=0.7)[0]
 H, num_qubits = dataset.hamiltonian, len(dataset.hamiltonian.wires)
 
 print("Required number of qubits:", num_qubits)
@@ -125,27 +125,27 @@
 # on hardware. Let's generate the cost function to check this.
 
 # Create a 4 qubit simulator
-dev = qml.device("default.qubit", seed=904932)
+dev = qp.device("default.qubit", seed=904932)
 
 # number of electrons
 electrons = 2
 
 # Define the Hartree-Fock initial state for our variational circuit
-initial_state = qml.qchem.hf_state(electrons, num_qubits)
+initial_state = qp.qchem.hf_state(electrons, num_qubits)
 
 # Construct the UCCSD ansatz
-singles, doubles = qml.qchem.excitations(electrons, num_qubits)
-s_wires, d_wires = qml.qchem.excitations_to_wires(singles, doubles)
+singles, doubles = qp.qchem.excitations(electrons, num_qubits)
+s_wires, d_wires = qp.qchem.excitations_to_wires(singles, doubles)
 ansatz = functools.partial(
-    qml.UCCSD, init_state=initial_state, s_wires=s_wires, d_wires=d_wires
+    qp.UCCSD, init_state=initial_state, s_wires=s_wires, d_wires=d_wires
 )
 
 # generate the cost function
-@qml.set_shots(1000)
-@qml.qnode(dev, interface="jax")
+@qp.set_shots(1000)
+@qp.qnode(dev, interface="jax")
 def cost_circuit(params):
     ansatz(params, wires=range(num_qubits))
-    return qml.expval(H)
+    return qp.expval(H)
 
 ##############################################################################
 # If we evaluate this cost function, we can see that it corresponds to 15 different
@@ -154,7 +154,7 @@ def cost_circuit(params):
 from jax import random as random
 key, scale = random.PRNGKey(0), jnp.pi
 params = random.normal(key, shape=(len(singles) + len(doubles),)) * scale
-with qml.Tracker(dev) as tracker:  # track the number of executions
+with qp.Tracker(dev) as tracker:  # track the number of executions
     print("Cost function value:", cost_circuit(params))
 
 print("Number of quantum evaluations:", tracker.totals['executions'])
@@ -163,7 +163,7 @@ def cost_circuit(params):
 # How about a larger molecule? Let's try the
 # `water molecule <https://pennylane.ai/datasets/qchem/h2o-molecule>`__:
 
-dataset = qml.data.load('qchem', molname="H2O")[0]
+dataset = qp.data.load('qchem', molname="H2O")[0]
 H, num_qubits = dataset.hamiltonian, len(dataset.hamiltonian.wires)
 
 print("Required number of qubits:", num_qubits)
@@ -359,8 +359,8 @@ def cost_circuit(params):
 
 
 obs = [
-    qml.PauliX(0) @ qml.PauliY(1),
-    qml.PauliX(0) @ qml.PauliZ(2)
+    qp.PauliX(0) @ qp.PauliY(1),
+    qp.PauliX(0) @ qp.PauliZ(2)
 ]
 
 
@@ -369,20 +369,20 @@ def cost_circuit(params):
 # the two QWC terms.
 
 
-dev = qml.device("default.qubit", wires=3)
+dev = qp.device("default.qubit", wires=3)
 
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def circuit1(weights):
-    qml.StronglyEntanglingLayers(weights, wires=range(3))
-    return qml.expval(obs[0])
+    qp.StronglyEntanglingLayers(weights, wires=range(3))
+    return qp.expval(obs[0])
 
 
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def circuit2(weights):
-    qml.StronglyEntanglingLayers(weights, wires=range(3))
-    return qml.expval(obs[1])
+    qp.StronglyEntanglingLayers(weights, wires=range(3))
+    return qp.expval(obs[1])
 
-param_shape = qml.templates.StronglyEntanglingLayers.shape(n_layers=3, n_wires=3)
+param_shape = qp.templates.StronglyEntanglingLayers.shape(n_layers=3, n_wires=3)
 key, scale = random.PRNGKey(192933), 0.1
 weights = scale * random.normal(key, shape=param_shape)
 
@@ -393,20 +393,20 @@ def circuit2(weights):
 # Now, let's use our QWC approach to reduce this down to a *single* measurement
 # of the probabilities in the shared eigenbasis of both QWC observables:
 
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def circuit_qwc(weights):
-    qml.StronglyEntanglingLayers(weights, wires=range(3))
+    qp.StronglyEntanglingLayers(weights, wires=range(3))
 
     # rotate wire 0 into the shared eigenbasis
-    qml.RY(-jnp.pi / 2, wires=0)
+    qp.RY(-jnp.pi / 2, wires=0)
 
     # rotate wire 1 into the shared eigenbasis
-    qml.RX(jnp.pi / 2, wires=1)
+    qp.RX(jnp.pi / 2, wires=1)
 
     # wire 2 does not require a rotation
 
     # measure probabilities in the computational basis
-    return qml.probs(wires=range(3))
+    return qp.probs(wires=range(3))
 
 
 rotated_probs = circuit_qwc(weights)
@@ -438,12 +438,12 @@ def circuit_qwc(weights):
 # Luckily, PennyLane automatically performs this QWC grouping under the hood. We simply
 # return the two QWC Pauli terms from the QNode:
 
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def circuit(weights):
-    qml.StronglyEntanglingLayers(weights, wires=range(3))
+    qp.StronglyEntanglingLayers(weights, wires=range(3))
     return [
-        qml.expval(qml.PauliX(0) @ qml.PauliY(1)),
-        qml.expval(qml.PauliX(0) @ qml.PauliZ(2))
+        qp.expval(qp.PauliX(0) @ qp.PauliY(1)),
+        qp.expval(qp.PauliX(0) @ qp.PauliZ(2))
     ]
 
 
@@ -452,10 +452,10 @@ def circuit(weights):
 
 ##############################################################################
 # Behind the scenes, PennyLane is making use of our built-in
-# :mod:`qml.pauli <pennylane.pauli>` module, which contains functions for diagonalizing QWC
+# :mod:`qp.pauli <pennylane.pauli>` module, which contains functions for diagonalizing QWC
 # terms:
 
-rotations, new_obs = qml.pauli.diagonalize_qwc_pauli_words(obs)
+rotations, new_obs = qp.pauli.diagonalize_qwc_pauli_words(obs)
 
 print(rotations)
 print(new_obs)
@@ -463,7 +463,7 @@ def circuit(weights):
 
 ##############################################################################
 # Here, the first line corresponds to the basis rotations that were discussed above, written in
-# terms of ``RX`` and ``RY`` rotations. Check out the :mod:`qml.pauli <pennylane.pauli>`
+# terms of ``RX`` and ``RY`` rotations. Check out the :mod:`qp.pauli <pennylane.pauli>`
 # documentation for more details on its provided functionality and how it works.
 #
 # Given a Hamiltonian containing a large number of Pauli terms,
@@ -543,19 +543,19 @@ def circuit(weights):
 from matplotlib import pyplot as plt
 
 terms = [
-    qml.PauliZ(0),
-    qml.PauliZ(0) @ qml.PauliZ(1),
-    qml.PauliZ(0) @ qml.PauliZ(1) @ qml.PauliZ(2),
-    qml.PauliZ(0) @ qml.PauliZ(1) @ qml.PauliZ(2) @ qml.PauliZ(3),
-    qml.PauliX(2) @ qml.PauliX(3),
-    qml.PauliY(0) @ qml.PauliX(2) @ qml.PauliX(3),
-    qml.PauliY(0) @ qml.PauliY(1) @ qml.PauliX(2) @ qml.PauliX(3)
+    qp.PauliZ(0),
+    qp.PauliZ(0) @ qp.PauliZ(1),
+    qp.PauliZ(0) @ qp.PauliZ(1) @ qp.PauliZ(2),
+    qp.PauliZ(0) @ qp.PauliZ(1) @ qp.PauliZ(2) @ qp.PauliZ(3),
+    qp.PauliX(2) @ qp.PauliX(3),
+    qp.PauliY(0) @ qp.PauliX(2) @ qp.PauliX(3),
+    qp.PauliY(0) @ qp.PauliY(1) @ qp.PauliX(2) @ qp.PauliX(3)
 ]
 
 def format_pauli_word(term):
     """Convenience function that nicely formats a PennyLane
     tensor observable as a Pauli word"""
-    if isinstance(term, qml.ops.Prod):
+    if isinstance(term, qp.ops.Prod):
         return " ".join([format_pauli_word(t) for t in term])
 
     return f"{term.name[-1]}{term.wires.tolist()[0]}"
@@ -675,9 +675,9 @@ def format_pauli_word(term):
 # the entire process using the provided grouping functions.
 #
 # Steps 1-3 (finding and grouping QWC terms in the Hamiltonian) can be done via the
-# :func:`qml.pauli.group_observables <pennylane.pauli.group_observables>` function:
+# :func:`qp.pauli.group_observables <pennylane.pauli.group_observables>` function:
 
-obs_groupings = qml.pauli.group_observables(terms, grouping_type='qwc', method='rlf')
+obs_groupings = qp.pauli.group_observables(terms, grouping_type='qwc', method='rlf')
 
 
 ##############################################################################
@@ -686,23 +686,23 @@ def format_pauli_word(term):
 # heuristic (in this case, ``"rlf"`` refers to Recursive Largest First, a variant of largest first colouring heuristic).
 #
 # If we want to see what the required rotations and measurements are, we can use the
-# :func:`qml.pauli.diagonalize_qwc_groupings <pennylane.pauli.diagonalize_qwc_groupings>`
+# :func:`qp.pauli.diagonalize_qwc_groupings <pennylane.pauli.diagonalize_qwc_groupings>`
 # function:
 
-rotations, measurements = qml.pauli.diagonalize_qwc_groupings(obs_groupings)
+rotations, measurements = qp.pauli.diagonalize_qwc_groupings(obs_groupings)
 
 ##############################################################################
 # However, this isn't strictly necessary—recall previously that the QNode
 # has the capability to *automatically* measure qubit-wise commuting observables!
 
-dev = qml.device("lightning.qubit", wires=4)
+dev = qp.device("lightning.qubit", wires=4)
 
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def circuit(weights, group=None, **kwargs):
-    qml.StronglyEntanglingLayers(weights, wires=range(4))
-    return [qml.expval(o) for o in group]
+    qp.StronglyEntanglingLayers(weights, wires=range(4))
+    return [qp.expval(o) for o in group]
 
-param_shape = qml.templates.StronglyEntanglingLayers.shape(n_layers=3, n_wires=4)
+param_shape = qp.templates.StronglyEntanglingLayers.shape(n_layers=3, n_wires=4)
 key = random.PRNGKey(1)
 weights = random.normal(key, shape=param_shape) * 0.1
 result = [jnp.array(circuit(weights, group=g)) for g in obs_groupings]
@@ -722,12 +722,12 @@ def circuit(weights, group=None, **kwargs):
 # problems), we can use the option ``grouping_type="qwc"`` in :class:`~.pennylane.Hamiltonian` to
 # automatically optimize the measurements.
 
-H = qml.Hamiltonian(coeffs=jnp.ones(len(terms)), observables=terms, grouping_type="qwc")
+H = qp.Hamiltonian(coeffs=jnp.ones(len(terms)), observables=terms, grouping_type="qwc")
 _, H_ops = H.terms()
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def cost_fn(weights):
-    qml.StronglyEntanglingLayers(weights, wires=range(4))
-    return qml.expval(H)
+    qp.StronglyEntanglingLayers(weights, wires=range(4))
+    return qp.expval(H)
 print(cost_fn(weights))
 
 ##############################################################################
@@ -738,12 +738,12 @@ def cost_fn(weights):
 # how this affects the number of measurements required to perform the VQE on :math:`\text{H}_2 \text{O}!`
 # Let's use our new-found knowledge to see what happens.
 
-dataset = qml.data.load('qchem', molname="H2O")[0]
+dataset = qp.data.load('qchem', molname="H2O")[0]
 H, num_qubits = dataset.hamiltonian, len(dataset.hamiltonian.wires)
 print("Number of Hamiltonian terms/required measurements:", len(H_ops))
 
 # grouping
-groups = qml.pauli.group_observables(H_ops, grouping_type='qwc', method='rlf')
+groups = qp.pauli.group_observables(H_ops, grouping_type='qwc', method='rlf')
 print("Number of required measurements after optimization:", len(groups))
 
 ##############################################################################
diff --git a/demonstrations_v2/tutorial_measurement_optimize/metadata.json b/demonstrations_v2/tutorial_measurement_optimize/metadata.json
index 3c1d70ff0b..c03f02dea2 100644
--- a/demonstrations_v2/tutorial_measurement_optimize/metadata.json
+++ b/demonstrations_v2/tutorial_measurement_optimize/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2021-01-18T00:00:00+00:00",
-    "dateOfLastModification": "2025-11-10T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Quantum Chemistry"
     ],
diff --git a/demonstrations_v2/tutorial_mitigation_advantage/demo.py b/demonstrations_v2/tutorial_mitigation_advantage/demo.py
index 037d4fa863..a2291e225b 100644
--- a/demonstrations_v2/tutorial_mitigation_advantage/demo.py
+++ b/demonstrations_v2/tutorial_mitigation_advantage/demo.py
@@ -79,7 +79,7 @@
 at the classical mixtures introduced by the Kraus operators of the noise channel. That is why we need to use the mixed state simulator.
 For more information see e.g. our :doc:`demo on simulating noisy circuits <demos/tutorial_noisy_circuits>`.
 """
-import pennylane as qml
+import pennylane as qp
 import jax
 import jax.numpy as jnp
 import matplotlib.pyplot as plt
@@ -88,12 +88,12 @@
 
 n_wires = 9
 # Describe noise
-noise_gate = qml.DepolarizingChannel
+noise_gate = qp.DepolarizingChannel
 p = 0.005
 
 # Load devices
-dev_ideal = qml.device("default.mixed", wires=n_wires)
-dev_noisy = qml.noise.insert(dev_ideal, noise_gate, p, position="all")
+dev_ideal = qp.device("default.mixed", wires=n_wires)
+dev_noisy = qp.noise.insert(dev_ideal, noise_gate, p, position="all")
 
 # 3x3 grid with nearest neighbors
 connections = [(0, 1), (1, 2),
@@ -103,16 +103,16 @@
                (1, 4), (4, 7),
                (2, 5), (5, 8)]
 
-def time_evolution(theta_h, n_layers = 10, obs = qml.PauliZ(4)):
+def time_evolution(theta_h, n_layers = 10, obs = qp.PauliZ(4)):
     for _ in range(n_layers):
         for i, j in connections:
-            qml.IsingZZ(-jnp.pi/2, wires=(i, j))
+            qp.IsingZZ(-jnp.pi/2, wires=(i, j))
         for i in range(n_wires):
-            qml.RX(theta_h, wires=i)
-    return qml.expval(obs)
+            qp.RX(theta_h, wires=i)
+    return qp.expval(obs)
 
-qnode_ideal = qml.QNode(time_evolution, dev_ideal, interface="jax")
-qnode_noisy = qml.QNode(time_evolution, dev_noisy, interface="jax")
+qnode_ideal = qp.QNode(time_evolution, dev_ideal, interface="jax")
+qnode_noisy = qp.QNode(time_evolution, dev_noisy, interface="jax")
 
 ##############################################################################
 # We can now simulate the final expectation value with and without noise. Note that the ``IsingZZ`` gate is not natively
@@ -167,11 +167,11 @@ def time_evolution(theta_h, n_layers = 10, obs = qml.PauliZ(4)):
 # our model by an appropriate gain factor. Here, :math:`G=(1, 1.2, 1.6)` in accordance with [#ibm]_. In order to do this in PennyLane, we simply
 # set up two new noisy devices with the appropriately attenuated noise parameters.
 
-dev_noisy1 = qml.noise.insert(dev_ideal, noise_gate, p*1.2, position="all")
-dev_noisy2 = qml.noise.insert(dev_ideal, noise_gate, p*1.6, position="all")
+dev_noisy1 = qp.noise.insert(dev_ideal, noise_gate, p*1.2, position="all")
+dev_noisy2 = qp.noise.insert(dev_ideal, noise_gate, p*1.6, position="all")
 
-qnode_noisy1 = qml.QNode(time_evolution, dev_noisy1, interface="jax")
-qnode_noisy2 = qml.QNode(time_evolution, dev_noisy2, interface="jax")
+qnode_noisy1 = qp.QNode(time_evolution, dev_noisy1, interface="jax")
+qnode_noisy2 = qp.QNode(time_evolution, dev_noisy2, interface="jax")
 
 res_noisy1 = jax.vmap(qnode_noisy1)(thetas)
 res_noisy2 = jax.vmap(qnode_noisy2)(thetas)
@@ -197,7 +197,7 @@ def time_evolution(theta_h, n_layers = 10, obs = qml.PauliZ(4)):
 # We now repeat this procedure for all values of :math:`\theta_h` and see how the results are much improved.
 # We can use :func:`~pennylane.noise.richardson_extrapolate` that performs a polynomial fit of a degree matching the input data size.
 
-res_mitigated = [qml.noise.richardson_extrapolate(Gs, [res_noisy[i], res_noisy1[i], res_noisy2[i]]) for i in range(len(res_ideal))]
+res_mitigated = [qp.noise.richardson_extrapolate(Gs, [res_noisy[i], res_noisy1[i], res_noisy2[i]]) for i in range(len(res_ideal))]
 
 plt.plot(thetas, res_ideal, label="exact")
 plt.plot(thetas, res_mitigated, label="mitigated")
diff --git a/demonstrations_v2/tutorial_mitigation_advantage/metadata.json b/demonstrations_v2/tutorial_mitigation_advantage/metadata.json
index 69626781c8..328335b6e3 100644
--- a/demonstrations_v2/tutorial_mitigation_advantage/metadata.json
+++ b/demonstrations_v2/tutorial_mitigation_advantage/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2023-06-16T00:00:00+00:00",
-    "dateOfLastModification": "2025-10-15T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Quantum Hardware",
         "Quantum Computing"
diff --git a/demonstrations_v2/tutorial_mol_geo_opt/demo.py b/demonstrations_v2/tutorial_mol_geo_opt/demo.py
index 91d6315eef..a8fc8c6092 100644
--- a/demonstrations_v2/tutorial_mol_geo_opt/demo.py
+++ b/demonstrations_v2/tutorial_mol_geo_opt/demo.py
@@ -123,12 +123,12 @@
 # of the trihydrogen cation using the
 # :func:`~.pennylane.qchem.molecular_hamiltonian` function.
 
-import pennylane as qml
+import pennylane as qp
 
 
 def H(x):
-    molecule = qml.qchem.Molecule(symbols, x, charge=1)
-    return qml.qchem.molecular_hamiltonian(molecule)[0]
+    molecule = qp.qchem.Molecule(symbols, x, charge=1)
+    return qp.qchem.molecular_hamiltonian(molecule)[0]
 
 
 ##############################################################################
@@ -185,7 +185,7 @@ def H(x):
 # :func:`~.pennylane.qchem.hf_state` function to generate the
 # occupation-number vector representing the Hartree-Fock state
 
-hf = qml.qchem.hf_state(electrons=2, orbitals=6)
+hf = qp.qchem.hf_state(electrons=2, orbitals=6)
 print(hf)
 
 ##############################################################################
@@ -194,16 +194,16 @@ def H(x):
 # First, we define the quantum device used to compute the expectation value.
 # In this example, we use the ``default.qubit`` simulator:
 num_wires = 6
-dev = qml.device("default.qubit", wires=num_wires)
+dev = qp.device("default.qubit", wires=num_wires)
 
 
-@qml.qnode(dev, interface="jax")
+@qp.qnode(dev, interface="jax")
 def circuit(params, obs, wires):
-    qml.BasisState(hf, wires=wires)
-    qml.DoubleExcitation(params[0], wires=[0, 1, 2, 3])
-    qml.DoubleExcitation(params[1], wires=[0, 1, 4, 5])
+    qp.BasisState(hf, wires=wires)
+    qp.DoubleExcitation(params[0], wires=[0, 1, 2, 3])
+    qp.DoubleExcitation(params[1], wires=[0, 1, 4, 5])
 
-    return qml.expval(obs)
+    return qp.expval(obs)
 
 
 ##############################################################################
diff --git a/demonstrations_v2/tutorial_mol_geo_opt/metadata.json b/demonstrations_v2/tutorial_mol_geo_opt/metadata.json
index 5d9c500c9b..bd1ab0a1bc 100644
--- a/demonstrations_v2/tutorial_mol_geo_opt/metadata.json
+++ b/demonstrations_v2/tutorial_mol_geo_opt/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2021-06-30T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Chemistry"
     ],
diff --git a/demonstrations_v2/tutorial_mps/demo.py b/demonstrations_v2/tutorial_mps/demo.py
index 444eb70734..f1eb037c36 100644
--- a/demonstrations_v2/tutorial_mps/demo.py
+++ b/demonstrations_v2/tutorial_mps/demo.py
@@ -571,24 +571,24 @@ def dense_to_mps(psi, bond_dim):
 # Let us run a VQE example from using the `PennyLane Datasets <https://pennylane.ai/datasets/>`__ data for the :math:`H_6` molecule. What you will see here is mostly boilerplate code in PennyLane for using ``default.tensor`,` and you can
 # see our :doc:`demo on how to use default.tensor <demos/tutorial_How_to_simulate_quantum_circuits_with_tensor_networks/>` for more details.
 
-import pennylane as qml
+import pennylane as qp
 
-[dataset] = qml.data.load("qchem", molname="H6", bondlength=1.3, basis="STO-3G")
+[dataset] = qp.data.load("qchem", molname="H6", bondlength=1.3, basis="STO-3G")
 
 H = dataset.hamiltonian # molecular Hamiltonian in qubit basis
 n_wires = len(H.wires)  # number of qubits
 
 def circuit():
-    qml.BasisState(dataset.hf_state, wires=H.wires) # Hartree–Fock initial state
+    qp.BasisState(dataset.hf_state, wires=H.wires) # Hartree–Fock initial state
     for op in dataset.vqe_gates:                    # Applying all pre-optimized VQE gates
-        qml.apply(op)               
-    return qml.expval(H)                            # expectation value of molecular Hamiltonian
+        qp.apply(op)               
+    return qp.expval(H)                            # expectation value of molecular Hamiltonian
 
 # set up device with hyper-parameters and kwargs
-mps = qml.device("default.tensor", wires=n_wires, method="mps", max_bond_dim=30, contract="auto-mps")
+mps = qp.device("default.tensor", wires=n_wires, method="mps", max_bond_dim=30, contract="auto-mps")
 
 # Create the QNode to execute the circuit on the device, and call it (w/o arguments)
-res = qml.QNode(circuit, mps)()
+res = qp.QNode(circuit, mps)()
 
 # Compare MPS simulation result with pre-optimized state-vector result
 res, dataset.vqe_energy
@@ -611,8 +611,8 @@ def circuit():
 ress = []
 
 for bond_dim in bond_dims:
-    mps = qml.device("default.tensor", wires=n_wires, method="mps", max_bond_dim=bond_dim, contract="auto-mps")
-    res = qml.QNode(circuit, mps)()
+    mps = qp.device("default.tensor", wires=n_wires, method="mps", max_bond_dim=bond_dim, contract="auto-mps")
+    res = qp.QNode(circuit, mps)()
     ress.append(res)
 
 
diff --git a/demonstrations_v2/tutorial_mps/metadata.json b/demonstrations_v2/tutorial_mps/metadata.json
index 5302ecb0dd..1a82840021 100644
--- a/demonstrations_v2/tutorial_mps/metadata.json
+++ b/demonstrations_v2/tutorial_mps/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2024-09-29T00:00:00+00:00",
-    "dateOfLastModification": "2026-01-14T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing",
         "Algorithms",
diff --git a/demonstrations_v2/tutorial_multiclass_classification/demo.py b/demonstrations_v2/tutorial_multiclass_classification/demo.py
index 46789a958b..a6669bf3f2 100644
--- a/demonstrations_v2/tutorial_multiclass_classification/demo.py
+++ b/demonstrations_v2/tutorial_multiclass_classification/demo.py
@@ -47,7 +47,7 @@
 of size 4.
 """
 
-import pennylane as qml
+import pennylane as qp
 import torch
 import numpy as np
 from torch.autograd import Variable
@@ -67,7 +67,7 @@
 num_layers = 6
 total_iterations = 100
 
-dev = qml.device("default.qubit", wires=num_qubits)
+dev = qp.device("default.qubit", wires=num_qubits)
 
 
 #################################################################################
@@ -81,12 +81,12 @@
 
 def layer(W):
     for i in range(num_qubits):
-        qml.Rot(W[i, 0], W[i, 1], W[i, 2], wires=i)
+        qp.Rot(W[i, 0], W[i, 1], W[i, 2], wires=i)
     for j in range(num_qubits - 1):
-        qml.CNOT(wires=[j, j + 1])
+        qp.CNOT(wires=[j, j + 1])
     if num_qubits >= 2:
         # Apply additional CNOT to entangle the last with the first qubit
-        qml.CNOT(wires=[num_qubits - 1, 0])
+        qp.CNOT(wires=[num_qubits - 1, 0])
 
 
 #################################################################################
@@ -104,17 +104,17 @@ def layer(W):
 
 
 def circuit(weights, feat=None):
-    qml.AmplitudeEmbedding(feat, range(num_qubits), pad_with=0.0, normalize=True)
+    qp.AmplitudeEmbedding(feat, range(num_qubits), pad_with=0.0, normalize=True)
 
     for W in weights:
         layer(W)
 
-    return qml.expval(qml.PauliZ(0))
+    return qp.expval(qp.PauliZ(0))
 
 
 qnodes = []
 for iq in range(num_classes):
-    qnode = qml.QNode(circuit, dev, interface="torch")
+    qnode = qp.QNode(circuit, dev, interface="torch")
     qnodes.append(qnode)
 
 
diff --git a/demonstrations_v2/tutorial_multiclass_classification/metadata.json b/demonstrations_v2/tutorial_multiclass_classification/metadata.json
index aadd811853..4fae5a0130 100644
--- a/demonstrations_v2/tutorial_multiclass_classification/metadata.json
+++ b/demonstrations_v2/tutorial_multiclass_classification/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2020-04-09T00:00:00+00:00",
-    "dateOfLastModification": "2025-12-19T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Algorithms"
     ],
diff --git a/demonstrations_v2/tutorial_neutral_atoms/demo.py b/demonstrations_v2/tutorial_neutral_atoms/demo.py
index 478dfaaf78..cac9ead1d6 100644
--- a/demonstrations_v2/tutorial_neutral_atoms/demo.py
+++ b/demonstrations_v2/tutorial_neutral_atoms/demo.py
@@ -304,7 +304,7 @@
 # this can be easily changed in programmable devices.
 # Let's plot this function to get an idea of what the pulse looks like. First, let's import all the relevant libraries.
 
-import pennylane as qml
+import pennylane as qp
 from pennylane import numpy as np
 import matplotlib.pyplot as plt
 import jax
@@ -359,16 +359,16 @@ def blackman_window(peak, time):
 detuning = 0
 
 # For now, let's act on only one neutral atom
-single_qubit_dev = qml.device("default.qubit", wires=1)
+single_qubit_dev = qp.device("default.qubit", wires=1)
 
 
-@qml.qnode(single_qubit_dev)
+@qp.qnode(single_qubit_dev)
 def state_evolution():
 
-    # Use qml.evolve to find the final state of the atom after interacting with the pulse
-    qml.evolve(H_d(blackman_window, phase, detuning, wires=[0]))([peak], t=[0, duration])
+    # Use qp.evolve to find the final state of the atom after interacting with the pulse
+    qp.evolve(H_d(blackman_window, phase, detuning, wires=[0]))([peak], t=[0, duration])
 
-    return qml.state()
+    return qp.state()
 
 
 print("The final state is {}".format(state_evolution().round(2)))
@@ -384,13 +384,13 @@ def state_evolution():
 detuning = 100  # Some large detuning to prove the point
 
 
-@qml.qnode(single_qubit_dev)
+@qp.qnode(single_qubit_dev)
 def state_evolution_detuned():
 
-    # Use qml.evolve to find the final state of the atom after interacting with the pulse
-    qml.evolve(H_d(blackman_window, phase, detuning, wires=[0]))([peak], t=[0, duration])
+    # Use qp.evolve to find the final state of the atom after interacting with the pulse
+    qp.evolve(H_d(blackman_window, phase, detuning, wires=[0]))([peak], t=[0, duration])
 
-    return qml.state()
+    return qp.state()
 
 
 print(
@@ -438,17 +438,17 @@ def neutral_atom_RX(theta):
     peak = theta / duration / 0.42 / (2 * jnp.pi)  # Recall that duration is 0.2
 
     # Set phase and detuning equal to zero for RX gate
-    qml.evolve(H_d(blackman_window, 0, 0, wires=[0]))([peak], t=[0, duration])
+    qp.evolve(H_d(blackman_window, 0, 0, wires=[0]))([peak], t=[0, duration])
 
 
 print(
     "For theta = pi/2, the matrix for the pulse-based RX gate is \n {} \n".format(
-        qml.matrix(neutral_atom_RX, wire_order=[0])(jnp.pi / 2).round(2)
+        qp.matrix(neutral_atom_RX, wire_order=[0])(jnp.pi / 2).round(2)
     )
 )
 print(
     "The matrix for the exact RX(pi/2) gate is \n {}".format(
-        qml.matrix(qml.RX(jnp.pi / 2, wires=0)).round(2)
+        qp.matrix(qp.RX(jnp.pi / 2, wires=0)).round(2)
     )
 )
 ##############################################################################
@@ -461,17 +461,17 @@ def neutral_atom_RY(theta):
     peak = theta / duration / 0.42 / (2 * jnp.pi)  # Recall that duration is 0.2
 
     # Set phase equal to pi/2 and detuning equal to zero for RY gate
-    qml.evolve(H_d(blackman_window, -jnp.pi / 2, 0, wires=[0]))([peak], t=[0, duration])
+    qp.evolve(H_d(blackman_window, -jnp.pi / 2, 0, wires=[0]))([peak], t=[0, duration])
 
 
 print(
     "For theta = pi/2, the matrix for the pulse-based RY gate is \n {} \n".format(
-        qml.matrix(neutral_atom_RY, wire_order=[0])(jnp.pi / 2).round(2)
+        qp.matrix(neutral_atom_RY, wire_order=[0])(jnp.pi / 2).round(2)
     )
 )
 print(
     "The matrix for the exact RY(pi/2) gate is \n {}".format(
-        qml.matrix(qml.RY(jnp.pi / 2, wires=0)).round(2)
+        qp.matrix(qp.RY(jnp.pi / 2, wires=0)).round(2)
     )
 )
 ##############################################################################
@@ -521,7 +521,7 @@ def H_i(distance, coupling):
     atomic_coordinates = [[0, 0], [0, distance]]
 
     # Return the interaction term for two atoms in terms of the distance
-    return qml.pulse.rydberg_interaction(
+    return qp.pulse.rydberg_interaction(
         atomic_coordinates, interaction_coeff=coupling, wires=[0, 1]
     )
 
@@ -551,7 +551,7 @@ def energy_gap(distance):
     H = (interaction_term + drive_term)([peak], time)
 
     # Calculate the eigenvalues for the full Hamiltonian
-    eigenvalues = jnp.linalg.eigvals(qml.matrix(H))
+    eigenvalues = jnp.linalg.eigvals(qp.matrix(H))
     return jnp.sort(eigenvalues - eigenvalues[0])
 
 
@@ -608,7 +608,7 @@ def two_pi_pulse(distance, coupling, wires=[0]):
     full_hamiltonian = H_d(blackman_window, 0, 0, wires) + H_i(distance, coupling)
 
     # Return the 2 pi pulse
-    qml.evolve(full_hamiltonian)([2 * jnp.pi / 0.42 / 0.2 / (2 * jnp.pi)], t=[0, 0.2])
+    qp.evolve(full_hamiltonian)([2 * jnp.pi / 0.42 / 0.2 / (2 * jnp.pi)], t=[0, 0.2])
 
 
 def pi_pulse(distance, coupling, wires=[0]):
@@ -616,7 +616,7 @@ def pi_pulse(distance, coupling, wires=[0]):
     full_hamiltonian = H_d(blackman_window, 0, 0, wires) + H_i(distance, coupling)
 
     # Return the pi pulse
-    qml.evolve(full_hamiltonian)([jnp.pi / 0.42 / 0.2 / (2 * jnp.pi)], t=[0, 0.2])
+    qp.evolve(full_hamiltonian)([jnp.pi / 0.42 / 0.2 / (2 * jnp.pi)], t=[0, 0.2])
 
 
 ##############################################################################
@@ -626,10 +626,10 @@ def pi_pulse(distance, coupling, wires=[0]):
 # Then, let's see the effect the sequence of pulses has on the :math:`\vert 00 \rangle` state when the atoms are close enough.
 #
 
-dev_two_qubits = qml.device("default.qubit", wires=2)
+dev_two_qubits = qp.device("default.qubit", wires=2)
 
 
-@qml.qnode(dev_two_qubits)
+@qp.qnode(dev_two_qubits)
 def neutral_atom_CZ(distance, coupling):
 
     pi_pulse(distance, coupling, wires=[0])
@@ -638,7 +638,7 @@ def neutral_atom_CZ(distance, coupling):
 
     pi_pulse(distance, coupling, wires=[0])
 
-    return qml.state()
+    return qp.state()
 
 
 print(
diff --git a/demonstrations_v2/tutorial_neutral_atoms/metadata.json b/demonstrations_v2/tutorial_neutral_atoms/metadata.json
index ffefd50101..09d7877354 100644
--- a/demonstrations_v2/tutorial_neutral_atoms/metadata.json
+++ b/demonstrations_v2/tutorial_neutral_atoms/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2023-05-30T00:00:00+00:00",
-    "dateOfLastModification": "2025-12-19T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Quantum Hardware",
         "Quantum Computing"
diff --git a/demonstrations_v2/tutorial_noisy_circuit_optimization/demo.py b/demonstrations_v2/tutorial_noisy_circuit_optimization/demo.py
index 434e3537ad..ccebd61814 100644
--- a/demonstrations_v2/tutorial_noisy_circuit_optimization/demo.py
+++ b/demonstrations_v2/tutorial_noisy_circuit_optimization/demo.py
@@ -79,53 +79,53 @@
 to carry out our noisy simulations.
 """
 
-import pennylane as qml
+import pennylane as qp
 from pennylane import numpy as np
 import matplotlib.pyplot as plt
 
-dev = qml.device("cirq.mixedsimulator", wires=2)
+dev = qp.device("cirq.mixedsimulator", wires=2)
 
 # CHSH observables
-A1 = qml.PauliZ(0)
-A2 = qml.PauliX(0)
-B1 = qml.Hermitian(np.array([[1, 1], [1, -1]]) / np.sqrt(2), wires=1)
-B2 = qml.Hermitian(np.array([[1, -1], [-1, -1]]) / np.sqrt(2), wires=1)
+A1 = qp.PauliZ(0)
+A2 = qp.PauliX(0)
+B1 = qp.Hermitian(np.array([[1, 1], [1, -1]]) / np.sqrt(2), wires=1)
+B2 = qp.Hermitian(np.array([[1, -1], [-1, -1]]) / np.sqrt(2), wires=1)
 CHSH_observables = [A1 @ B1, A1 @ B2, A2 @ B1, A2 @ B2]
 
 
 # subcircuit for creating an entangled pair of qubits
 def bell_pair():
-    qml.Hadamard(wires=0)
-    qml.CNOT(wires=[0, 1])
+    qp.Hadamard(wires=0)
+    qp.CNOT(wires=[0, 1])
 
 
 # circuits for measuring each distinct observable
-@qml.set_shots(1000)
-@qml.qnode(dev)
+@qp.set_shots(1000)
+@qp.qnode(dev)
 def measure_A1B1():
     bell_pair()
-    return qml.expval(A1 @ B1)
+    return qp.expval(A1 @ B1)
 
 
-@qml.set_shots(1000)
-@qml.qnode(dev)
+@qp.set_shots(1000)
+@qp.qnode(dev)
 def measure_A1B2():
     bell_pair()
-    return qml.expval(A1 @ B2)
+    return qp.expval(A1 @ B2)
 
 
-@qml.set_shots(1000)
-@qml.qnode(dev)
+@qp.set_shots(1000)
+@qp.qnode(dev)
 def measure_A2B1():
     bell_pair()
-    return qml.expval(A2 @ B1)
+    return qp.expval(A2 @ B1)
 
 
-@qml.set_shots(1000)
-@qml.qnode(dev)
+@qp.set_shots(1000)
+@qp.qnode(dev)
 def measure_A2B2():
     bell_pair()
-    return qml.expval(A2 @ B2)
+    return qp.expval(A2 @ B2)
 
 
 # now we measure each circuit and construct the CHSH inequality
@@ -181,8 +181,8 @@ def measure_A2B2():
     # we overwrite the bell_pair() subcircuit to add
     # extra noisy channels after the entangled state is created
     def bell_pair():
-        qml.Hadamard(wires=0)
-        qml.CNOT(wires=[0, 1])
+        qp.Hadamard(wires=0)
+        qp.CNOT(wires=[0, 1])
         cirq_ops.BitFlip(p, wires=0)
         cirq_ops.BitFlip(p, wires=1)
 
@@ -261,13 +261,13 @@ def bell_pair():
 #           situations.
 
 
-@qml.set_shots(1000)
-@qml.qnode(dev)
+@qp.set_shots(1000)
+@qp.qnode(dev)
 def circuit(gate_params, noise_param=0.0):
-    qml.RX(gate_params[0], wires=0)
-    qml.RY(gate_params[1], wires=0)
+    qp.RX(gate_params[0], wires=0)
+    qp.RY(gate_params[1], wires=0)
     cirq_ops.Depolarize(noise_param, wires=0)
-    return qml.expval(qml.PauliZ(0))
+    return qp.expval(qp.PauliZ(0))
 
 
 gate_pars = [0.54, 0.12]
@@ -316,7 +316,7 @@ def noisy_cost(x):
 
 
 # initialize the optimizer
-opt = qml.GradientDescentOptimizer(stepsize=0.4)
+opt = qp.GradientDescentOptimizer(stepsize=0.4)
 
 # set the number of steps
 steps = 100
diff --git a/demonstrations_v2/tutorial_noisy_circuit_optimization/metadata.json b/demonstrations_v2/tutorial_noisy_circuit_optimization/metadata.json
index 23bd8bc4af..733d12d692 100644
--- a/demonstrations_v2/tutorial_noisy_circuit_optimization/metadata.json
+++ b/demonstrations_v2/tutorial_noisy_circuit_optimization/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2020-06-01T00:00:00+00:00",
-    "dateOfLastModification": "2025-10-15T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Devices and Performance"
     ],
diff --git a/demonstrations_v2/tutorial_noisy_circuits/demo.py b/demonstrations_v2/tutorial_noisy_circuits/demo.py
index 2009b347a6..d80b979d2b 100644
--- a/demonstrations_v2/tutorial_noisy_circuits/demo.py
+++ b/demonstrations_v2/tutorial_noisy_circuits/demo.py
@@ -57,7 +57,7 @@
 # Bell state :math:`|\psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle).` We ask the QNode to
 # return the expectation value of :math:`Z_0\otimes Z_1:`
 #
-import pennylane as qml
+import pennylane as qp
 from jax import numpy as np
 import jax
 import jaxopt
@@ -65,13 +65,13 @@
 jax.config.update("jax_platform_name", "cpu")
 jax.config.update('jax_enable_x64', True)
 
-dev = qml.device('default.mixed', wires=2)
+dev = qp.device('default.mixed', wires=2)
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def circuit():
-    qml.Hadamard(wires=0)
-    qml.CNOT(wires=[0, 1])
-    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))
+    qp.Hadamard(wires=0)
+    qp.CNOT(wires=[0, 1])
+    return qp.expval(qp.PauliZ(0) @ qp.PauliZ(1))
 
 
 print(f"QNode output = {circuit():.4f}")
@@ -81,11 +81,11 @@ def circuit():
 # equal to :math:`|\psi\rangle\langle\psi|,`
 # where :math:`|\psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle).`
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def density_matrix_circuit():
-    qml.Hadamard(wires=0)
-    qml.CNOT(wires=[0, 1])
-    return qml.state()
+    qp.Hadamard(wires=0)
+    qp.CNOT(wires=[0, 1])
+    return qp.state()
 
 matrix = density_matrix_circuit()
 print(f"Output density matrix is = \n{np.real(matrix)}")
@@ -128,7 +128,7 @@ def density_matrix_circuit():
 #     K_0 &= \sqrt{1-p}\begin{pmatrix}1 & 0\\ 0 & 1\end{pmatrix}, \\
 #     K_1 &= \sqrt{p}\begin{pmatrix}0 & 1\\ 1 & 0\end{pmatrix}.
 #
-# This channel can be implemented in PennyLane using the :class:`qml.BitFlip <pennylane.BitFlip>`
+# This channel can be implemented in PennyLane using the :class:`qp.BitFlip <pennylane.BitFlip>`
 # operation.
 #
 # Let's see what happens when we simulate this type of noise acting on
@@ -136,13 +136,13 @@ def density_matrix_circuit():
 #
 
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def bitflip_circuit(p):
-    qml.Hadamard(wires=0)
-    qml.CNOT(wires=[0, 1])
-    qml.BitFlip(p, wires=0)
-    qml.BitFlip(p, wires=1)
-    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1)), qml.state()
+    qp.Hadamard(wires=0)
+    qp.CNOT(wires=[0, 1])
+    qp.BitFlip(p, wires=0)
+    qp.BitFlip(p, wires=1)
+    return qp.expval(qp.PauliZ(0) @ qp.PauliZ(1)), qp.state()
 
 
 ps = [0.001, 0.01, 0.1, 0.2]
@@ -185,13 +185,13 @@ def bitflip_circuit(p):
 # below.
 
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def depolarizing_circuit(p):
-    qml.Hadamard(wires=0)
-    qml.CNOT(wires=[0, 1])
-    qml.DepolarizingChannel(p, wires=0)
-    qml.DepolarizingChannel(p, wires=1)
-    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))
+    qp.Hadamard(wires=0)
+    qp.CNOT(wires=[0, 1])
+    qp.DepolarizingChannel(p, wires=0)
+    qp.DepolarizingChannel(p, wires=1)
+    return qp.expval(qp.PauliZ(0) @ qp.PauliZ(1))
 
 
 ps = [0.001, 0.01, 0.1, 0.2]
@@ -245,13 +245,13 @@ def depolarizing_circuit(p):
 def sigmoid(x):
     return 1/(1+np.exp(-x))
 
-@qml.qnode(dev)
+@qp.qnode(dev)
 def damping_circuit(x):
-    qml.Hadamard(wires=0)
-    qml.CNOT(wires=[0, 1])
-    qml.AmplitudeDamping(sigmoid(x), wires=0)  # p = sigmoid(x)
-    qml.AmplitudeDamping(sigmoid(x), wires=1)
-    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))
+    qp.Hadamard(wires=0)
+    qp.CNOT(wires=[0, 1])
+    qp.AmplitudeDamping(sigmoid(x), wires=0)  # p = sigmoid(x)
+    qp.AmplitudeDamping(sigmoid(x), wires=1)
+    return qp.expval(qp.PauliZ(0) @ qp.PauliZ(1))
 
 ######################################################################
 # We optimize the circuit with respect to a simple cost function that attains its minimum when
diff --git a/demonstrations_v2/tutorial_noisy_circuits/metadata.json b/demonstrations_v2/tutorial_noisy_circuits/metadata.json
index f3881691b8..1db86074e2 100644
--- a/demonstrations_v2/tutorial_noisy_circuits/metadata.json
+++ b/demonstrations_v2/tutorial_noisy_circuits/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2021-02-22T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Getting Started"
     ],
diff --git a/demonstrations_v2/tutorial_odegen/demo.py b/demonstrations_v2/tutorial_odegen/demo.py
index c7b4cc1e21..3c7232729a 100644
--- a/demonstrations_v2/tutorial_odegen/demo.py
+++ b/demonstrations_v2/tutorial_odegen/demo.py
@@ -164,7 +164,7 @@
 
 """
 
-import pennylane as qml
+import pennylane as qp
 import numpy as np
 import jax.numpy as jnp
 import jax
@@ -176,7 +176,7 @@
 
 import matplotlib.pyplot as plt
 
-H_obj = qml.sum(qml.PauliX(0) @ qml.PauliX(1), qml.PauliY(0) @ qml.PauliY(1), qml.PauliZ(0) @ qml.PauliZ(1))
+H_obj = qp.sum(qp.PauliX(0) @ qp.PauliX(1), qp.PauliY(0) @ qp.PauliY(1), qp.PauliZ(0) @ qp.PauliZ(1))
 E_exact = -3.
 wires = H_obj.wires
 
@@ -209,17 +209,17 @@ def drive_field(T, wdrive):
     """Set the evolution time ``T`` and drive frequency ``wdrive``"""
     def wrapped(p, t):
         # The first len(p) values of the trainable params p characterize the pwc function
-        amp = qml.pulse.pwc(T)(p[:len(p)//2], t)
-        phi = qml.pulse.pwc(T)(p[len(p)//2:], t)
+        amp = qp.pulse.pwc(T)(p[:len(p)//2], t)
+        phi = qp.pulse.pwc(T)(p[len(p)//2:], t)
         return amp * jnp.sin(wdrive * t + phi)
 
     return wrapped
 
-H_pulse = qml.dot(-0.5*qubit_freq, [qml.PauliZ(i) for i in wires])
-H_pulse += g * (qml.PauliX(wires[0]) @ qml.PauliX(wires[1]) + qml.PauliY(wires[0]) @ qml.PauliY(wires[1]))
+H_pulse = qp.dot(-0.5*qubit_freq, [qp.PauliZ(i) for i in wires])
+H_pulse += g * (qp.PauliX(wires[0]) @ qp.PauliX(wires[1]) + qp.PauliY(wires[0]) @ qp.PauliY(wires[1]))
 
-H_pulse += drive_field(T_CR, qubit_freq[0]) * qml.PauliY(wires[0]) # on-resonance driving of qubit 0
-H_pulse += drive_field(T_CR, qubit_freq[0]) * qml.PauliY(wires[1]) # off-resonance driving of qubit 1
+H_pulse += drive_field(T_CR, qubit_freq[0]) * qp.PauliY(wires[0]) # on-resonance driving of qubit 0
+H_pulse += drive_field(T_CR, qubit_freq[0]) * qp.PauliY(wires[1]) # off-resonance driving of qubit 1
 
 ##############################################################################
 # We can now define the cost function that computes the expectation value of
@@ -227,26 +227,26 @@ def wrapped(p, t):
 # We then define the two separate qnodes with ODEgen and SPS as their differentiation methods, respectively.
 
 def qnode0(params):
-    qml.evolve(H_pulse)((params[0], params[1]), t=T_CR)
-    return qml.expval(H_obj)
+    qp.evolve(H_pulse)((params[0], params[1]), t=T_CR)
+    return qp.expval(H_obj)
 
-dev = qml.device("default.qubit", wires=range(2))
+dev = qp.device("default.qubit", wires=range(2))
 
-qnode_jax = qml.QNode(qnode0, dev, interface="jax")
+qnode_jax = qp.QNode(qnode0, dev, interface="jax")
 value_and_grad_jax = jax.jit(jax.value_and_grad(qnode_jax))
 
 num_split_times = 8
 gradient_kwargs = {"use_broadcasting": True, "num_split_times": num_split_times}
-qnode_sps = qml.QNode(
+qnode_sps = qp.QNode(
     qnode0,
     dev,
     interface="jax",
-    diff_method=qml.gradients.stoch_pulse_grad,
+    diff_method=qp.gradients.stoch_pulse_grad,
     gradient_kwargs=gradient_kwargs,
 )
 value_and_grad_sps = jax.value_and_grad(qnode_sps)
 
-qnode_odegen = qml.QNode(qnode0, dev, interface="jax", diff_method=qml.gradients.pulse_odegen)
+qnode_odegen = qp.QNode(qnode0, dev, interface="jax", diff_method=qp.gradients.pulse_odegen)
 value_and_grad_odegen = jax.value_and_grad(qnode_odegen)
 
 ##############################################################################
diff --git a/demonstrations_v2/tutorial_odegen/metadata.json b/demonstrations_v2/tutorial_odegen/metadata.json
index 607e117847..047482eaca 100644
--- a/demonstrations_v2/tutorial_odegen/metadata.json
+++ b/demonstrations_v2/tutorial_odegen/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2023-12-12T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Optimization",
         "Quantum Computing",