Retropropagación del operador (OBP) para la estimación de valores esperados

Tiempo estimado de ejecución: 4 minutos en un procesador Heron r3 (NOTA: Se trata únicamente de una estimación). (El tiempo de ejecución puede variar.)

Resultados del aprendizaje

Una vez completado este tutorial, los usuarios deberían comprender:

• Cómo utilizar qiskit-addon-obp para reducir la profundidad del circuito cuántico a costa de aumentar el número de ejecuciones del circuito
• Cómo utilizar qiskit-addon-utils para construir hamiltonianos XYZ y sus circuitos de evolución temporal

Requisitos previos

Recomendamos a los usuarios que se familiaricen con los siguientes temas antes de seguir este tutorial:

• Uso de la primitiva «Estimator» para calcular los valores esperados de una variable observable

En segundo plano

La retropropagación de operadores es una técnica que consiste en incorporar operaciones del final de un circuito cuántico en el observable medido, lo que suele reducir la profundidad del circuito a costa de introducir términos adicionales en el observable. El objetivo es realizar la retropropagación en la mayor parte posible del circuito sin que la magnitud observable aumente demasiado. En el complemento OBP Qiskit hay disponible una implementación basada en Qiskit. Consulte la documentación correspondiente para obtener más información.

Consideremos un circuito de ejemplo en el que se va a medir una observable $O = \sum_P c_P P$, donde $P$ son los estados de Pauli y $c_P$ son los coeficientes. Denotemos el circuito como un único « $U$ » unitario, que puede dividirse lógicamente en un « $U = U_C U_Q$ », tal y como se muestra en la siguiente figura.

La retropropagación de operadores absorbe el operador unitario $U_C$ en el observable, haciéndolo evolucionar como $O' = U_C^{\dagger}OU_C = \sum_P c_P U_C^{\dagger}PU_C$. En otras palabras, parte del cálculo se realiza de forma clásica mediante la evolución del observable desde $O$ hasta $O'$. El problema original puede reformularse ahora como la medición del observable $O'$ para el nuevo circuito de menor profundidad cuyo operador unitario es $U_Q$.

El espacio unitario $U_C$ se representa como un conjunto de subespacios $U_C = U_S U_{S-1}...U_2U_1$. Existen varias formas de definir un subespacio. Por ejemplo, en el circuito del ejemplo anterior, cada capa de puertas « $R_{zz}$ » y cada capa de puertas « $R_x$ » pueden considerarse como una sección individual. La retropropagación implica, tradicionalmente, el cálculo de un $O' = \Pi_{s=1}^S \sum_P c_P U_s^{\dagger} P U_s$. Cada corte $U_s$ puede representarse como $U_s = exp(\frac{-i\theta_s P_s}{2})$, donde $P_s$ es un operador de Pauli de -qubits $n$ y $\theta_s$ es un escalar. Es fácil comprobar que

En el ejemplo anterior, si se cumple que « $\{P,P_s\} = 0$ », entonces tendremos que ejecutar dos circuitos cuánticos, en lugar de uno, para calcular el valor esperado. Por lo tanto, la retropropagación podría aumentar el número de términos en el observable, lo que daría lugar a un mayor número de ejecuciones del circuito. Una forma de permitir una retropropagación más profunda en el circuito, evitando al mismo tiempo que el operador se haga demasiado grande, consiste en truncar los términos con coeficientes pequeños, en lugar de añadirlos al operador. Por ejemplo, en el ejemplo anterior, se podría optar por truncar el término que incluye « $P_sP$ », siempre que « $\theta_s$ » sea lo suficientemente pequeño. El truncamiento de términos puede reducir el número de circuitos cuánticos que hay que ejecutar, pero ello da lugar a un error en el cálculo del valor esperado final que es proporcional a la magnitud de los coeficientes de los términos truncados.

Requisitos

Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:

• Qiskit SDK v2.0 o posterior, con soporte para visualización
• Qiskit Runtime v0.22 o posterior (pip install qiskit-ibm-runtime)
• Complemento OBP Qiskit 0.3 o posterior (pip install qiskit-addon-obp)
• Utilidades del complemento Qiskit: 0.3 o posterior (pip install qiskit-addon-utils)

Configuración

import numpy as np
import matplotlib.pyplot as plt

from qiskit.primitives import StatevectorEstimator
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.quantum_info import SparsePauliOp
from qiskit.transpiler import CouplingMap
from qiskit.synthesis import LieTrotter

