diff --git a/doc/introduction/data.rst b/doc/introduction/data.rst
index f3bda0d32a6..64131c22144 100644
--- a/doc/introduction/data.rst
+++ b/doc/introduction/data.rst
@@ -133,8 +133,7 @@ To create a dataset, we can do the following:
 >>> dataset.data_name
 "Example"
 >>> dataset.hamiltonian
-(0.5) [X1]
-+ (1) [Z0]
+1.0 * Z(0) + 0.5 * X(1)
 >>> dataset.energies
 array([-1.5, -0.5,  0.5,  1.5])
 
@@ -147,8 +146,7 @@ We can then write this :class:`~pennylane.data.Dataset` to storage and read it a
 >>> read_dataset.data_name
 "Example"
 >>> read_dataset.hamiltonian
-(0.5) [X1]
-+ (1) [Z0]
+1.0 * Z(0) + 0.5 * X(1)
 >>> read_dataset.energies
 array([-1.5, -0.5,  0.5,  1.5])
 
diff --git a/doc/introduction/operations.rst b/doc/introduction/operations.rst
index 2be84f077ed..9666fc58184 100644
--- a/doc/introduction/operations.rst
+++ b/doc/introduction/operations.rst
@@ -132,10 +132,9 @@ For example:
 >>> mat = np.array([[1, 1], [1, -1]])
 >>> h = qml.pauli_decompose(mat)
 >>> type(h)
-<class 'pennylane.ops.qubit.hamiltonian.Hamiltonian'>
+pennylane.ops.op_math.linear_combination.LinearCombination
 >>> print(h)
-(1.0) [X0]
-+ (1.0) [Z0]
+1.0 * X(0) + 1.0 * Z(0)
 
 .. _intro_ref_ops_qubit:
 
diff --git a/pennylane/data/__init__.py b/pennylane/data/__init__.py
index c6d93d7e02d..baa5bdf75df 100644
--- a/pennylane/data/__init__.py
+++ b/pennylane/data/__init__.py
@@ -89,8 +89,7 @@
 ...   eigen = {"eigvals": eigvals, "eigvecs": eigvecs}
 ... )
 >>> dataset.hamiltonian
