Dinámica de espín y de carga en el modelo de Fermi-Hubbard

Estimación de uso: 5 minutos en un procesador Heron r3. (NOTA: Es solo una estimación. Tu tiempo de ejecución puede variar.)

Contexto

El modelo de Fermi-Hubbard describe electrones que saltan en una red con una penalización de energía por doble ocupación. En una cadena unidimensional de $L$ sitios con bordes abiertos, el hamiltoniano se define de la siguiente manera:

$$H = -\sum_{i=0}^{L-2} \sum_{\sigma \in \{\uparrow,\downarrow\}} t_{h,i} \left( c^\dagger_{i\sigma} c_{i+1,\sigma} + \mathrm{h.c.} \right) + \sum_{i=0}^{L-1} U_i\, n_{i\uparrow} n_{i\downarrow} - \sum_{i=0}^{L-1} \sum_\sigma \mu_i\, n_{i\sigma},$$

donde $c^\dagger_{i\sigma}$ crea un fermión de espín $\sigma$ en el sitio $i$ y $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ los cuenta. Los tres términos son la amplitud de salto $t_{h,i}$ en el enlace entre los sitios $i$ e $i+1$, la interacción en el sitio $U_i$ y el potencial químico $\mu_i$. El tutorial mide la energía en unidades de $t_h$ y el tiempo en unidades de $1/t_h$.

La dinámica en tiempo real de este modelo es difícil de calcular de forma clásica. La diagonalización exacta se detiene cerca de $L \lesssim 14$ sitios y, después de un quench, la entropía de entrelazamiento de un solucionador de redes tensoriales crece linealmente con el tiempo, por lo que la dimensión de enlace necesaria para una precisión fija crece exponencialmente.

En una dimensión, las excitaciones elementales son los espinones, que portan espín y no carga, y los holones, que portan carga y no espín. Ambos se propagan a velocidades diferentes, por lo que una perturbación localizada se divide en un frente rápido y otro lento. Este efecto se llama separación espín-carga, y este tutorial lo mide [1].

Este tutorial usa la Qiskit Function Dynamics Simulator, que toma el hamiltoniano, un estado inicial, un esquema de Trotter y una lista de observables, y devuelve valores esperados. Se encarga internamente del mapeo de Jordan-Wigner, la trotterización, la disposición de los qubits y la supresión de errores, por lo que en este tutorial no hay circuitos que construir ni transpilar.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener instalados los siguientes paquetes:

• Qiskit SDK v2.3.0 o posterior (pip install qiskit)

• Qiskit Runtime v0.50.0 o posterior (pip install qiskit-ibm-runtime)

• Catálogo de Qiskit Functions v0.21.0 o posterior (pip install qiskit-ibm-catalog)

• NumPy v2.5.1 o posterior (pip install numpy)