from qiskit_addon_utils.problem_generators import generate_xyz_hamiltonian
from qiskit_addon_utils.problem_generators import (
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth, combine_slices
from qiskit_addon_obp.utils.simplify import OperatorBudget
from qiskit_addon_obp import backpropagate
from qiskit_addon_obp.utils.truncating import setup_budget

from rustworkx.visualization import graphviz_draw

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import EstimatorV2, EstimatorOptions

Ejemplo de simulador a pequeña escala

Este tutorial implementa un patrón de Qiskit para simular la dinámica cuántica de una cadena de espines de Heisenberg utilizando el complemento OBP de Qiskit. Ten en cuenta que, en un simulador sin ruido, el valor esperado obtenido con y sin retropropagación será el mismo.

Paso 1: Asignar entradas clásicas a un problema cuántico

Representar la evolución temporal de un modelo cuántico de Heisenberg en un experimento cuántico

En primer lugar, utilizaremos la generate_xyz_hamiltonian función de qiskit-addon-utils para generar un hamiltoniano de tipo Heisenberg en un grafo de conectividad dado. Este gráfico puede ser tanto un gráfico de « rustworkx.PyGraph » como uno de « CouplingMap ». A continuación, utilizaremos una cadena CouplingMap lineal de 10 qubits.

num_qubits = 10
layout = [(i - 1, i) for i in range(1, num_qubits)]

# Instantiate a CouplingMap object
coupling_map = CouplingMap(layout)
graphviz_draw(coupling_map.graph, method="circo")

Output:

A continuación, generamos un operador de Pauli que modela un hamiltoniano XYZ de Heisenberg:

donde $G(V,E)$ es la gráfica de la aplicación de acoplamiento. En este tutorial, hemos utilizado $J_x, J_y, J_z$ y $\frac{\pi}{8}, \frac{\pi}{4}, \frac{\pi}{2}$, respectivamente, y $h_x, h_y, h_z$ y $\frac{\pi}{3}, \frac{\pi}{6}, \frac{\pi}{9}$, respectivamente.

# Get a qubit operator describing the Heisenberg XYZ model
hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(np.pi / 8, np.pi / 4, np.pi / 2),
    ext_magnetic_field=(np.pi / 3, np.pi / 6, np.pi / 9),
)
print(hamiltonian)

Output:

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII', 'IIIIIIIIIX', 'IIIIIIIIIY', 'IIIIIIIIIZ', 'IIIIIIIIXI', 'IIIIIIIIYI', 'IIIIIIIIZI', 'IIIIIIIXII', 'IIIIIIIYII', 'IIIIIIIZII', 'IIIIIIXIII', 'IIIIIIYIII', 'IIIIIIZIII', 'IIIIIXIIII', 'IIIIIYIIII', 'IIIIIZIIII', 'IIIIXIIIII', 'IIIIYIIIII', 'IIIIZIIIII', 'IIIXIIIIII', 'IIIYIIIIII', 'IIIZIIIIII', 'IIXIIIIIII', 'IIYIIIIIII', 'IIZIIIIIII', 'IXIIIIIIII', 'IYIIIIIIII', 'IZIIIIIIII', 'XIIIIIIIII', 'YIIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j,
 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j,
 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j,
 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j,
 0.78539816+0.j, 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j,
 1.57079633+0.j, 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j,
 0.39269908+0.j, 0.78539816+0.j, 1.57079633+0.j, 1.04719755+0.j,
 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j,
 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j,
 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j,
 1.04719755+0.j, 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j,
 0.52359878+0.j, 0.34906585+0.j, 1.04719755+0.j, 0.52359878+0.j,
 0.34906585+0.j])

A partir del operador del qubit, podemos generar un circuito cuántico que modele su evolución temporal. Hemos utilizado generate_time_evolution_circuit la descomposición de Lie-Trotter para construir el circuito de evolución temporal.

circuit = generate_time_evolution_circuit(
    hamiltonian,
    time=0.2,
    synthesis=LieTrotter(reps=2),
)
circuit.draw("mpl", style="iqp", fold=-1)

Output:

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

Crear segmentos de circuito para retropropagación

La backpropagate función realiza la retropropagación de secciones completas del circuito de una vez. Por lo tanto, la elección del método de segmentación puede influir en el rendimiento de la retropropagación para un problema concreto. Aquí agruparemos las puertas del mismo tipo en secciones utilizando la slice_by_depth función.

Si quieres obtener información más detallada sobre el «circuit slicing», consulta esta guía práctica del paquete qiskit-addon-utils .

slices = slice_by_depth(circuit, max_slice_depth=1)
print(f"Separated the circuit into {len(slices)} slices.")

