diff --git a/demonstrations_v2/tutorial_vqls/demo.py b/demonstrations_v2/tutorial_vqls/demo.py
index 28cceff0e1..32e1e68fd5 100644
--- a/demonstrations_v2/tutorial_vqls/demo.py
+++ b/demonstrations_v2/tutorial_vqls/demo.py
@@ -183,7 +183,7 @@
 """
 
 # Pennylane
-import pennylane as qml
+import pennylane as qp
 from pennylane import numpy as np
 
 # Plotting
@@ -219,7 +219,7 @@
 def U_b():
     """Unitary matrix rotating the ground state to the problem vector |b> = U_b |0>."""
     for idx in range(n_qubits):
-        qml.Hadamard(wires=idx)
+        qp.Hadamard(wires=idx)
 
 def CA(idx):
     """Controlled versions of the unitary components A_l of the problem matrix A."""
@@ -228,11 +228,11 @@ def CA(idx):
         None
 
     elif idx == 1:
-        qml.CNOT(wires=[ancilla_idx, 0])
-        qml.CZ(wires=[ancilla_idx, 1])
+        qp.CNOT(wires=[ancilla_idx, 0])
+        qp.CZ(wires=[ancilla_idx, 1])
 
     elif idx == 2:
-        qml.CNOT(wires=[ancilla_idx, 0])
+        qp.CNOT(wires=[ancilla_idx, 0])
 
 
 ##############################################################################
@@ -256,11 +256,11 @@ def variational_block(weights):
     """Variational circuit mapping the ground state |0> to the ansatz state |x>."""
     # We first prepare an equal superposition of all the states of the computational basis.
     for idx in range(n_qubits):
-        qml.Hadamard(wires=idx)
+        qp.Hadamard(wires=idx)
 
     # A very minimal variational circuit.
     for idx, element in enumerate(weights):
-        qml.RY(element, wires=idx)
+        qp.RY(element, wires=idx)
 
 
 ##############################################################################
@@ -276,18 +276,18 @@ def variational_block(weights):
 # and will be used to estimate the coefficients :math:`\mu_{l,l',j}` defined in the introduction.
 # A graphical representation of this circuit is shown at the top of this tutorial.
 
-dev_mu = qml.device("lightning.qubit", wires=tot_qubits)
+dev_mu = qp.device("lightning.qubit", wires=tot_qubits)
 
-@qml.qnode(dev_mu, interface="autograd")
+@qp.qnode(dev_mu, interface="autograd")
 def local_hadamard_test(weights, l=None, lp=None, j=None, part=None):
 
     # First Hadamard gate applied to the ancillary qubit.
-    qml.Hadamard(wires=ancilla_idx)
+    qp.Hadamard(wires=ancilla_idx)
 
     # For estimating the imaginary part of the coefficient "mu", we must add a "-i"
     # phase gate.
     if part == "Im" or part == "im":
-        qml.PhaseShift(-np.pi / 2, wires=ancilla_idx)
+        qp.PhaseShift(-np.pi / 2, wires=ancilla_idx)
 
     # Variational circuit generating a guess for the solution vector |x>
     variational_block(weights)
@@ -301,7 +301,7 @@ def local_hadamard_test(weights, l=None, lp=None, j=None, part=None):
 
     # Controlled Z operator at position j. If j = -1, apply the identity.
     if j != -1:
-        qml.CZ(wires=[ancilla_idx, j])
+        qp.CZ(wires=[ancilla_idx, j])
 
     # Unitary U_b associated to the problem vector |b>.
     U_b()
@@ -311,10 +311,10 @@ def local_hadamard_test(weights, l=None, lp=None, j=None, part=None):
     CA(lp)
 
     # Second Hadamard gate applied to the ancillary qubit.
-    qml.Hadamard(wires=ancilla_idx)
+    qp.Hadamard(wires=ancilla_idx)
 
     # Expectation value of Z for the ancillary qubit.
-    return qml.expval(qml.PauliZ(wires=ancilla_idx))
+    return qp.expval(qp.PauliZ(wires=ancilla_idx))
 
 
 ##############################################################################################
@@ -382,7 +382,7 @@ def cost_loc(weights):
 
 ##############################################################################
 # To minimize the cost function we use the gradient-descent optimizer.
-opt = qml.GradientDescentOptimizer(eta)
+opt = qp.GradientDescentOptimizer(eta)
 
 
 ##############################################################################
@@ -463,10 +463,10 @@ def cost_loc(weights):
 # For this task, we initialize a new PennyLane device and define the associated
 # *qnode* circuit.
 
-dev_x = qml.device("lightning.qubit", wires=n_qubits)
+dev_x = qp.device("lightning.qubit", wires=n_qubits)
 
-@qml.set_shots(n_shots)
-@qml.qnode(dev_x, interface="autograd")
+@qp.set_shots(n_shots)
+@qp.qnode(dev_x, interface="autograd")
 def prepare_and_sample(weights):
 
     # Variational circuit generating a guess for the solution vector |x>
@@ -475,7 +475,7 @@ def prepare_and_sample(weights):
     # We assume that the system is measured in the computational basis.
     # then sampling the device will give us a value of 0 or 1 for each qubit (n_qubits)
     # this will be repeated for the total number of shots provided (n_shots)
-    return qml.sample()
+    return qp.sample()
 
 
 ##############################################################################
diff --git a/demonstrations_v2/tutorial_vqls/metadata.json b/demonstrations_v2/tutorial_vqls/metadata.json
index e51c4d7111..a7d97e876e 100644
--- a/demonstrations_v2/tutorial_vqls/metadata.json
+++ b/demonstrations_v2/tutorial_vqls/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2019-11-04T00: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_vqt/demo.py b/demonstrations_v2/tutorial_vqt/demo.py
index 69dd331ac6..fcc3618823 100644
--- a/demonstrations_v2/tutorial_vqt/demo.py
+++ b/demonstrations_v2/tutorial_vqt/demo.py
@@ -107,7 +107,7 @@
 #
 
 
-import pennylane as qml
+import pennylane as qp
 from matplotlib import pyplot as plt
 import numpy as np
 import scipy
@@ -148,9 +148,9 @@ def create_hamiltonian_matrix(n, graph):
         x = y = z = 1
         for j in range(0, n):
             if j == i[0] or j == i[1]:
-                x = np.kron(x, qml.matrix(qml.PauliX(0)))
-                y = np.kron(y, qml.matrix(qml.PauliY(0)))
-                z = np.kron(z, qml.matrix(qml.PauliZ(0)))
+                x = np.kron(x, qp.matrix(qp.PauliX(0)))
+                y = np.kron(y, qp.matrix(qp.PauliY(0)))
+                z = np.kron(z, qp.matrix(qp.PauliZ(0)))
             else:
                 x = np.kron(x, np.identity(2))
                 y = np.kron(y, np.identity(2))
@@ -237,7 +237,7 @@ def single_rotation(phi_params, qubits):
 
     rotations = ["Z", "Y", "X"]
     for i in range(0, len(rotations)):
-        qml.AngleEmbedding(phi_params[i], wires=qubits, rotation=rotations[i])
+        qp.AngleEmbedding(phi_params[i], wires=qubits, rotation=rotations[i])
 
 
 ######################################################################
@@ -254,17 +254,17 @@ def CRX_ring(parameters, wires):
     n_wires = len(wires)
 
     for param, w in zip(parameters, wires):
-        qml.CRX(param, wires=[w % n_wires, (w + 1) % n_wires])
+        qp.CRX(param, wires=[w % n_wires, (w + 1) % n_wires])
 
 
 depth = 4
-dev = qml.device("lightning.qubit", wires=nr_qubits)
+dev = qp.device("lightning.qubit", wires=nr_qubits)
 
 
 def quantum_circuit(rotation_params, coupling_params, sample=None, return_state=False):
 
     # Prepares the initial basis state corresponding to the sample
-    qml.BasisState(sample, wires=range(nr_qubits))
+    qp.BasisState(sample, wires=range(nr_qubits))
 
     # Prepares the variational ansatz for the circuit
     for i in range(0, depth):
@@ -272,14 +272,14 @@ def quantum_circuit(rotation_params, coupling_params, sample=None, return_state=
         CRX_ring(coupling_params[i], list(range(nr_qubits)))
 
     if return_state:
-        return qml.state()
+        return qp.state()
 
     # Calculates the expectation value of the Hamiltonian with respect to the prepared states
-    return qml.expval(qml.Hermitian(ham_matrix, wires=range(nr_qubits)))
+    return qp.expval(qp.Hermitian(ham_matrix, wires=range(nr_qubits)))
 
 
 # Constructs the QNode
-qnode = qml.QNode(quantum_circuit, dev, interface="autograd")
+qnode = qp.QNode(quantum_circuit, dev, interface="autograd")
 
 
 ######################################################################
@@ -291,7 +291,7 @@ def quantum_circuit(rotation_params, coupling_params, sample=None, return_state=
 rotation_params = [[[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]] for i in range(0, depth)]
 coupling_params = [[1, 1, 1, 1] for i in range(0, depth)]
 print(
-    qml.draw(qnode, level="device", show_matrices=True)(
+    qp.draw(qnode, level="device", show_matrices=True)(
         rotation_params, coupling_params, sample=[1, 0, 1, 0]
     )
 )
diff --git a/demonstrations_v2/tutorial_vqt/metadata.json b/demonstrations_v2/tutorial_vqt/metadata.json
index f1aa07f6c3..8b03fcda50 100644
--- a/demonstrations_v2/tutorial_vqt/metadata.json
+++ b/demonstrations_v2/tutorial_vqt/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2020-07-07T00:00:00+00:00",
-    "dateOfLastModification": "2025-09-22T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Optimization"
     ],
diff --git a/demonstrations_v2/tutorial_xas/demo.py b/demonstrations_v2/tutorial_xas/demo.py
index abb6336cb8..c67d5a9a52 100644
--- a/demonstrations_v2/tutorial_xas/demo.py
+++ b/demonstrations_v2/tutorial_xas/demo.py
@@ -167,12 +167,12 @@
 # Since we will be using PennyLane for other aspects of this calculation, we want to make sure the molecular orbital coefficients are consistent between our PennyLane and PySCF calculations. 
 # To do this, we can obtain the molecular orbital coefficients from PennyLane using :func:`~.pennylane.qchem.scf`, and change the coefficients in the ``hf`` instance to match.
 
-import pennylane as qml
+import pennylane as qp
 
-mole = qml.qchem.Molecule(symbols, geometry, basis_name="sto-3g", unit="angstrom")
+mole = qp.qchem.Molecule(symbols, geometry, basis_name="sto-3g", unit="angstrom")
 
 # Run self-consistent fields method to get molecular orbital coefficients.
-_, coeffs, _, _, _ = qml.qchem.scf(mole)()
+_, coeffs, _, _, _ = qp.qchem.scf(mole)()
 
 # Change MO coefficients in hf object to PennyLane calculated values.
 hf.mo_coeff = coeffs  
@@ -258,13 +258,13 @@
 # This is fine for our simple example, but for more complex instances you may want to modify the dipole operator to restrict the final states (see the Appendix for more details).
 
 # Get core and active orbital indices.
-core, active = qml.qchem.active_space(
+core, active = qp.qchem.active_space(
     mole.n_electrons, mole.n_orbitals, 
     active_electrons=n_electron_cas, 
     active_orbitals=n_cas
 )
 
-m_rho = qml.qchem.dipole_moment(mole, cutoff=1e-8, core=core, active=active)()
+m_rho = qp.qchem.dipole_moment(mole, cutoff=1e-8, core=core, active=active)()
 rhos = range(len(m_rho))  # [0, 1, 2] are [x, y, z].
 
 wf_dipole = []
@@ -272,7 +272,7 @@
 
 # Loop over cartesian coordinates and calculate m_rho|I>.
 for rho in rhos:
-    dipole_matrix_rho = qml.matrix(m_rho[rho], wire_order=range(2 * n_cas))
+    dipole_matrix_rho = qp.matrix(m_rho[rho], wire_order=range(2 * n_cas))
     wf = dipole_matrix_rho.dot(wf_casci)  # Multiply state into dipole matrix.
 
     if np.allclose(wf, np.zeros_like(wf)):  # If wf is zero, then set norm as zero.
@@ -294,15 +294,15 @@
 # We will also add one auxiliary wire for the measurement circuit, which we will prepare as the 0 wire with an applied Hadamard gate.
 
 device_type = "lightning.qubit"
-dev_prop = qml.device(device_type, wires=int(2*n_cas) + 1)
+dev_prop = qp.device(device_type, wires=int(2*n_cas) + 1)
 
 
-@qml.qnode(dev_prop)
+@qp.qnode(dev_prop)
 def initial_circuit(wf):
     """Circuit to load initial state and prepare auxiliary qubit."""
-    qml.StatePrep(wf, wires=dev_prop.wires.tolist()[1:])
-    qml.Hadamard(wires=0)
-    return qml.state()
+    qp.StatePrep(wf, wires=dev_prop.wires.tolist()[1:])
+    qp.Hadamard(wires=0)
+    return qp.state()
 
 
 ######################################################################
@@ -331,14 +331,14 @@ def initial_circuit(wf):
 #
 # The core constant and the one- and two-electron integrals can be computed in PennyLane using :func:`~pennylane.qchem.electron_integrals`.
 
-core_constant, one, two = qml.qchem.electron_integrals(mole, core=core, active=active)()
+core_constant, one, two = qp.qchem.electron_integrals(mole, core=core, active=active)()
 core_constant = core_constant[0]
 
 ######################################################################
 # We will have to convert these to chemists' notation [#Sherrill2005]_.
 
-two_chemist = qml.math.einsum("prsq->pqrs", two)
-one_chemist = one - qml.math.einsum("pqrr->pq", two) / 2.0
+two_chemist = qp.math.einsum("prsq->pqrs", two)
+one_chemist = one - qp.math.einsum("pqrr->pq", two) / 2.0
 
 ######################################################################
 # Next, we will perform compressed double factorization of the Hamiltonian's two-electron integrals to approximate them as a sum of :math:`L` fragments
@@ -355,7 +355,7 @@ def initial_circuit(wf):
 config.update("jax_enable_x64", True)  # Required for factorize consistency.
 
 # Factorize Hamiltonian, producing matrices Z and U for each fragment.
-_, Z, U = qml.qchem.factorize(two_chemist, compressed=True)
+_, Z, U = qp.qchem.factorize(two_chemist, compressed=True)
 
 print("Shape of the factors: ")
 print("two_chemist", two_chemist.shape)
@@ -363,8 +363,8 @@ def initial_circuit(wf):
 print("Z", Z.shape)
 
 # Compare factorized two-electron fragment sum to the original.
-approx_two_chemist = qml.math.einsum("tpk,tqk,tkl,trl,tsl->pqrs", U, U, Z, U, U)
-assert qml.math.allclose(approx_two_chemist, two_chemist, atol=0.1)
+approx_two_chemist = qp.math.einsum("tpk,tqk,tkl,trl,tsl->pqrs", U, U, Z, U, U)
+assert qp.math.allclose(approx_two_chemist, two_chemist, atol=0.1)
 
 ######################################################################
 # Note there are some terms in this decomposition of the two-electron integrals that are exactly diagonalizable, and can be combined with the one-electron integrals to simplify how the Hamiltonian time evolution is implemented. 
@@ -373,7 +373,7 @@ def initial_circuit(wf):
 
 # Calculate the one-electron extra terms.
 Z_prime = np.stack([np.diag(np.sum(Z[i], axis=-1)) for i in range(Z.shape[0])], axis=0)
-one_electron_extra = qml.math.einsum("tpk,tkk,tqk->pq", U, Z_prime, U)
+one_electron_extra = qp.math.einsum("tpk,tkk,tqk->pq", U, Z_prime, U)
 
 # Diagonalize the one-electron integral matrix while adding the one-electron extra.
 eigenvals, U0 = np.linalg.eigh(one_chemist + one_electron_extra)
@@ -411,8 +411,8 @@ def initial_circuit(wf):
 
 def U_rotations(U, control_wires):
     """Circuit implementing the basis rotations of the Hamiltonian fragments."""
-    U_spin = qml.math.kron(U, qml.math.eye(2))  # Apply to both spins.
-    qml.BasisRotation(
+    U_spin = qp.math.kron(U, qp.math.eye(2))  # Apply to both spins.
+    qp.BasisRotation(
         unitary_matrix=U_spin, wires=[int(i + control_wires) for i in range(2 * n_cas)]
     )
 
@@ -439,14 +439,14 @@ def Z_rotations(Z, step, is_one_electron_term, control_wires):
     if is_one_electron_term:
         for sigma in range(2):
             for i in range(n_cas):
-                qml.ctrl(
-                    qml.X(wires=int(2*i + sigma + control_wires)),
+                qp.ctrl(
+                    qp.X(wires=int(2*i + sigma + control_wires)),
                     control=range(control_wires),
                     control_values=0,
                 )
-                qml.RZ(-Z[i, i] * step / 2, wires=int(2*i + sigma + control_wires))
-                qml.ctrl(
-                    qml.X(wires=int(2*i + sigma + control_wires)),
+                qp.RZ(-Z[i, i] * step / 2, wires=int(2*i + sigma + control_wires))
+                qp.ctrl(
+                    qp.X(wires=int(2*i + sigma + control_wires)),
                     control=range(control_wires),
                     control_values=0,
                 )
@@ -456,16 +456,16 @@ def Z_rotations(Z, step, is_one_electron_term, control_wires):
         for sigma, tau in product(range(2), repeat=2):
             for i, k in product(range(n_cas), repeat=2):
                 if i != k or sigma != tau:  # Skip the one-electron correction terms.
-                    qml.ctrl(qml.X(wires=int(2*i + sigma + control_wires)),
+                    qp.ctrl(qp.X(wires=int(2*i + sigma + control_wires)),
                             control=range(control_wires), control_values=0)
-                    qml.MultiRZ(Z[i, k] / 8.0 * step,
+                    qp.MultiRZ(Z[i, k] / 8.0 * step,
                         wires=[int(2*i + sigma + control_wires),
                                int(2*k + tau + control_wires)])
-                    qml.ctrl(qml.X(wires=int(2 * i + sigma + control_wires)),
+                    qp.ctrl(qp.X(wires=int(2 * i + sigma + control_wires)),
                         control=range(control_wires), control_values=0)
         globalphase = np.trace(Z)/4.0*step - np.sum(Z)*step/2.0
 
-    qml.PhaseShift(-globalphase, wires=0)
+    qp.PhaseShift(-globalphase, wires=0)
 
 
 ######################################################################
@@ -503,7 +503,7 @@ def first_order_trotter(step, prior_U, final_rotation, reverse=False):
         U_rotations(prior_U, 1)
 
     # Global phase adjustment from core constant.
-    qml.PhaseShift(-core_constant * step, wires=0)
+    qp.PhaseShift(-core_constant * step, wires=0)
 
     return prior_U
 
@@ -519,7 +519,7 @@ def second_order_trotter(dev, state, step):
 
     def circuit():
         # State preparation -- set as the final state from the previous iteration.
-        qml.StatePrep(state, wires=qubits)
+        qp.StatePrep(state, wires=qubits)
 
         prior_U = np.eye(n_cas)  # No initial prior U, so set as identity matrix.
         prior_U = first_order_trotter(step / 2, prior_U=prior_U, 
@@ -527,9 +527,9 @@ def circuit():
         prior_U = first_order_trotter(step / 2, prior_U=prior_U, 
                                       final_rotation=True, reverse=True)
 
-        return qml.state()
+        return qp.state()
 
-    return qml.QNode(circuit, dev)
+    return qp.QNode(circuit, dev)
 
 
 ######################################################################
@@ -541,9 +541,9 @@ def circuit():
 
 
 def meas_circuit(state):
-    qml.StatePrep(state, wires=range(int(2*n_cas) + 1))
+    qp.StatePrep(state, wires=range(int(2*n_cas) + 1))
     # Measure in PauliX/PauliY to get the real/imaginary parts.
-    return [qml.expval(op) for op in [qml.PauliX(wires=0), qml.PauliY(wires=0)]]
+    return [qp.expval(op) for op in [qp.PauliX(wires=0), qp.PauliY(wires=0)]]
 
 
 ######################################################################
@@ -630,11 +630,11 @@ def L_j(t_j):
 
         # Define measurement circuit device with shots.
         shots = shots_list[i]  # Kernel-aware number of shots.
-        dev_est = qml.device(device_type, wires=int(2*n_cas) + 1)
+        dev_est = qp.device(device_type, wires=int(2*n_cas) + 1)
 
         # Update state and then measure expectation values.
         state = circuit()
-        measurement = qml.QNode(meas_circuit, dev_est, shots=shots)(state)
+        measurement = qp.QNode(meas_circuit, dev_est, shots=shots)(state)
 
         expvals[:, i] += dipole_norm[rho]**2 * np.array(measurement).real
 
@@ -706,7 +706,7 @@ def L_j(t_j):
 mo_coeffs = coeffs[:, n_core : n_core + n_cas]
 
 # Convert to molecular orbital basis.
-dip_ints_mo = qml.math.einsum("ik,xkl,lj->xij", mo_coeffs.T, dip_ints_ao, mo_coeffs)
+dip_ints_mo = qp.math.einsum("ik,xkl,lj->xij", mo_coeffs.T, dip_ints_ao, mo_coeffs)
 
 
 def final_state_overlap(ci_id):
@@ -715,7 +715,7 @@ def final_state_overlap(ci_id):
         mycasci.ci[0], mycasci.ci[ci_id], n_cas, n_electron_cas
     )
     # Transition dipole moments.
-    return qml.math.einsum("xij,ji->x", dip_ints_mo, t_dm1)
+    return qp.math.einsum("xij,ji->x", dip_ints_mo, t_dm1)
 
 
 # Compute overlaps.
diff --git a/demonstrations_v2/tutorial_xas/metadata.json b/demonstrations_v2/tutorial_xas/metadata.json
index a078eb70cb..0ce7d069ab 100644
--- a/demonstrations_v2/tutorial_xas/metadata.json
+++ b/demonstrations_v2/tutorial_xas/metadata.json
@@ -11,7 +11,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2025-08-29T09:00:00+00:00",
-    "dateOfLastModification": "2025-10-15T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Quantum Chemistry",
         "Algorithms"
diff --git a/demonstrations_v2/tutorial_zx_calculus/demo.py b/demonstrations_v2/tutorial_zx_calculus/demo.py
index d8a292adb3..6b6a5b28ce 100644
--- a/demonstrations_v2/tutorial_zx_calculus/demo.py
+++ b/demonstrations_v2/tutorial_zx_calculus/demo.py
@@ -513,19 +513,19 @@
 
 import matplotlib.pyplot as plt
 
-import pennylane as qml
+import pennylane as qp
 import pyzx
 
-dev = qml.device("default.qubit", wires=2)
+dev = qp.device("default.qubit", wires=2)
 
 
-@qml.transforms.to_zx
-@qml.qnode(device=dev)
+@qp.transforms.to_zx
+@qp.qnode(device=dev)
 def circuit():
-    qml.PauliX(wires=0),
-    qml.PauliY(wires=1),
-    qml.CNOT(wires=[0, 1]),
-    return qml.expval(qml.PauliZ(wires=0))
+    qp.PauliX(wires=0),
+    qp.PauliY(wires=1),
+    qp.CNOT(wires=[0, 1]),
+    return qp.expval(qp.PauliZ(wires=0))
 
 
 g = circuit()
@@ -572,7 +572,7 @@ def circuit():
 
 graph = random_circuit.to_graph()
 
-tape = qml.transforms.from_zx(graph)
+tape = qp.transforms.from_zx(graph)
 print(tape.operations)
 
 #############################################################################
@@ -629,76 +629,76 @@ def circuit():
 #
 
 
-dev = qml.device("default.qubit", wires=5)
+dev = qp.device("default.qubit", wires=5)
 
 
-@qml.transforms.to_zx
-@qml.qnode(device=dev)
+@qp.transforms.to_zx
+@qp.qnode(device=dev)
 def mod_5_4():
-    qml.PauliX(wires=4),
-    qml.Hadamard(wires=4),
-    qml.CNOT(wires=[3, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[0, 4]),
-    qml.T(wires=[4]),
-    qml.CNOT(wires=[3, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[0, 4]),
-    qml.T(wires=[3]),
-    qml.T(wires=[4]),
-    qml.CNOT(wires=[0, 3]),
-    qml.T(wires=[0]),
-    qml.adjoint(qml.T(wires=[3]))
-    qml.CNOT(wires=[0, 3]),
-    qml.CNOT(wires=[3, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[2, 4]),
-    qml.T(wires=[4]),
-    qml.CNOT(wires=[3, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[2, 4]),
-    qml.T(wires=[3]),
-    qml.T(wires=[4]),
-    qml.CNOT(wires=[2, 3]),
-    qml.T(wires=[2]),
-    qml.adjoint(qml.T(wires=[3]))
-    qml.CNOT(wires=[2, 3]),
-    qml.Hadamard(wires=[4]),
-    qml.CNOT(wires=[3, 4]),
-    qml.Hadamard(wires=4),
-    qml.CNOT(wires=[2, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[1, 4]),
-    qml.T(wires=[4]),
-    qml.CNOT(wires=[2, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[1, 4]),
-    qml.T(wires=[4]),
-    qml.T(wires=[2]),
-    qml.CNOT(wires=[1, 2]),
-    qml.T(wires=[1]),
-    qml.adjoint(qml.T(wires=[2]))
-    qml.CNOT(wires=[1, 2]),
-    qml.Hadamard(wires=[4]),
-    qml.CNOT(wires=[2, 4]),
-    qml.Hadamard(wires=4),
-    qml.CNOT(wires=[1, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[0, 4]),
-    qml.T(wires=[4]),
-    qml.CNOT(wires=[1, 4]),
-    qml.adjoint(qml.T(wires=[4])),
-    qml.CNOT(wires=[0, 4]),
-    qml.T(wires=[4]),
-    qml.T(wires=[1]),
-    qml.CNOT(wires=[0, 1]),
-    qml.T(wires=[0]),
-    qml.adjoint(qml.T(wires=[1])),
-    qml.CNOT(wires=[0, 1]),
-    qml.Hadamard(wires=[4]),
-    qml.CNOT(wires=[1, 4]),
-    qml.CNOT(wires=[0, 4]),
-    return qml.expval(qml.PauliZ(wires=0))
+    qp.PauliX(wires=4),
+    qp.Hadamard(wires=4),
+    qp.CNOT(wires=[3, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[0, 4]),
+    qp.T(wires=[4]),
+    qp.CNOT(wires=[3, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[0, 4]),
+    qp.T(wires=[3]),
+    qp.T(wires=[4]),
+    qp.CNOT(wires=[0, 3]),
+    qp.T(wires=[0]),
+    qp.adjoint(qp.T(wires=[3]))
+    qp.CNOT(wires=[0, 3]),
+    qp.CNOT(wires=[3, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[2, 4]),
+    qp.T(wires=[4]),
+    qp.CNOT(wires=[3, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[2, 4]),
+    qp.T(wires=[3]),
+    qp.T(wires=[4]),
+    qp.CNOT(wires=[2, 3]),
+    qp.T(wires=[2]),
+    qp.adjoint(qp.T(wires=[3]))
+    qp.CNOT(wires=[2, 3]),
+    qp.Hadamard(wires=[4]),
+    qp.CNOT(wires=[3, 4]),
+    qp.Hadamard(wires=4),
+    qp.CNOT(wires=[2, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[1, 4]),
+    qp.T(wires=[4]),
+    qp.CNOT(wires=[2, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[1, 4]),
+    qp.T(wires=[4]),
+    qp.T(wires=[2]),
+    qp.CNOT(wires=[1, 2]),
+    qp.T(wires=[1]),
+    qp.adjoint(qp.T(wires=[2]))
+    qp.CNOT(wires=[1, 2]),
+    qp.Hadamard(wires=[4]),
+    qp.CNOT(wires=[2, 4]),
+    qp.Hadamard(wires=4),
+    qp.CNOT(wires=[1, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[0, 4]),
+    qp.T(wires=[4]),
+    qp.CNOT(wires=[1, 4]),
+    qp.adjoint(qp.T(wires=[4])),
+    qp.CNOT(wires=[0, 4]),
+    qp.T(wires=[4]),
+    qp.T(wires=[1]),
+    qp.CNOT(wires=[0, 1]),
+    qp.T(wires=[0]),
+    qp.adjoint(qp.T(wires=[1])),
+    qp.CNOT(wires=[0, 1]),
+    qp.Hadamard(wires=[4]),
+    qp.CNOT(wires=[1, 4]),
+    qp.CNOT(wires=[0, 4]),
+    return qp.expval(qp.PauliZ(wires=0))
 
 
 g = mod_5_4()
@@ -770,20 +770,20 @@ def mod_5_4():
 # and which is made possible because we used `pyzx.teleport_reduce` and do not need to extract
 # the circuit.
 
-qscript_opt = qml.transforms.from_zx(g)
+qscript_opt = qp.transforms.from_zx(g)
 
-wires = qml.wires.Wires([4, 3, 0, 2, 1])
+wires = qp.wires.Wires([4, 3, 0, 2, 1])
 wires_map = dict(zip(qscript_opt.wires, wires))
-qscript_opt_reorder, processing = qml.map_wires(qscript_opt, wire_map=wires_map)
+qscript_opt_reorder, processing = qp.map_wires(qscript_opt, wire_map=wires_map)
 
-@qml.qnode(device=dev)
+@qp.qnode(device=dev)
 def mod_5_4():
     for o in processing(qscript_opt_reorder):
-        qml.apply(o)
-    return qml.expval(qml.PauliZ(wires=0))
+        qp.apply(o)
+    return qp.expval(qp.PauliZ(wires=0))
 
 
-specs = qml.specs(mod_5_4)()
+specs = qp.specs(mod_5_4)()
 
 print("Number of quantum gates:", specs["resources"].num_gates)
 print("Circuit gates:", specs["resources"].gate_types)
diff --git a/demonstrations_v2/tutorial_zx_calculus/metadata.json b/demonstrations_v2/tutorial_zx_calculus/metadata.json
index c244f4fbe8..d5c0df619c 100644
--- a/demonstrations_v2/tutorial_zx_calculus/metadata.json
+++ b/demonstrations_v2/tutorial_zx_calculus/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": true,
     "executable_latest": true,
     "dateOfPublication": "2023-06-06T00:00:00+00:00",
-    "dateOfLastModification": "2026-01-14T15:48:14+00:00",
+    "dateOfLastModification": "2026-04-17T15:48:14+00:00",
     "categories": [
         "Quantum Computing",
         "Compilation"
diff --git a/demonstrations_v2/vqe_parallel/demo.py b/demonstrations_v2/vqe_parallel/demo.py
index 1223044da9..44f515e3e9 100644
--- a/demonstrations_v2/vqe_parallel/demo.py
+++ b/demonstrations_v2/vqe_parallel/demo.py
@@ -35,7 +35,7 @@
 
 import matplotlib.pyplot as plt
 from pennylane import numpy as np
-import pennylane as qml
+import pennylane as qp
 from pennylane import qchem
 
 ##############################################################################
@@ -66,7 +66,7 @@
 # `PennyLane Datasets library <https://pennylane.ai/datasets/qchem/h2-molecule>`__:
 
 bonds = [0.5, 0.58, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1]
-datasets = qml.data.load("qchem", molname="H2", bondlength=bonds, basis="STO-3G")
+datasets = qp.data.load("qchem", molname="H2", bondlength=bonds, basis="STO-3G")
 
 ##############################################################################
 # We can now extract the qubit Hamiltonians from these datasets for each bond length:
@@ -136,8 +136,8 @@
 #
 # To do this, start by instantiating a device for each term:
 
-dev1 = [qml.device("rigetti.qvm", device="4q-qvm") for _ in range(8)]
-dev2 = [qml.device("rigetti.qvm", device="9q-square-qvm") for _ in range(7)]
+dev1 = [qp.device("rigetti.qvm", device="4q-qvm") for _ in range(8)]
+dev2 = [qp.device("rigetti.qvm", device="9q-square-qvm") for _ in range(7)]
 devs = dev1 + dev2
 
 ##############################################################################
@@ -168,12 +168,12 @@
 
 
 def circuit(param, H):
-    qml.BasisState(np.array([1, 1, 0, 0], requires_grad=False), wires=[0, 1, 2, 3])
-    qml.RY(param, wires=2)
-    qml.CNOT(wires=[2, 3])
-    qml.CNOT(wires=[2, 0])
-    qml.CNOT(wires=[3, 1])
-    return qml.expval(H)
+    qp.BasisState(np.array([1, 1, 0, 0], requires_grad=False), wires=[0, 1, 2, 3])
+    qp.RY(param, wires=2)
+    qp.CNOT(wires=[2, 3])
+    qp.CNOT(wires=[2, 0])
+    qp.CNOT(wires=[3, 1])
+    return qp.expval(H)
 
 
 ##############################################################################
@@ -200,7 +200,7 @@ def circuit(param, H):
     print(
         f"{i+1} / {len(bonds)}: Sequential execution; Running for inter-atomic distance {bonds[i]} Å"
     )
-    energies_seq.append(qml.QNode(circuit, devs[0])(param, h))
+    energies_seq.append(qp.QNode(circuit, devs[0])(param, h))
 
 dt_seq = time.time() - t0
 
@@ -218,7 +218,7 @@ def compute_energy_parallel(H, devs, param):
     results = []
 
     for i in range(len(H_ops)):
-        qnode = qml.QNode(circuit, devs[i])
+        qnode = qp.QNode(circuit, devs[i])
         results.append(dask.delayed(qnode)(param, H_ops[i]))
 
     results = dask.compute(*results, scheduler="threads")
@@ -229,7 +229,7 @@ def compute_energy_parallel(H, devs, param):
 ##############################################################################
 # We can now compute all 10 samples from the energy surface sequentially, where each execution is making use of
 # parallel device execution. Curiously, in this example the overhead from doing so outweighs the speed-up
-# and the execution is slower than standard execution using ``qml.expval``. For different circuits and
+# and the execution is slower than standard execution using ``qp.expval``. For different circuits and
 # different Hamiltonians, however, parallelization may provide significant speed-ups.
 
 print("Evaluating the potential energy surface in parallel")
@@ -259,16 +259,16 @@ def compute_energy_parallel_optimized(H, devs, param):
     assert len(H_ops) == len(devs)
     results = []
 
-    obs_groupings, coeffs_groupings = qml.pauli.group_observables(
+    obs_groupings, coeffs_groupings = qp.pauli.group_observables(
         H_ops, H_coeffs, "qwc"
     )
 
     for i, (obs, coeffs) in enumerate(zip(obs_groupings, coeffs_groupings)):
-        H_part = qml.Hamiltonian(coeffs, obs)
-        qnode = qml.QNode(circuit, devs[i])
+        H_part = qp.Hamiltonian(coeffs, obs)
+        qnode = qp.QNode(circuit, devs[i])
         results.append(dask.delayed(qnode)(param, H_part))
 
-    result = qml.math.sum(dask.compute(*results, scheduler="threads"))
+    result = qp.math.sum(dask.compute(*results, scheduler="threads"))
     return result
 
 print(
diff --git a/demonstrations_v2/vqe_parallel/metadata.json b/demonstrations_v2/vqe_parallel/metadata.json
index 99a308023e..00333508c3 100644
--- a/demonstrations_v2/vqe_parallel/metadata.json
+++ b/demonstrations_v2/vqe_parallel/metadata.json
@@ -8,7 +8,7 @@
     "executable_stable": false,
     "executable_latest": false,
     "dateOfPublication": "2020-02-14T00: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/zne_catalyst/demo.py b/demonstrations_v2/zne_catalyst/demo.py
index 1db090fed4..74d536d79e 100644
--- a/demonstrations_v2/zne_catalyst/demo.py
+++ b/demonstrations_v2/zne_catalyst/demo.py
@@ -96,7 +96,7 @@
 # equal to 1.
 
 import numpy as np
-import pennylane as qml
+import pennylane as qp
 from catalyst import mitigate_with_zne
 
 import warnings
@@ -108,21 +108,21 @@
 np.random.seed(42)
 
 n_layers = 5
-template = qml.SimplifiedTwoDesign
+template = qp.SimplifiedTwoDesign
 weights_shape = template.shape(n_layers, n_wires)
 w1, w2 = [2 * np.pi * np.random.random(s) for s in weights_shape]
 
 def circuit(w1, w2):
     template(w1, w2, wires=range(n_wires))
-    qml.adjoint(template)(w1, w2, wires=range(n_wires))
-    return qml.expval(qml.PauliZ(0))
+    qp.adjoint(template)(w1, w2, wires=range(n_wires))
+    return qp.expval(qp.PauliZ(0))
 
 ##############################################################################
 # As a sanity check, we first execute the circuit on the Qrack simulator without any noise.
 
-noiseless_device = qml.device("qrack.simulator", n_wires, noise=0)
+noiseless_device = qp.device("qrack.simulator", n_wires, noise=0)
 
-ideal_value = qml.QNode(circuit, device=noiseless_device)(w1, w2)
+ideal_value = qp.QNode(circuit, device=noiseless_device)(w1, w2)
 print(f"Ideal value: {ideal_value}")
 
 ##############################################################################
@@ -136,9 +136,9 @@ def circuit(w1, w2):
 # The probability of error is specified by the value of the ``noise`` constructor argument.
 
 NOISE_LEVEL = 0.01
-noisy_device = qml.device("qrack.simulator", n_wires, shots=1000, noise=NOISE_LEVEL)
+noisy_device = qp.device("qrack.simulator", n_wires, shots=1000, noise=NOISE_LEVEL)
 
-noisy_qnode = qml.QNode(circuit, device=noisy_device, mcm_method="one-shot")
+noisy_qnode = qp.QNode(circuit, device=noisy_device, mcm_method="one-shot")
 noisy_value = noisy_qnode(w1, w2)
 print(f"Error without mitigation: {abs(ideal_value - noisy_value):.3f}")
 
@@ -166,7 +166,7 @@ def circuit(w1, w2):
 
 ##############################################################################
 # Finally, we'll choose the extrapolation technique. Both exponential and polynomial extrapolation
-# is available in the :mod:`qml.noise <pennylane.noise>` module, and both of these functions can be passed directly
+# is available in the :mod:`qp.noise <pennylane.noise>` module, and both of these functions can be passed directly
 # into Catalyst's :func:`catalyst.mitigate_with_zne` function. In this tutorial we use polynomial extrapolation,
 # which we hypothesize best models the behavior of the noise scenario we are considering.
 
@@ -180,7 +180,7 @@ def circuit(w1, w2):
 # to define a very simple :func:`~.QNode`, which represents the mitigated version of the original circuit.
 
 
-@qml.qjit
+@qp.qjit
 def mitigated_circuit_qjit(w1, w2):
     return mitigate_with_zne(
         noisy_qnode,
@@ -209,11 +209,11 @@ def mitigated_circuit_qjit(w1, w2):
 
 
 def mitigated_circuit(w1, w2):
-    return qml.noise.mitigate_with_zne(
+    return qp.noise.mitigate_with_zne(
         noisy_qnode,
         scale_factors=scale_factors,
         extrapolate=extrapolation_method,
-        folding=qml.noise.fold_global,
+        folding=qp.noise.fold_global,
     )(w1, w2)
 
 
@@ -230,10 +230,10 @@ def mitigated_circuit(w1, w2):
 # reduce the running time of this tutorial, while still showcasing the performance differences.  
 import timeit
 
-noisy_device = qml.device("qrack.simulator", n_wires, shots=100, noise=NOISE_LEVEL)
-noisy_qnode = qml.QNode(circuit, device=noisy_device, mcm_method="one-shot")
+noisy_device = qp.device("qrack.simulator", n_wires, shots=100, noise=NOISE_LEVEL)
+noisy_qnode = qp.QNode(circuit, device=noisy_device, mcm_method="one-shot")
 
-@qml.qjit
+@qp.qjit
 def mitigated_circuit_qjit(w1, w2):
     return mitigate_with_zne(
         noisy_qnode,
diff --git a/demonstrations_v2/zne_catalyst/metadata.json b/demonstrations_v2/zne_catalyst/metadata.json
index 9376a3d40a..9211308089 100644
--- a/demonstrations_v2/zne_catalyst/metadata.json
+++ b/demonstrations_v2/zne_catalyst/metadata.json
@@ -11,7 +11,7 @@
     "executable_stable": false,
     "executable_latest": true,
     "dateOfPublication": "2024-11-15T00:00:00+00:00",
-    "dateOfLastModification": "2025-10-15T00:00:00+00:00",
+    "dateOfLastModification": "2026-04-17T00:00:00+00:00",
     "categories": [
         "Algorithms",
         "Quantum Computing"