• Matplotlib v3.10.3 o posterior (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 o posterior (pip install qctrl-visualizer)

Configuración

Primero, autentícate con tu clave de API de IBM Quantum. Luego, selecciona la Qiskit Function de la siguiente manera. (Este código supone que ya guardaste tu cuenta en tu entorno local).

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qctrlvisualizer": "qctrlvisualizer", "qiskit_ibm_catalog": "qiskit-ibm-catalog", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from qiskit_ibm_catalog import QiskitFunctionsCatalog
from qiskit_ibm_runtime import QiskitRuntimeService
import qctrlvisualizer as qv

plt.style.use(qv.get_qctrl_style())

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")

# Access Function
sim_qf = catalog.load("q-ctrl/dynamics-simulator")

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

Define el hamiltoniano

La función toma el modelo como un diccionario con tres claves. Cada una acepta un único valor float, que se aplica de manera uniforme a lo largo de la cadena, o una lista que especifica el valor sitio por sitio. Las listas tienen $L-1$ entradas para hopping_amplitude, porque una cadena abierta tiene un enlace menos que sitios, y $L$ entradas para on_site_interaction y chemical_potential.

Este tutorial usa una cadena uniforme de 31 sitios con una interacción repulsiva $U/t_h = 8$. Un potencial químico uniforme solo aporta una fase global cuando el número de partículas es fijo, por lo que este ejemplo lo establece en cero. El potencial químico importa solo cuando varía de un sitio a otro, por ejemplo, para añadir desorden o una única impureza.

site_count = 31
hopping_amplitude = 1.0
on_site_interaction = 8.0
chemical_potential = 0.0

model = {
    "hopping_amplitude": [hopping_amplitude] * (site_count - 1),
    "on_site_interaction": [on_site_interaction] * site_count,
    "chemical_potential": [chemical_potential] * site_count,
}

print(
    f"{site_count} sites, {2 * site_count} qubits, U/t_h = {on_site_interaction:.0f}"
)

31 sites, 62 qubits, U/t_h = 8

Establece el esquema de Trotter

El diccionario simulation contiene el paso de tiempo ($\Delta t$) y el número de pasos de Trotter, y el tiempo total de evolución es su producto. Pasos más grandes alcanzan tiempos posteriores con una profundidad de circuito fija, pero incurren en más error de Trotter, que crece con $U$ porque los términos del hamiltoniano conmutan menos. El estudio de hardware publicado usó $\Delta t = 0.15$ para $U/t_h = 8$ [2], y este tutorial usa el mismo valor. La profundidad del circuito crece linealmente con el número de pasos y no depende en absoluto de $L$.

time_step = 0.15
step_count = 30

simulation = {"time_step": time_step, "step_count": step_count}

print(
    f"Total evolution time: t = {time_step * step_count:.2f} (units of 1/t_h)"
)

Total evolution time: t = 4.50 (units of 1/t_h)

Prepara el estado inicial

Cada sitio contiene dos orbitales de espín, por lo que una cadena de $L$ sitios usa $2L$ qubits. El estado inicial es una cadena de ocupación de Fock con un carácter por sitio:

$$\begin{aligned} \texttt{empty (e)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (0,0), \\ \texttt{up (u)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (1,0), \\ \texttt{down (d)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (0,1), \\ \texttt{both (b)} &\longleftrightarrow \left(n_{i,\uparrow}, n_{i,\downarrow}\right) = (1,1). \end{aligned}$$

La cadena se lee empezando por el sitio menos significativo, por lo que el carácter situado más a la derecha es el sitio 0 y el índice físico del sitio del carácter $k$ contado desde la izquierda es $i = L - 1 - k$. Por lo tanto, la cadena udud significa que el sitio 0 es d, el sitio 1 es u, el sitio 2 es d y el sitio 3 es u.

A continuación, prepara un estado de Néel a medio llenar, alternando u y d, con el sitio central reemplazado por una vacante. La vacante es la perturbación cuya propagación se sigue.

vacancy_site = site_count // 2

# Build in increasing site order, then reverse: the leftmost character of the
# string is the highest site index.
site_occupations = [
    "u" if site % 2 == 0 else "d" for site in range(site_count)
]
site_occupations[vacancy_site] = "e"
initial_state = "".join(reversed(site_occupations))

print(f"Vacancy at site {vacancy_site}")
print(initial_state)

Vacancy at site 15
udududududududueudududududududu

Construye los observables

Cada observable es un diccionario que representa una suma ponderada de productos de operadores de ocupación:

$$O = \sum_{\alpha} c_{\alpha} \prod_{(\sigma, j) \in T_{\alpha}} n_{j,\sigma}.$$

Cada entrada de operators es un término $T_{\alpha}$, escrito como una lista de pares [name, site] que se multiplican entre sí, y la entrada correspondiente de coefficients es su peso. Los dos nombres son "Nup" para $n_{j,\uparrow}$ y "Ndn" para $n_{j,\downarrow}$. Un único observable de un punto $n_{3,\uparrow}$ se escribe de la siguiente manera:

{"operators": [[["Nup", 3]]], "coefficients": [1.0]}

La siguiente celda imprime el conjunto completo más pequeño de entradas que acepta la función, lo cual sirve como referencia de estructura mientras lees el código que sigue.

minimal_inputs = {
    "model": {
        "hopping_amplitude": 1.0,
        "on_site_interaction": [2.0, 2.5],
        "chemical_potential": 0.5,
    },
    "simulation": {"step_count": 5, "time_step": 0.1},
    "observables": [{"operators": [[["Nup", 0]]], "coefficients": [1.0]}],
    "initial_state": "ud",
    "shots": 4096,
}

for key, value in minimal_inputs.items():
    print(f"{key}: {value}")

model: {'hopping_amplitude': 1.0, 'on_site_interaction': [2.0, 2.5], 'chemical_potential': 0.5}
simulation: {'step_count': 5, 'time_step': 0.1}
observables: [{'operators': [[['Nup', 0]]], 'coefficients': [1.0]}]
initial_state: ud
shots: 4096

Se necesitan tres grupos de observables. La densidad de carga $n_i = n_{i\uparrow} + n_{i\downarrow}$ y el espín local $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ son sumas de dos términos sobre los mismos operadores, que difieren solo en sus coeficientes. El correlador de espín $S^z_i S^z_{i_*}$ entre el sitio $i$ y el sitio de la vacante se expande en cuatro términos de dos puntos con coeficientes $\pm 1/4$.

Los operadores de ocupación conmutan, por lo que el orden de los factores dentro de un término es irrelevante, pero el orden de la lista observables fija el orden de los valores devueltos. El código registra un índice para cada grupo mientras construye la lista. El correlador omite el sitio de la vacante, donde el producto repetiría un factor.

observables = []
charge_indices = []
spin_indices = []
correlator_indices = []

# Charge density: n_i = n_{i,up} + n_{i,down}
for site in range(site_count):
    charge_indices.append(len(observables))
    observables.append(
        {
            "operators": [[["Nup", site]], [["Ndn", site]]],
            "coefficients": [1.0, 1.0],
        }
    )

# Local spin: S^z_i = (n_{i,up} - n_{i,down}) / 2
for site in range(site_count):
    spin_indices.append(len(observables))
    observables.append(
        {
            "operators": [[["Nup", site]], [["Ndn", site]]],
            "coefficients": [0.5, -0.5],
        }
    )

# Spin correlator: S^z_i S^z_{i_*}, expanded into four two-point terms.
correlator_sites = [
    site for site in range(site_count) if site != vacancy_site
]
for site in correlator_sites:
    correlator_indices.append(len(observables))
    observables.append(
        {
            "operators": [
                [["Nup", site], ["Nup", vacancy_site]],
                [["Nup", site], ["Ndn", vacancy_site]],
                [["Ndn", site], ["Nup", vacancy_site]],
                [["Ndn", site], ["Ndn", vacancy_site]],
            ],
            "coefficients": [0.25, -0.25, -0.25, 0.25],
        }
    )

print(f"{len(observables)} observables")
print("Example correlator:", observables[correlator_indices[0]])

92 observables
Example correlator: {'operators': [[['Nup', 0], ['Nup', 15]], [['Nup', 0], ['Ndn', 15]], [['Ndn', 0], ['Nup', 15]], [['Ndn', 0], ['Ndn', 15]]], 'coefficients': [0.25, -0.25, -0.25, 0.25]}

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

Este tutorial no tiene paso de transpilación. La función realiza internamente los siguientes pasos [2], y no puedes configurarlos desde el notebook:

• Asigna los modos fermiónicos a qubits con un orden elegido para mantener bajos los pesos de Pauli en una red heavy-hex

• Usa una red SWAP fermiónica para que cada término de salto se convierta en una rotación de dos qubits

• Incrusta la cadena a lo largo de un camino que evita los qubits débiles

• Aplica desacoplamiento dinámico y twirling de Pauli durante la compilación

Tú eliges el backend y el número de shots. La función toma el nombre de un backend como cadena de texto, por lo que este tutorial selecciona el dispositivo menos ocupado lo bastante grande para la cadena y pasa su nombre.

La señal más pequeña que necesitas resolver determina el número de shots. El trazador de carga alcanza amplitudes cercanas a $0.15$, pero el trazador de espín conectado tiene un pico de alrededor de $0.04$, aproximadamente cuatro veces menor. El estudio publicado usó 20 000 shots por circuito por este motivo [2].

service = QiskitRuntimeService()
backend = service.backend("ibm_boston")
backend_name = backend.name
shot_count = 20_000

print(f"{backend_name}, {shot_count} shots per circuit")

ibm_boston, 20000 shots per circuit

Paso 3: Ejecutar con las primitivas de Qiskit

A continuación, envía la simulación. run devuelve un identificador del trabajo de inmediato en lugar de bloquear.

job = sim_qf.run(
    model=model,
    simulation=simulation,
    observables=observables,
    initial_state=initial_state,
    shots=shot_count,
    backend_name=backend_name,
)

print(job.job_id, job.status())

6261b073-be19-48a1-be60-2f369b25613c QUEUED

Usa job.status() para seguir el progreso y job.result() para bloquear hasta que el trabajo termine. El resultado es un mapeo cuya entrada observable_expectation_values se indexa primero por observable, en el orden en que se proporcionaron, y luego por tiempo. Hay step_count + 1 puntos de tiempo porque se incluye el valor en $t=0$.

result = job.result()

print("Result keys:", sorted(result.keys()))

Result keys: ['observable_expectation_values', 'observables', 'provider_job_ids', 'warnings']

Paso 4: Posprocesar y devolver el resultado en el formato clásico deseado

Reconstruir los correladores

Divide el arreglo plano de nuevo en los tres grupos usando los índices registrados anteriormente y luego forma los dos correladores trazadores. El trazador de carga es el cambio en la densidad por sitio relativo a $t=0$:

$$C^c_i(t) = \langle n_{i\uparrow}(t) + n_{i\downarrow}(t) \rangle - \langle n_{i\uparrow}(0) + n_{i\downarrow}(0) \rangle.$$

El trazador de espín es el correlador conectado entre cada sitio y el sitio de la vacante $i_*$:

$$C^s_i(t) = 4 \left( \langle S^z_i S^z_{i_*} \rangle - \langle S^z_i \rangle \langle S^z_{i_*} \rangle \right).$$

Ambos se construyen a partir de operadores en la base de Fock, por lo que provienen de las mismas cadenas de bits medidas y el segundo no cuesta ninguna ejecución adicional en el hardware.

Cuando construiste los observables, el código registró en la lista la posición de cada operador de carga (charge_indices), de espín (spin_indices) y de correlación de espín (correlator_indices). Como el simulador devuelve los valores esperados en el mismo orden, puedes usar estas posiciones para extraer los resultados correspondientes. Para las correlaciones de espín, el código también usa correlator_sites para colocar cada resultado en la posición correcta de la red y deja la vacante como NaN.

expectation_values = np.asarray(
    result["observable_expectation_values"], dtype=float
)
times = np.arange(expectation_values.shape[1]) * time_step

charge_density = expectation_values[charge_indices]
local_spin = expectation_values[spin_indices]
spin_correlator = np.full_like(charge_density, np.nan)
spin_correlator[correlator_sites] = expectation_values[correlator_indices]

# Charge tracer: change in per-site density relative to t = 0.
charge_tracer = charge_density - charge_density[:, [0]]

# Spin tracer: connected correlator with the vacancy site.
spin_tracer = 4.0 * (spin_correlator - local_spin * local_spin[vacancy_site])
spin_tracer[vacancy_site] = np.nan

print(
    f"{expectation_values.shape[0]} observables, "
    f"{expectation_values.shape[1]} time points up to t = {times[-1]:.2f}"
)
print(
    f"Peak amplitudes: charge {np.nanmax(np.abs(charge_tracer)):.3f}, "
    f"spin {np.nanmax(np.abs(spin_tracer)):.3f}"
)

92 observables, 31 time points up to t = 4.50
Peak amplitudes: charge 0.978, spin 0.168

Visualizar los frentes de onda

Si representas cada correlador como un mapa de calor del espacio y el tiempo, el cono de luz se ve directamente. El tiempo fluye hacia abajo y el eje horizontal es la posición a lo largo de la cadena.

def plot_tracer(axis, tracer, title):
    """Draw one tracer correlator as a space-time heatmap, time flowing down."""
    # Scale from the off-centre sites. The vacancy site carries an order-one
    # change in charge density and is left to clip rather than set the scale.
    off_center = np.ones(site_count, dtype=bool)
    off_center[vacancy_site] = False
    limit = np.nanpercentile(np.abs(tracer[off_center]), 98)

    colormap = qv.QCTRL_DIVERGENT_COLORMAP.reversed().copy()
    colormap.set_bad("0.82")  # the spin tracer has no value at the vacancy

    image = axis.imshow(
        tracer.T,
        aspect="auto",
        origin="upper",
        cmap=colormap,
        vmin=-limit,
        vmax=limit,
        extent=(-0.5, site_count - 0.5, times[-1], times[0]),
        interpolation="nearest",
    )

    axis.annotate(
        "vacancy",
        xy=(vacancy_site, times[-1]),
        xytext=(vacancy_site, times[-1] + 0.55),
        ha="center",
        va="top",
        fontsize=8,
        annotation_clip=False,
        arrowprops=dict(arrowstyle="-", lw=0.8, color="0.4"),
    )

    axis.set_title(title)
    axis.set_xlabel("Site")
    return image

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)

charge_image = plot_tracer(
    axes[0], charge_tracer, r"Charge tracer $C^{c}_{i}(t)$"
)
spin_image = plot_tracer(axes[1], spin_tracer, r"Spin tracer $C^{s}_{i}(t)$")

axes[0].set_ylabel(r"Time $t$  $(1/t_h)$")
fig.colorbar(charge_image, ax=axes[0], fraction=0.046)
fig.colorbar(spin_image, ax=axes[1], fraction=0.046)
fig.suptitle(f"$L$ = {site_count},  $U/t_h$ = {on_site_interaction:.0f}")
fig.tight_layout()
plt.show()

Extraer las velocidades de los frentes de onda

Convierte los mapas de calor en números siguiendo el protocolo de detección del estudio publicado [2]. En cada instante de tiempo, simetriza el perfil espacial respecto a la vacante, lo que promedia la parte antisimétrica del ruido del hardware, normaliza por la mayor amplitud fuera del centro y registra el sitio más externo por encima de un umbral. Luego las detecciones se dividen en tramos de pasos de tiempo consecutivos y solo se ajusta el tramo más largo, lo que mantiene fuera del resultado el transitorio de tiempos tempranos.

def front_distance(profile, center, threshold):
    """Distance from `center` to the outermost site above a relative threshold."""
    values = np.nan_to_num(np.asarray(profile, dtype=float))
    length = values.size

    symmetrized = values.copy()
    for site in range(length):
        mirror = 2 * center - site
        if 0 <= mirror < length:
            symmetrized[site] = 0.5 * (values[site] + values[mirror])

    amplitude = np.abs(symmetrized)
    off_center = np.ones(length, dtype=bool)
    off_center[center] = False

    peak = amplitude[off_center].max()
    if peak <= 0.0:
        return np.nan
    above = np.flatnonzero((amplitude / peak > threshold) & off_center)
    return float(np.abs(above - center).max()) if above.size else np.nan

def longest_run(mask):
    """Indices of the longest run of consecutive True values."""
    best, current = [], []
    for index, keep in enumerate(mask):
        current = current + [index] if keep else []
        if len(current) > len(best):
            best = current
    return np.array(best, dtype=int)

def front_velocity(tracer, center, threshold=0.3, quantile=0.75, min_block=5):
    """Fit front position against time, using only the longest contiguous run."""
    distances = np.array(
        [
            front_distance(tracer[:, k], center, threshold)
            for k in range(tracer.shape[1])
        ]
    )
    edge = max(center, tracer.shape[0] - 1 - center)
    usable = (
        np.isfinite(distances) & (distances >= 2) & (distances <= edge - 2)
    )

    block = longest_run(usable)
    if block.size < min_block:
        return np.nan, distances, np.array([], dtype=int), np.nan

    candidates = [
        stats.theilslopes(distances[block[:k]], times[block[:k]]).slope
        for k in range(min_block, block.size + 1)
    ]
    cutoff = np.quantile(candidates, quantile)
    keep = max(i for i, value in enumerate(candidates) if value >= cutoff)

    fit_block = block[: keep + min_block]
    slope, intercept = np.polyfit(times[fit_block], distances[fit_block], 1)
    return float(slope), distances, fit_block, float(intercept)

La velocidad ajustada depende del umbral de detección y del cuantil de recorte, y esa sensibilidad domina sobre el ruido de disparo. Barre ambos e informa la mediana con la envolvente completa, como hace el análisis publicado.

thresholds = [0.15, 0.225, 0.30, 0.375, 0.45]
quantiles = [0.60, 0.675, 0.75, 0.825, 0.90]

def velocity_envelope(tracer, center):
    """Median, minimum, and maximum velocity over the parameter scan."""
    values = [
        front_velocity(tracer, center, threshold=p, quantile=q)[0]
        for p in thresholds
        for q in quantiles
    ]
    values = np.array([v for v in values if np.isfinite(v)])
    if values.size == 0:
        return np.nan, np.nan, np.nan
    return float(np.median(values)), float(values.min()), float(values.max())

charge_velocity, charge_low, charge_high = velocity_envelope(
    charge_tracer, vacancy_site
)
spin_velocity, spin_low, spin_high = velocity_envelope(
    spin_tracer, vacancy_site
)

_, charge_front, charge_block, charge_intercept = front_velocity(
    charge_tracer, vacancy_site
)
_, spin_front, spin_block, spin_intercept = front_velocity(
    spin_tracer, vacancy_site
)

print(f"v_c = {charge_velocity:.2f}  [{charge_low:.2f}, {charge_high:.2f}]")
print(f"v_s = {spin_velocity:.2f}  [{spin_low:.2f}, {spin_high:.2f}]")
print(f"v_s / v_c = {spin_velocity / charge_velocity:.2f}")

v_c = 1.93  [1.81, 2.18]
v_s = 1.61  [1.23, 2.12]
v_s / v_c = 0.83

PUBLISHED = {"charge": 1.78, "spin": 1.46}

edge = max(vacancy_site, site_count - 1 - vacancy_site)
fig, axis = plt.subplots(figsize=(7, 4.6))

axis.axhspan(0, 2, color="0.92", zorder=0)
axis.axhspan(edge - 2, edge + 2, color="0.92", zorder=0)

series = [
    (
        "charge",
        "$v_c$",
        charge_front,
        charge_block,
        charge_velocity,
        charge_intercept,
        charge_low,
        charge_high,
        "C0",
        "o",
    ),
    (
        "spin",
        "$v_s$",
        spin_front,
        spin_block,
        spin_velocity,
        spin_intercept,
        spin_low,
        spin_high,
        "C1",
        "s",
    ),
]

for name, sym, front, block, v, b, lo, hi, color, marker in series:
    excluded = np.setdiff1d(np.arange(front.size), block)
    axis.plot(
        times[excluded],
        front[excluded],
        marker,
        color=color,
        markersize=5,
        markerfacecolor="none",
        alpha=0.45,
        zorder=2,
    )
    axis.plot(
        times[block],
        front[block],
        marker,
        color=color,
        markersize=5,
        label=f"{name}, {sym} = {v:.2f} [{lo:.2f}, {hi:.2f}]",
        zorder=3,
    )

    span = times[block]
    anchor = v * span[0] + b
    axis.plot(
        span,
        v * (span - span[0]) + anchor,
        "-",
        color=color,
        lw=1.6,
        zorder=4,
    )
    axis.fill_between(
        span,
        lo * (span - span[0]) + anchor,
        hi * (span - span[0]) + anchor,
        color=color,
        alpha=0.15,
        zorder=1,
    )
    axis.plot(
        span,
        PUBLISHED[name] * (span - span[0]) + anchor,
        "--",
        color=color,
        lw=1.2,
        alpha=0.8,
        zorder=4,
    )

axis.plot([], [], "--", color="0.4", lw=1.2, label="published [2]")
axis.plot(
    [],
    [],
    "o",
    color="0.5",
    markerfacecolor="none",
    label="excluded from fit",
)

axis.set_xlabel(r"Time $t$  $(1/t_h)$")
axis.set_ylabel("Front distance from vacancy (sites)")
axis.set_xlim(0, times[-1])
axis.set_ylim(
    0, min(edge + 1, np.nanmax(np.r_[charge_front, spin_front]) + 3)
)
axis.legend(loc="upper left", frameon=False, fontsize=9)
axis.set_title(
    f"L = {site_count}, $U/t_h$ = {on_site_interaction:.0f}, front propagation",
    fontsize=11,
)

fig.tight_layout()
plt.show()

Dibujar las detecciones sobre los mapas de calor muestra si el algoritmo siguió el borde delantero. Los marcadores rellenos son las detecciones que entraron en el ajuste. Los marcadores abiertos fueron descartados por las reglas de enmascaramiento, ya sea porque están demasiado cerca de la vacante, demasiado cerca del extremo de la cadena, o porque son lo bastante tardíos como para que el frente se haya ensanchado y su posición ya no sea nítida. La línea discontinua es el ajuste sobre los tiempos que usó, y su pendiente es la velocidad indicada en el título de cada panel.

La columna de la vacante está marcada debajo del eje. La carga se satura allí porque su valor es del orden de uno, y el espín no tiene valor allí, por lo que se dibuja en gris.

El cono de carga se abre más que el cono de espín. Ambas perturbaciones dejan el mismo defecto en el mismo momento y viajan hacia afuera a velocidades diferentes.

def plot_tracer_front(axis, tracer, front, block, velocity, intercept, title):
    """Overlay the detected wavefront and the fitted cone on a tracer heatmap."""
    image = plot_tracer(axis, tracer, title)

    excluded = np.setdiff1d(np.arange(front.size), block)
    for index, filled in ((block, True), (excluded, False)):
        good = index[np.isfinite(front[index])]
        if good.size == 0:
            continue
        for sign in (-1, 1):
            axis.plot(
                vacancy_site + sign * front[good],
                times[good],
                "o",
                markersize=4,
                color="k",
                markerfacecolor="k" if filled else "none",
                markeredgewidth=0.9,
                alpha=1.0 if filled else 0.4,
            )

    span = times[block]
    for sign in (-1, 1):
        axis.plot(
            vacancy_site + sign * (velocity * span + intercept),
            span,
            "--",
            color="k",
            lw=1.5,
        )

    axis.set_title(f"{title}    $v$ = {velocity:.2f}")
    return image

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharey=True)

charge_image = plot_tracer_front(
    axes[0],
    charge_tracer,
    charge_front,
    charge_block,
    charge_velocity,
    charge_intercept,
    r"Charge tracer $C^{c}_{i}(t)$",
)
spin_image = plot_tracer_front(
    axes[1],
    spin_tracer,
    spin_front,
    spin_block,
    spin_velocity,
    spin_intercept,
    r"Spin tracer $C^{s}_{i}(t)$",
)

axes[0].set_ylabel(r"Time $t$  $(1/t_h)$")
fig.colorbar(charge_image, ax=axes[0], fraction=0.046)
fig.colorbar(spin_image, ax=axes[1], fraction=0.046)
fig.suptitle(f"$L$ = {site_count},  $U/t_h$ = {on_site_interaction:.0f}")
fig.tight_layout()
plt.show()

El frente de carga adelanta al frente de espín. Un hueco en un fondo de Mott unidimensional no viaja como un electrón. En cambio, se fracciona en un holón que lleva la carga y un espinón que lleva la perturbación de espín, y con $U/t_h = 8$ el espinón es el más lento de los dos.

El orden depende de la interacción. Con $U = 0$, el modelo son fermiones libres con $v_c = 2$ y $v_s = 4$, por lo que el espín es más rápido. La repulsión suprime el transporte de espín dejando casi intacto el transporte de carga, y para $U/t_h = 8$ los valores publicados en hardware son $v_c = 1.78$ y $v_s = 1.46$ [2]. Repetir la ejecución con varios valores de $U$ y graficar $v_s / v_c$ permite mapear el cruce.

Ten en cuenta estas salvedades sobre los números. La detección basada en umbral tiene un sesgo ligeramente a la baja, por lo que es de esperar una velocidad de carga unos pocos puntos porcentuales por debajo del valor publicado. La envolvente de espín es ancha porque el correlador conectado está más cerca del piso de ruido, y por eso se cita la mediana del barrido en lugar de un único ajuste.

Referencias

[1] Essler, F. H. L., Frahm, H., Gohmann, F., Klumper, A., & Korepin, V. E. (2005). The One-Dimensional Hubbard Model. Cambridge University Press.

[2] Hartnett, G. S., Najafi, K. S., Khindanov, A., Liao, H., Schutzman, M., Hush, M. R., Biercuk, M. J., & Baum, Y. (2026). Fast, accurate, high-resolution simulation of large-scale Fermi-Hubbard models on a digital quantum processor. arXiv:2605.04025.