Output:

Separated the circuit into 18 slices.

Limitar el tamaño que puede alcanzar el operador durante la retropropagación

Durante la retropropagación, el número de términos del operador tenderá, por lo general, a acercarse rápidamente a $2^L$, donde $L$ es el número de cortes. Cuando dos términos del operador no conmutan a nivel de qubit, necesitamos circuitos independientes para obtener los valores esperados que les corresponden. Por ejemplo, si tenemos un observable de dos qubits $O = 0.1 XX + 0.3 IZ - 0.5 IX$, dado que $[XX,IX] = 0$, basta con realizar una medición en una sola base para calcular los valores esperados de estos dos términos. Sin embargo, $IZ$ es anticomutativo con los otros dos términos, por lo que necesitamos una medición de base independiente para calcular el valor esperado de $IZ$. En otras palabras, necesitamos dos circuitos en lugar de uno para calcular $\langle O \rangle$. A medida que aumenta el número de términos del operador, existe la posibilidad de que también aumente el número de ejecuciones de circuitos necesarias.

El tamaño del operador puede limitarse especificando el kwarg operator_budget de la función backpropagate , que acepta una OperatorBudget instancia.

Para controlar la cantidad de recursos adicionales (número de ejecuciones de circuitos y, por lo tanto, el tiempo de QPU necesario) que se asignan, limitamos el número máximo de grupos de Pauli conmutativos a nivel de qubit que puede tener el observable retropropagado. Aquí especificamos que la retropropagación debe detenerse cuando el número de grupos de Pauli conmutativos por qubit en el operador supere los ocho.

op_budget = OperatorBudget(max_qwc_groups=8)

Retropropagación de segmentos del circuito

Primero especificamos que el observable es $M_Z = \frac{1}{N} \sum_{i=1}^N \langle Z_i \rangle$, siendo $N$ el número de qubits. Retropropagaremos cortes del circuito de evolución temporal hasta que los términos del observable ya no puedan combinarse en ocho o menos grupos de Pauli conmutables en qubit.

observable = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
    num_qubits=num_qubits,
)
observable

Output:

SparsePauliOp(['IIIIIIIIIZ', 'IIIIIIIIZI', 'IIIIIIIZII', 'IIIIIIZIII', 'IIIIIZIIII', 'IIIIZIIIII', 'IIIZIIIIII', 'IIZIIIIIII', 'IZIIIIIIII', 'ZIIIIIIIII'],
              coeffs=[0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j, 0.1+0.j,
 0.1+0.j, 0.1+0.j])

A continuación verás que hemos aplicado la retropropagación a seis cortes, y que los términos se han agrupado en seis grupos, en lugar de ocho. Esto implica que, si se aplicara la retropropagación a una capa más, el número de grupos de Pauli superaría los ocho. Podemos comprobar que así es examinando los metadatos devueltos. Ten en cuenta también que, en esta parte, la transformación del circuito es exacta. Es decir, no se ha truncado ningún término de la nueva observable $O’$. El circuito de retropropagación y el operador de retropropagación dan el mismo resultado que el circuito y el operador originales.

# Backpropagate slices onto the observable
bp_obs, remaining_slices, metadata = backpropagate(
    observable, slices, operator_budget=op_budget
)
# Recombine the slices remaining after backpropagation
bp_circuit = combine_slices(remaining_slices)

print(f"Backpropagated {metadata.num_backpropagated_slices} slices.")
print(
    f"New observable has {len(bp_obs.paulis)} terms, which can be combined into "
    f"{len(bp_obs.group_commuting(qubit_wise=True))} groups."
)
print(
    f"Note that backpropagating one more slice would result in "
    f"{metadata.backpropagation_history[-1].num_paulis[0]} terms "
    f"across {metadata.backpropagation_history[-1].num_qwc_groups} groups."
)
print("The remaining circuit after backpropagation looks as follows:")
bp_circuit.draw("mpl", fold=-1, scale=0.6)

Output:

Backpropagated 6 slices.
New observable has 60 terms, which can be combined into 6 groups.
Note that backpropagating one more slice would result in 114 terms across 12 groups.
The remaining circuit after backpropagation looks as follows:

Para el ejemplo a pequeña escala en un simulador, no utilizaremos el truncamiento. Esto se debe a que, en ausencia de ruido, el circuito con y sin retropropagación da el mismo resultado, y el truncamiento empeora el resultado debido a la aproximación añadida.

Transpilar los circuitos al conjunto de puertas básico