-(1.0) [Z0]
-+ (1.0) [Z1]
+1.0 * Z(0) + 1.0 * Z(1)
 >>> dataset.eigen
 {'eigvals': array([-2.,  0.,  0.,  2.]),
 'eigvecs': array([[0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
diff --git a/pennylane/devices/preprocess.py b/pennylane/devices/preprocess.py
index 6c753519226..7abdc3fc0b8 100644
--- a/pennylane/devices/preprocess.py
+++ b/pennylane/devices/preprocess.py
@@ -389,7 +389,7 @@ def validate_observables(
     ...    return obj.name in {"PauliX", "PauliY", "PauliZ"}
     >>> tape = qml.tape.QuantumScript([], [qml.expval(qml.Z(0) + qml.Y(0))])
     >>> validate_observables(tape, accepted_observable)
-    DeviceError: Observable <Hamiltonian: terms=2, wires=[0]> not supported on device
+    DeviceError: Observable Z(0) + Y(0) not supported on device
 
     Note that if the observable is a :class:`~.Tensor`, the validation is run on each object in the
     ``Tensor`` instead.
diff --git a/pennylane/measurements/__init__.py b/pennylane/measurements/__init__.py
index 1600c0fbfce..eb6ca2e54ec 100644
--- a/pennylane/measurements/__init__.py
+++ b/pennylane/measurements/__init__.py
@@ -227,11 +227,11 @@ def circuit(x):
     >>> H = 2.0 * qml.X(0)
     >>> mp = qml.expval(H)
     >>> mp._flatten()
-    ((<Hamiltonian: terms=1, wires=[0]>, None), ())
+    ((2.0 * X(0), None), (('wires', None),))
     >>> type(mp)._unflatten(*mp._flatten())
-    expval(  (2) [X0])
+    expval(2.0 * X(0))
     >>> jax.tree_util.tree_leaves(mp)
-    [2]
+    [2.0]
 
 Adding your new measurement to PennyLane
 ----------------------------------------
diff --git a/pennylane/ops/functions/generator.py b/pennylane/ops/functions/generator.py
index b39986731b1..2fe0062d91e 100644
--- a/pennylane/ops/functions/generator.py
+++ b/pennylane/ops/functions/generator.py
@@ -169,8 +169,12 @@ def generator(op: qml.operation.Operator, format="prefactor"):
     >>> op = qml.RX(0.2, wires=0)
     >>> qml.generator(op, format="prefactor")  # output will always be (obs, prefactor)
     (X(0), -0.5)
-    >>> qml.generator(op, format="hamiltonian")  # output will always be a Hamiltonian
-    (-0.5) [X0]
+    >>> qml.generator(op, format="hamiltonian")  # output will always be a Hamiltonian/LinearCombination
+    -0.5 * X(0)
+    >>> with qml.operation.disable_new_opmath_cm():
+    ...     gen = qml.generator(op, format="hamiltonian")) # legacy Hamiltonian class
+    ...     print(gen, type(gen))
+    (-0.5) [X0] <class 'pennylane.ops.qubit.hamiltonian.Hamiltonian'>
     >>> qml.generator(qml.PhaseShift(0.1, wires=0), format="observable")  # ouput will be a simplified obs where possible
     Projector([1], wires=[0])
     >>> qml.generator(op, format="arithmetic")  # output is an instance of `SProd`
diff --git a/pennylane/qaoa/cost.py b/pennylane/qaoa/cost.py
index bb6afeece2b..f1c5ef5c050 100644
--- a/pennylane/qaoa/cost.py
+++ b/pennylane/qaoa/cost.py
@@ -53,9 +53,7 @@ def bit_driver(wires: Union[Iterable, qaoa.Wires], b: int):
     >>> wires = range(3)
     >>> hamiltonian = qaoa.bit_driver(wires, 1)
     >>> print(hamiltonian)
-      (1) [Z0]
-    + (1) [Z1]
-    + (1) [Z2]
+    1 * Z(0) + 1 * Z(1) + 1 * Z(2)
     """
     if b == 0:
         coeffs = [-1 for _ in wires]
@@ -96,12 +94,7 @@ def edge_driver(graph: Union[nx.Graph, rx.PyGraph], reward: list):
     >>> graph = nx.Graph([(0, 1), (1, 2)])
     >>> hamiltonian = qaoa.edge_driver(graph, ["11", "10", "01"])
     >>> print(hamiltonian)
-      (0.25) [Z0]
-    + (0.25) [Z1]
-    + (0.25) [Z1]
-    + (0.25) [Z2]
-    + (0.25) [Z0 Z1]
-    + (0.25) [Z1 Z2]
+    0.25 * (Z(0) @ Z(1)) + 0.25 * Z(0) + 0.25 * Z(1) + 0.25 * (Z(1) @ Z(2)) + 0.25 * Z(1) + 0.25 * Z(2)
 
     >>> import rustworkx as rx
     >>> graph = rx.PyGraph()
@@ -109,12 +102,7 @@ def edge_driver(graph: Union[nx.Graph, rx.PyGraph], reward: list):
     >>> graph.add_edges_from([(0, 1,""), (1,2,"")])
     >>> hamiltonian = qaoa.edge_driver(graph, ["11", "10", "01"])
     >>> print(hamiltonian)
-      (0.25) [Z0]
-    + (0.25) [Z1]
-    + (0.25) [Z1]
-    + (0.25) [Z2]
-    + (0.25) [Z0 Z1]
-    + (0.25) [Z1 Z2]
+    0.25 * (Z(0) @ Z(1)) + 0.25 * Z(0) + 0.25 * Z(1) + 0.25 * (Z(1) @ Z(2)) + 0.25 * Z(1) + 0.25 * Z(2)
 
     In the above example, ``"11"``, ``"10"``, and ``"01"`` are assigned a lower
     energy than ``"00"``. For example, a quick calculation of expectation values gives us:
@@ -275,13 +263,9 @@ def maxcut(graph: Union[nx.Graph, rx.PyGraph]):
     >>> graph = nx.Graph([(0, 1), (1, 2)])
     >>> cost_h, mixer_h = qml.qaoa.maxcut(graph)
     >>> print(cost_h)
-      (-1.0) [I0]
-    + (0.5) [Z0 Z1]
-    + (0.5) [Z1 Z2]
+    0.5 * (Z(0) @ Z(1)) + 0.5 * (Z(1) @ Z(2)) + -0.5 * (I(0) @ I(1)) + -0.5 * (I(1) @ I(2))
     >>> print(mixer_h)
-      (1) [X0]
-    + (1) [X1]
-    + (1) [X2]
+    1 * X(0) + 1 * X(1) + 1 * X(2)
 
     >>> import rustworkx as rx
     >>> graph = rx.PyGraph()
@@ -289,13 +273,9 @@ def maxcut(graph: Union[nx.Graph, rx.PyGraph]):
     >>> graph.add_edges_from([(0, 1,""), (1,2,"")])
     >>> cost_h, mixer_h = qml.qaoa.maxcut(graph)
     >>> print(cost_h)
-      (-1.0) [I0]
-    + (0.5) [Z0 Z1]
-    + (0.5) [Z1 Z2]
+    0.5 * (Z(0) @ Z(1)) + 0.5 * (Z(1) @ Z(2)) + -0.5 * (I(0) @ I(1)) + -0.5 * (I(1) @ I(2))
     >>> print(mixer_h)
-      (1) [X0]
-    + (1) [X1]
-    + (1) [X2]
+    1 * X(0) + 1 * X(1) + 1 * X(2)
     """
 
     if not isinstance(graph, (nx.Graph, rx.PyGraph)):
diff --git a/pennylane/qaoa/cycle.py b/pennylane/qaoa/cycle.py
index 5d6a3de3cad..0c46e5eb6bd 100644
--- a/pennylane/qaoa/cycle.py
+++ b/pennylane/qaoa/cycle.py
@@ -342,13 +342,15 @@ def loss_hamiltonian(graph: Union[nx.Graph, rx.PyGraph, rx.PyDiGraph]) -> qml.op
     >>> for k, v in edge_weight_data.items():
             g[k[0]][k[1]]["weight"] = v
     >>> h = loss_hamiltonian(g)
-    >>> print(h)
-      (-0.6931471805599453) [Z0]
-    + (0.0) [Z1]
-    + (0.4054651081081644) [Z2]
-    + (0.6931471805599453) [Z3]
-    + (0.9162907318741551) [Z4]
-    + (1.0986122886681098) [Z5]
+    >>> h
+    (
+        -0.6931471805599453 * Z(0)
+      + 0.0 * Z(1)
+      + 0.4054651081081644 * Z(2)
+      + 0.6931471805599453 * Z(3)
+      + 0.9162907318741551 * Z(4)
+      + 1.0986122886681098 * Z(5)
+    )
 
     >>> import rustworkx as rx
     >>> g = rx.generators.directed_mesh_graph(3)
@@ -357,12 +359,14 @@ def loss_hamiltonian(graph: Union[nx.Graph, rx.PyGraph, rx.PyDiGraph]) -> qml.op
             g.update_edge(k[0], k[1], {"weight": v})
     >>> h = loss_hamiltonian(g)
     >>> print(h)
-      (-0.6931471805599453) [Z0]
-    + (0.0) [Z1]
-    + (0.4054651081081644) [Z2]
-    + (0.6931471805599453) [Z3]
-    + (0.9162907318741551) [Z4]
-    + (1.0986122886681098) [Z5]
+    (
+        -0.6931471805599453 * Z(0)
+      + 0.0 * Z(1)
+      + 0.4054651081081644 * Z(2)
+      + 0.6931471805599453 * Z(3)
+      + 0.9162907318741551 * Z(4)
+      + 1.0986122886681098 * Z(5)
+    )
 
     Args:
         graph (nx.Graph or rx.PyGraph or rx.PyDiGraph): the graph specifying possible edges
diff --git a/pennylane/qaoa/mixers.py b/pennylane/qaoa/mixers.py
index 881a4dcb7f2..565f60001b9 100644
--- a/pennylane/qaoa/mixers.py
+++ b/pennylane/qaoa/mixers.py
@@ -53,9 +53,7 @@ def x_mixer(wires: Union[Iterable, Wires]):
     >>> wires = range(3)
     >>> mixer_h = qaoa.x_mixer(wires)
     >>> print(mixer_h)
-      (1) [X0]
-    + (1) [X1]
-    + (1) [X2]
+    1 * X(0) + 1 * X(1) + 1 * X(2)
     """
 
     wires = Wires(wires)
@@ -172,15 +170,17 @@ def bit_flip_mixer(graph: Union[nx.Graph, rx.PyGraph], b: int):
     >>> from networkx import Graph
     >>> graph = Graph([(0, 1), (1, 2)])
     >>> mixer_h = qaoa.bit_flip_mixer(graph, 0)
-    >>> print(mixer_h)
-      (0.25) [X1]
-    + (0.5) [X0]
-    + (0.5) [X2]
-    + (0.25) [X1 Z2]
-    + (0.25) [X1 Z0]
-    + (0.5) [X0 Z1]
-    + (0.5) [X2 Z1]
-    + (0.25) [X1 Z0 Z2]
+    >>> mixer_h
+    (
+        0.5 * X(0)
+      + 0.5 * (X(0) @ Z(1))
+      + 0.25 * X(1)
+      + 0.25 * (X(1) @ Z(2))
+      + 0.25 * (X(1) @ Z(0))
+      + 0.25 * (X(1) @ Z(0) @ Z(2))
+      + 0.5 * X(2)
+      + 0.5 * (X(2) @ Z(1))
+    )
 
     >>> import rustworkx as rx
     >>> graph = rx.PyGraph()
@@ -188,14 +188,16 @@ def bit_flip_mixer(graph: Union[nx.Graph, rx.PyGraph], b: int):
     >>> graph.add_edges_from([(0, 1, ""), (1, 2, "")])
     >>> mixer_h = qaoa.bit_flip_mixer(graph, 0)
     >>> print(mixer_h)
-      (0.25) [X1]
-    + (0.5) [X0]
-    + (0.5) [X2]
-    + (0.25) [X1 Z0]
-    + (0.25) [X1 Z2]
-    + (0.5) [X0 Z1]
-    + (0.5) [X2 Z1]
-    + (0.25) [X1 Z2 Z0]
+    (
+        0.5 * X(0)
+      + 0.5 * (X(0) @ Z(1))
+      + 0.25 * X(1)
+      + 0.25 * (X(1) @ Z(2))
+      + 0.25 * (X(1) @ Z(0))
+      + 0.25 * (X(1) @ Z(0) @ Z(2))
+      + 0.5 * X(2)
+      + 0.5 * (X(2) @ Z(1))
+    )
     """
 
     if not isinstance(graph, (nx.Graph, rx.PyGraph)):
diff --git a/pennylane/qchem/convert.py b/pennylane/qchem/convert.py
index bebf1e083b5..7794c2b54f2 100644
--- a/pennylane/qchem/convert.py
+++ b/pennylane/qchem/convert.py
@@ -137,7 +137,7 @@ def _openfermion_to_pennylane(qubit_operator, wires=None):
 
     **Example**
 
-    >>> q_op = 0.1*QubitOperator('X0') + 0.2*QubitOperator('Y0 Z2')
+    >>> q_op = 0.1 * QubitOperator('X0') + 0.2 * QubitOperator('Y0 Z2')
     >>> q_op
     0.1 [X0] +
     0.2 [Y0 Z2]
diff --git a/pennylane/qchem/number.py b/pennylane/qchem/number.py
index d31a1d51429..df75ba01753 100644
--- a/pennylane/qchem/number.py
+++ b/pennylane/qchem/number.py
@@ -47,11 +47,13 @@ def particle_number(orbitals):
 
     >>> orbitals = 4
     >>> print(particle_number(orbitals))
-    (2.0) [I0]
-    + (-0.5) [Z0]
-    + (-0.5) [Z1]
-    + (-0.5) [Z2]
-    + (-0.5) [Z3]
+    (
+        2.0 * I(0)
+      + -0.5 * Z(0)
+      + -0.5 * Z(1)
+      + -0.5 * Z(2)
+      + -0.5 * Z(3)
+    )
     """
 
     if orbitals <= 0:
diff --git a/pennylane/qchem/openfermion_obs.py b/pennylane/qchem/openfermion_obs.py
index d7fba9e7f1f..e1716b72f9a 100644
--- a/pennylane/qchem/openfermion_obs.py
+++ b/pennylane/qchem/openfermion_obs.py
@@ -124,13 +124,15 @@ def observable(fermion_ops, init_term=0, mapping="jordan_wigner", wires=None):
 
     >>> t = FermionOperator("0^ 0", 0.5) + FermionOperator("1^ 1", 0.25)
     >>> v = FermionOperator("1^ 0^ 0 1", -0.15) + FermionOperator("2^ 0^ 2 0", 0.3)
-    >>> print(observable([t, v], mapping="jordan_wigner"))
-    (0.2625) [I0]
-    + (-0.1375) [Z0]
-    + (-0.0875) [Z1]
-    + (-0.0375) [Z0 Z1]
-    + (0.075) [Z2]
-    + (-0.075) [Z0 Z2]
+    >>> observable([t, v], mapping="jordan_wigner")
+    (
+        0.2625 * I(0)
+      + -0.1375 * Z(0)
+      + -0.0875 * Z(1)
+      + -0.0375 * (Z(0) @ Z(1))
+      + 0.075 * Z(2)
+      + -0.075 * (Z(0) @ Z(2))
+    )
     """
     openfermion, _ = _import_of()
 
@@ -563,30 +565,30 @@ def dipole_of(
     >>> symbols = ["H", "H", "H"]
     >>> coordinates = np.array([0.028, 0.054, 0.0, 0.986, 1.610, 0.0, 1.855, 0.002, 0.0])
     >>> dipole_obs = dipole_of(symbols, coordinates, charge=1)
-    >>> print(dipole_obs)
-    [<Hamiltonian: terms=18, wires=[0, 1, 2, 3, 4, 5]>,
-    <Hamiltonian: terms=18, wires=[0, 1, 2, 3, 4, 5]>,
-    <Hamiltonian: terms=1, wires=[0]>]
-
-    >>> print(dipole_obs[0]) # x-component of D
-    (0.24190977644628117) [Z4]
-    + (0.24190977644628117) [Z5]
-    + (0.4781123173263878) [Z0]
-    + (0.4781123173263878) [Z1]
-    + (0.714477906181248) [Z2]
-    + (0.714477906181248) [Z3]
-    + (-0.3913638489487808) [Y0 Z1 Y2]
-    + (-0.3913638489487808) [X0 Z1 X2]
-    + (-0.3913638489487808) [Y1 Z2 Y3]
-    + (-0.3913638489487808) [X1 Z2 X3]
-    + (-0.1173495878099553) [Y2 Z3 Y4]
-    + (-0.1173495878099553) [X2 Z3 X4]
-    + (-0.1173495878099553) [Y3 Z4 Y5]
-    + (-0.1173495878099553) [X3 Z4 X5]
-    + (0.26611147045300276) [Y0 Z1 Z2 Z3 Y4]
-    + (0.26611147045300276) [X0 Z1 Z2 Z3 X4]
-    + (0.26611147045300276) [Y1 Z2 Z3 Z4 Y5]
-    + (0.26611147045300276) [X1 Z2 Z3 Z4 X5]
+    >>> print([(h.wires) for h in dipole_obs])
+    [<Wires = [0, 1, 2, 3, 4, 5]>, <Wires = [0, 1, 2, 3, 4, 5]>, <Wires = [0]>]
+
+    >>> dipole_obs[0] # x-component of D
+    (
+        0.4781123173263876 * Z(0)
+      + 0.4781123173263876 * Z(1)
+      + -0.3913638489489803 * (Y(0) @ Z(1) @ Y(2))
+      + -0.3913638489489803 * (X(0) @ Z(1) @ X(2))
+      + -0.3913638489489803 * (Y(1) @ Z(2) @ Y(3))
+      + -0.3913638489489803 * (X(1) @ Z(2) @ X(3))
+      + 0.2661114704527088 * (Y(0) @ Z(1) @ Z(2) @ Z(3) @ Y(4))
+      + 0.2661114704527088 * (X(0) @ Z(1) @ Z(2) @ Z(3) @ X(4))
+      + 0.2661114704527088 * (Y(1) @ Z(2) @ Z(3) @ Z(4) @ Y(5))
+      + 0.2661114704527088 * (X(1) @ Z(2) @ Z(3) @ Z(4) @ X(5))
+      + 0.7144779061810713 * Z(2)
+      + 0.7144779061810713 * Z(3)
+      + -0.11734958781031017 * (Y(2) @ Z(3) @ Y(4))
+      + -0.11734958781031017 * (X(2) @ Z(3) @ X(4))
+      + -0.11734958781031017 * (Y(3) @ Z(4) @ Y(5))
+      + -0.11734958781031017 * (X(3) @ Z(4) @ X(5))
+      + 0.24190977644645698 * Z(4)
+      + 0.24190977644645698 * Z(5)
+    )
     """
     openfermion, _ = _import_of()
 
diff --git a/pennylane/qchem/spin.py b/pennylane/qchem/spin.py
index a7ada7ae89c..01a67de2320 100644
--- a/pennylane/qchem/spin.py
+++ b/pennylane/qchem/spin.py
@@ -140,26 +140,28 @@ def spin2(electrons, orbitals):
 
     >>> electrons = 2
     >>> orbitals = 4
-    >>> print(spin2(electrons, orbitals))
-    (0.75) [I0]
-    + (0.375) [Z1]
-    + (-0.375) [Z0 Z1]
-    + (0.125) [Z0 Z2]
-    + (0.375) [Z0]
-    + (-0.125) [Z0 Z3]
-    + (-0.125) [Z1 Z2]
-    + (0.125) [Z1 Z3]
-    + (0.375) [Z2]
-    + (0.375) [Z3]
-    + (-0.375) [Z2 Z3]
-    + (0.125) [Y0 X1 Y2 X3]
-    + (0.125) [Y0 Y1 X2 X3]
-    + (0.125) [Y0 Y1 Y2 Y3]
-    + (-0.125) [Y0 X1 X2 Y3]
-    + (-0.125) [X0 Y1 Y2 X3]
-    + (0.125) [X0 X1 X2 X3]
-    + (0.125) [X0 X1 Y2 Y3]
-    + (0.125) [X0 Y1 X2 Y3]
+    >>> spin2(electrons, orbitals)
+    (
+        0.75 * I(0)
+      + 0.375 * Z(0)
+      + 0.375 * Z(1)
+      + -0.375 * (Z(0) @ Z(1))
+      + 0.375 * Z(2)
+      + 0.125 * (Z(0) @ Z(2))
+      + 0.375 * Z(3)
+      + -0.125 * (Z(0) @ Z(3))
+      + -0.125 * (Z(1) @ Z(2))
+      + 0.125 * (Z(1) @ Z(3))
+      + -0.375 * (Z(2) @ Z(3))
+      + 0.125 * (Y(0) @ Y(2) @ X(3) @ X(1))
+      + 0.125 * (Y(0) @ X(2) @ X(3) @ Y(1))
+      + 0.125 * (Y(0) @ Y(2) @ Y(3) @ Y(1))
+      + -0.125 * (Y(0) @ X(2) @ Y(3) @ X(1))
+      + -0.125 * (X(0) @ Y(2) @ X(3) @ Y(1))
+      + 0.125 * (X(0) @ X(2) @ X(3) @ X(1))
+      + 0.125 * (X(0) @ Y(2) @ Y(3) @ X(1))
+      + 0.125 * (X(0) @ X(2) @ Y(3) @ Y(1))
+    )
     """
 
     if electrons <= 0:
@@ -221,10 +223,12 @@ def spinz(orbitals):
 
     >>> orbitals = 4
     >>> print(spinz(orbitals))
-    (-0.25) [Z0]
-    + (0.25) [Z1]
-    + (-0.25) [Z2]
-    + (0.25) [Z3]
+    (
+        -0.25 * Z(0)
+      + 0.25 * Z(1)
+      + -0.25 * Z(2)
+      + 0.25 * Z(3)
+    )
     """
 
     if orbitals <= 0:
diff --git a/pennylane/qchem/tapering.py b/pennylane/qchem/tapering.py
index 0f3f94692d3..8dccd068f0a 100644
--- a/pennylane/qchem/tapering.py
+++ b/pennylane/qchem/tapering.py
@@ -143,11 +143,7 @@ def symmetry_generators(h):
     >>> H, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)
     >>> t = symmetry_generators(H)
     >>> t
-    [<Hamiltonian: terms=1, wires=[0, 1]>,
-     <Hamiltonian: terms=1, wires=[0, 2]>,
-     <Hamiltonian: terms=1, wires=[0, 3]>]
-    >>> print(t[0])
-    (1.0) [Z0 Z1]
+    [Z(0) @ Z(1), Z(0) @ Z(2), Z(0) @ Z(3)]
     """
     num_qubits = len(h.wires)
 
@@ -372,10 +368,12 @@ def taper(h, generators, paulixops, paulix_sector):
     >>> paulixops = paulix_ops(generators, 4)
     >>> paulix_sector = [1, -1, -1]
     >>> H_tapered = taper(H, generators, paulixops, paulix_sector)
-    >>> print(H_tapered)
-      ((-0.321034397355757+0j)) [I0]
-    + ((0.1809270275619003+0j)) [X0]
-    + ((0.7959678503869626+0j)) [Z0]
+    >>> H_tapered
+    (
+        (-0.3210343973331179-2.0816681711721685e-17j) * I(0)
+      + (0.7959678504583807+0j) * Z(0)
+      + (0.18092702760702645+0j) * X(0)
+    )
     """
 
     ps_h = pauli_sentence(h)
diff --git a/pennylane/workflow/construct_batch.py b/pennylane/workflow/construct_batch.py
index 7154a754839..79f3de45f65 100644
--- a/pennylane/workflow/construct_batch.py
+++ b/pennylane/workflow/construct_batch.py
@@ -229,9 +229,8 @@ def circuit(x):
         >>> batch, fn = construct_batch(circuit, level="device")(1.23)
         >>> batch[0].circuit
         [RY(tensor(1., requires_grad=True), wires=[1]),
-        RX(tensor(2., requires_grad=True), wires=[0]),
-        expval(  (1) [X0]
-        + (1) [Y0])]
+         RX(tensor(2., requires_grad=True), wires=[0]),
+         expval(X(0) + Y(0))]
 
         These tapes can be natively executed by the device, though with non-backprop devices the parameters
         will need to be converted to numpy with :func:`~.convert_to_numpy_parameters`.
@@ -244,22 +243,20 @@ def circuit(x):
         >>> batch, fn = construct_batch(circuit, level="gradient")(1.23)
         >>> batch[0].circuit
         [RY(tensor(1., requires_grad=True), wires=[1]),
-        RX(tensor(2., requires_grad=True), wires=[0]),
-        expval(  (1) [X0]
-        + (1) [Y0])]
+         RX(tensor(2., requires_grad=True), wires=[0]),
+         expval(X(0) + Y(0))]
 
         We can inspect what was directly captured from the qfunc with ``level=0``.
 
         >>> batch, fn = construct_batch(circuit, level=0)(1.23)
         >>> batch[0].circuit
         [RandomLayers(tensor([[1., 2.]], requires_grad=True), wires=[0, 1]),
-        RX(1.23, wires=[0]),
-        RX(-1.23, wires=[0]),
-        SWAP(wires=[0, 1]),
-        X(0),
-        X(0),
-        expval(  (1) [X0]
-        + (1) [Y0])]
+         RX(1.23, wires=[0]),
+         RX(-1.23, wires=[0]),
+         SWAP(wires=[0, 1]),
+         X(0),
+         X(0),
+         expval(X(0) + Y(0))]
 
         And iterate though stages in the transform program with different integers.
         If we request ``level=1``, the ``cancel_inverses`` transform has been applied.
@@ -267,11 +264,10 @@ def circuit(x):
         >>> batch, fn = construct_batch(circuit, level=1)(1.23)
         >>> batch[0].circuit
         [RandomLayers(tensor([[1., 2.]], requires_grad=True), wires=[0, 1]),
-        RX(1.23, wires=[0]),
-        RX(-1.23, wires=[0]),
-        SWAP(wires=[0, 1]),
-        expval(  (1) [X0]
-        + (1) [Y0])]
+         RX(1.23, wires=[0]),
+         RX(-1.23, wires=[0]),
+         SWAP(wires=[0, 1]),
+         expval(X(0) + Y(0))]
 
         We can also slice into a subset of the transform program.  ``slice(1, None)`` would skip the first user
         transform ``cancel_inverses``:
@@ -279,11 +275,10 @@ def circuit(x):
         >>> batch, fn = construct_batch(circuit, level=slice(1,None))(1.23)
         >>> batch[0].circuit
         [RY(tensor(1., requires_grad=True), wires=[1]),
-        RX(tensor(2., requires_grad=True), wires=[0]),
-        X(0),
-        X(0),
-        expval(  (1) [X0]
-        + (1) [Y0])]
+         RX(tensor(2., requires_grad=True), wires=[0]),
+         X(0),
+         X(0),
+         expval(X(0) + Y(0))]
 
     """