Ahora compilamos tanto el circuito original como el resultante de la retropropagación en la puerta base del backend. No es necesario realizar la transpilación en el backend real, ya que vamos a ejecutar el código en un simulador para la instancia pequeña.

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=133
)
print(backend)

Output:

<IBMBackend('ibm_kingston')>

pm_basis = generate_preset_pass_manager(
    optimization_level=3, basis_gates=backend.configuration().basis_gates
)
isa_circuit = pm_basis.run(circuit)
isa_bp_circuit = pm_basis.run(bp_circuit)

Paso 3: Ejecutar utilizando Qiskit primitives

En primer lugar, creamos dos bloques primitivos unificados (PUB) que corresponden al circuito original y al circuito retropropagado. A continuación, aplicamos los intervalos de confianza a un estimador ideal para obtener los valores esperados.

pubs = [(isa_circuit, observable), (isa_bp_circuit, bp_obs)]

rng = np.random.default_rng()
estimator = StatevectorEstimator(seed=rng)
job = estimator.run(pubs)

Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado

Ahora obtenemos los valores esperados de los circuitos originales y de los circuitos con retropropagación.

primitive_result = job.result()
circuit_expval = primitive_result[0].data.evs.item()
bp_circuit_expval = primitive_result[1].data.evs.item()

methods = [
    "No backpropagation",
    "Backpropagation",
]
values = [circuit_expval, bp_circuit_expval]

ax = plt.gca()
plt.bar(methods, values, color="#a56eff", width=0.4, edgecolor="#8a3ffc")
ax.set_ylim([0.6, 0.92])
ax.set_ylabel(r"$M_Z$", fontsize=12)

Output:

Text(0, 0.5, '$M_Z$')

Como era de esperar, los dos valores esperados coinciden. Dado que estamos trabajando con un simulador de vectores de estado sin ruido, la retropropagación es una transformación exacta del par circuito-observable, por lo que los flujos de trabajo original y el de retropropagación deben producir el mismo valor de « $M_Z$ ». La ventaja de la retropropagación solo se hace evidente en hardware con ruido, donde el circuito de retropropagación, al ser más corto, acumula menos error, tal y como se ilustra en el ejemplo de hardware a gran escala que se muestra a continuación.

Ejemplo de hardware a gran escala

A la hora de diseñar un experimento, resulta útil empezar con un circuito sencillo para facilitar las visualizaciones y las simulaciones. Ahora analizamos la retropropagación de operadores para un hamiltoniano de Heisenberg de 50 qubits con el mismo conjunto de valores para los parámetros $J$ y $h$ y el mismo observable $M_Z$, pero para cuatro pasos de Trotter. El valor esperado ideal a esta escala no se puede calcular mediante un método de fuerza bruta, por lo que utilizamos una red tensorial y obtenemos que el valor esperado ideal es $\simeq 0.89$

Además de la retropropagación, en este ejemplo a gran escala también presentamos la retropropagación con truncamiento. Lo ideal es realizar la retropropagación lo máximo posible para reducir la profundidad del circuito efectivo. Sin embargo, esto suele dar lugar a un gran número de términos no conmutativos en el observable actualizado, lo que aumenta la sobrecarga cuántica. Por lo tanto, podemos eliminar los términos observables con coeficientes pequeños mediante una técnica denominada truncamiento. Si bien el truncamiento permite una mayor propagación al reducir el número de términos en el observable actualizado, también introduce cierta aproximación. Por lo tanto, es necesario limitar el truncamiento dentro de ciertos límites para que el error de aproximación no supere la reducción del ruido obtenida gracias a una retropropagación más profunda.

Para limitar el grado de truncamiento, asignamos un margen de error a cada segmento, así como un margen de error total para todo el circuito de retropropagación, utilizando la setup_budget función. Esto garantiza que el truncamiento se controle tanto para cada sección como para el circuito completo. Consulte también esta guía para conocer otras formas de distribuir el presupuesto.

num_qubits = 50
layout = [(i - 1, i) for i in range(1, num_qubits)]

# Instantiate a CouplingMap object
coupling_map = CouplingMap(layout)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(np.pi / 8, np.pi / 4, np.pi / 2),
    ext_magnetic_field=(np.pi / 3, np.pi / 6, np.pi / 9),
)

# Generate a time evolution circuit for the Hamiltonian
circuit = generate_time_evolution_circuit(
    hamiltonian,
    time=0.2,
    synthesis=LieTrotter(reps=4),
)

# Define the observable to measure
observable = SparsePauliOp.from_sparse_list(
    [("Z", [i], 1 / num_qubits) for i in range(num_qubits)],
    num_qubits,
)

slices = slice_by_depth(circuit, max_slice_depth=1)

# Define the maximum number of qwc groups allowed in the
# backpropagated observable,
# and the truncation error budget
op_budget = OperatorBudget(max_qwc_groups=15)
truncation_error_budget = setup_budget(
    max_error_total=0.03, max_error_per_slice=0.005
)

# First backpropagation without truncation
bp_obs, remaining_slices, metadata = backpropagate(
    observable, slices, operator_budget=op_budget
)
bp_circuit = combine_slices(remaining_slices)

# Now backpropagate with truncation, using the same operator budget and
# the defined truncation error budget
bp_obs_trunc, remaining_slices_trunc, metadata = backpropagate(
    observable,
    slices,
    operator_budget=op_budget,
    truncation_error_budget=truncation_error_budget,
)
bp_circuit_trunc = combine_slices(
    remaining_slices_trunc, include_barriers=False
)

# Now we transpile the original circuit and the two backpropagated circuits,
# and apply the layout to the corresponding observables
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)

isa_circuit = pm.run(circuit)
isa_bp_circuit = pm.run(bp_circuit)
isa_bp_circuit_trunc = pm.run(bp_circuit_trunc)

isa_observable = observable.apply_layout(isa_circuit.layout)
isa_bp_observable = bp_obs.apply_layout(isa_bp_circuit.layout)
isa_bp_observable_trunc = bp_obs_trunc.apply_layout(
    isa_bp_circuit_trunc.layout
)

# Compare the 2-qubit depth of each transpiled circuit to see how much
# depth backpropagation saved
print(
    f"2-qubit depth without backpropagation: "
    f"{isa_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
    f"2-qubit depth with backpropagation: "
    f"{isa_bp_circuit.depth(lambda x: x.operation.num_qubits == 2)}"
)
print(
    f"2-qubit depth with backpropagation and truncation: "
    f"{isa_bp_circuit_trunc.depth(lambda x: x.operation.num_qubits == 2)}"
)

pubs = [
    (isa_circuit, isa_observable),
    (isa_bp_circuit, isa_bp_observable),
    (isa_bp_circuit_trunc, isa_bp_observable_trunc),
]

# Now we instantiate the Estimator primitive for the hardware with
# ZNE and measurement error
# mitigation and compute the three circuits and observables
options = EstimatorOptions()
options.default_precision = 0.01
options.resilience_level = 2
options.resilience.zne.noise_factors = [1, 1.2, 1.4]
options.resilience.zne.extrapolator = ["linear"]
estimator = EstimatorV2(mode=backend, options=options)

estimator.options.environment.job_tags = ["TUT_OBP"]
job = estimator.run(pubs)

# Retrieve the results and the standard deviations
result_no_bp = job.result()[0].data.evs.item()
result_bp = job.result()[1].data.evs.item()
result_bp_trunc = job.result()[2].data.evs.item()

std_no_bp = job.result()[0].data.stds.item()
std_bp = job.result()[1].data.stds.item()
std_bp_trunc = job.result()[2].data.stds.item()

Output:

2-qubit depth without backpropagation: 24
2-qubit depth with backpropagation: 20
2-qubit depth with backpropagation and truncation: 18

print(f"Expectation value without backpropagation: {result_no_bp}")
print(f"Backpropagated expectation value: {result_bp}")
print(f"Backpropagated expectation value with truncation: {result_bp_trunc}")

Output:

Expectation value without backpropagation: 0.9543907942381811
Backpropagated expectation value: 0.9445337385406468
Backpropagated expectation value with truncation: 0.934050286970965

# Plot the results
methods = [
    "No backpropagation",
    "Backpropagation",
    "Backpropagation w/ truncation",
]
values = [result_no_bp, result_bp, result_bp_trunc]
error_bars = [std_no_bp, std_bp, std_bp_trunc]

ax = plt.gca()
plt.bar(methods, values, color="#a56eff", width=0.4, edgecolor="#8a3ffc")
plt.errorbar(methods, values, yerr=error_bars, fmt="o", color="r", capsize=5)
plt.axhline(0.89)
ax.set_ylim([0.8, 0.98])
plt.text(0.25, 0.895, "Exact result")
ax.set_ylabel(r"$M_Z$", fontsize=12)

Output:

Text(0, 0.5, '$M_Z$')

Próximos pasos

Si te ha parecido interesante este trabajo, quizá te interese el siguiente material: