Dinâmica de spin e de carga no modelo de Fermi-Hubbard

Estimativa de uso: 5 minutos em um processador Heron r3. (OBSERVAÇÃO: Esta é apenas uma estimativa. Seu tempo de execução pode variar.)

Contexto

O modelo de Fermi-Hubbard descreve elétrons saltando em uma rede, com uma penalidade de energia para dupla ocupação. Em uma cadeia unidimensional de $L$ sítios com fronteiras abertas, o Hamiltoniano é definido da seguinte forma:

$$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},$$

onde $c^\dagger_{i\sigma}$ cria um férmion de spin $\sigma$ no sítio $i$ e $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ os conta. Os três termos são a amplitude de salto $t_{h,i}$ na ligação entre os sítios $i$ e $i+1$, a interação local $U_i$ e o potencial químico $\mu_i$. O tutorial mede a energia em unidades de $t_h$ e o tempo em unidades de $1/t_h$.

A dinâmica em tempo real deste modelo é difícil de calcular classicamente. A diagonalização exata para por volta de $L \lesssim 14$ sítios e, depois de um quench, a entropia de emaranhamento de um solucionador de rede tensorial cresce linearmente no tempo, de modo que a dimensão de ligação necessária para uma acurácia fixa cresce exponencialmente.

Em uma dimensão, as excitações elementares são os spinons, que carregam spin e nenhuma carga, e os holons, que carregam carga e nenhum spin. Os dois se propagam a velocidades diferentes, então uma perturbação localizada se divide em uma frente rápida e uma lenta. Esse efeito é chamado de separação spin-carga, e este tutorial o mede [1].

Este tutorial usa a Qiskit Function Dynamics Simulator, que recebe o Hamiltoniano, um estado inicial, um cronograma de Trotter e uma lista de observáveis, e retorna valores esperados. Ela cuida internamente do mapeamento de Jordan-Wigner, da trotterização, do layout de qubits e da supressão de erros, então não há circuitos para construir ou transpilar neste tutorial.

Requisitos

Antes de começar este tutorial, certifique-se de que você tem os seguintes pacotes instalados:

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

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

• Qiskit Functions Catalog v0.21.0 ou posterior (pip install qiskit-ibm-catalog)

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

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

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

Configuração

Primeiro, autentique-se usando sua chave de API do IBM Quantum. Depois, selecione a Qiskit Function da seguinte forma. (Este código pressupõe que você já salvou sua conta no seu ambiente 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")

Etapa 1: Mapear entradas clássicas para um problema quântico

Definir o Hamiltoniano

A função recebe o modelo como um dicionário com três chaves. Cada uma aceita um único float, que se aplica uniformemente ao longo da cadeia, ou uma lista que especifica o valor sítio a sítio. As listas têm $L-1$ entradas para hopping_amplitude, porque uma cadeia aberta tem uma ligação a menos que sítios, e $L$ entradas para on_site_interaction e chemical_potential.

Este tutorial usa uma cadeia uniforme de 31 sítios com uma interação repulsiva $U/t_h = 8$. Um potencial químico uniforme contribui apenas com uma fase global quando o número de partículas é fixo, então este exemplo o define como zero. O potencial químico importa apenas quando varia de sítio para sítio, por exemplo, para adicionar desordem ou uma ú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

Definir o cronograma de Trotter

O dicionário simulation contém o passo de tempo ($\Delta t$) e o número de passos de Trotter, e o tempo total de evolução é o produto deles. Passos maiores alcançam tempos posteriores com profundidade de circuito fixa, mas incorrem em mais erro de Trotter, que cresce com $U$ porque os termos do Hamiltoniano comutam menos. O estudo de hardware publicado usou $\Delta t = 0.15$ para $U/t_h = 8$ [2], e este tutorial usa o mesmo valor. A profundidade do circuito cresce linearmente com o número de passos e não depende de $L$ de forma alguma.

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)

Preparar o estado inicial

Cada sítio contém dois orbitais de spin, então uma cadeia de $L$ sítios usa $2L$ qubits. O estado inicial é uma string de ocupação de Fock com um caractere por sítio:

$$\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}$$

A string é lida começando pelo sítio menos significativo, então o caractere mais à direita é o sítio 0 e o índice físico do sítio do caractere $k$ contado a partir da esquerda é $i = L - 1 - k$. A string udud significa, portanto, que o sítio 0 é d, o sítio 1 é u, o sítio 2 é d e o sítio 3 é u.

Em seguida, prepare um estado de Néel de preenchimento médio, alternando u e d, com o sítio central substituído por uma vacância. A vacância é a perturbação cuja propagação é acompanhada.

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

Construir os observáveis

Cada observável é um dicionário que representa uma soma ponderada de produtos de operadores de ocupação:

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

Cada entrada de operators é um termo $T_{\alpha}$, escrito como uma lista de pares [name, site] que se multiplicam, e a entrada correspondente de coefficients é o seu peso. Os dois nomes são "Nup" para $n_{j,\uparrow}$ e "Ndn" para $n_{j,\downarrow}$. Um único observável de um ponto $n_{3,\uparrow}$ é escrito da seguinte forma:

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

A próxima célula imprime o menor conjunto completo de entradas que a função aceita, o que é útil como referência de formato ao ler o código a seguir.

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

São necessários três grupos de observáveis. A densidade de carga $n_i = n_{i\uparrow} + n_{i\downarrow}$ e o spin local $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ são ambos somas de dois termos sobre os mesmos operadores, diferindo apenas nos coeficientes. O correlador de spin $S^z_i S^z_{i_*}$ entre o sítio $i$ e o sítio da vacância se expande em quatro termos de dois pontos com coeficientes $\pm 1/4$.

Os operadores de ocupação comutam, então a ordem dos fatores dentro de um termo é irrelevante, mas a ordem da lista observables fixa a ordem dos valores retornados. O código registra um índice para cada grupo enquanto constrói a lista. O correlador pula o sítio da vacância, onde o produto repetiria um fator.

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]}

Etapa 2: Otimizar o problema para execução em hardware quântico

Este tutorial não tem etapa de transpilação. A função realiza internamente as seguintes etapas [2], e você não pode configurá-las a partir do notebook:

• Mapeia os modos fermiônicos para qubits sob uma ordenação escolhida para manter baixos os pesos de Pauli em uma rede heavy-hex

• Usa uma rede de SWAP fermiônica para que cada termo de salto se torne uma rotação de dois qubits

• Incorpora a cadeia ao longo de um caminho que evita qubits fracos

• Aplica desacoplamento dinâmico e Pauli twirling durante a compilação

Você escolhe o backend e o número de shots. A função recebe o nome de um backend como string, então este tutorial seleciona o dispositivo menos ocupado grande o bastante para a cadeia e passa o nome dele.

O menor sinal que você precisa resolver define o número de shots. O traçador de carga atinge amplitudes próximas de $0.15$, mas o traçador de spin conectado tem pico em torno de $0.04$, cerca de quatro vezes menor. O estudo publicado usou 20.000 shots por circuito por esse 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

Etapa 3: Executar usando primitivas do Qiskit

Em seguida, envie a simulação. run retorna um identificador do job imediatamente, em vez 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

Use job.status() para acompanhar o progresso e job.result() para bloquear até o job terminar. O resultado é um mapeamento cuja entrada observable_expectation_values é indexada primeiro por observável, na ordem em que foram fornecidos, e depois por tempo. Há step_count + 1 pontos de tempo porque o valor em $t=0$ está incluído.

result = job.result()

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

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

Etapa 4: Pós-processar e retornar o resultado no formato clássico desejado

Reconstruir os correlatores

Divida o array plano de volta nos três grupos usando os índices registrados anteriormente e, em seguida, forme os dois correlatores traçadores. O traçador de carga é a variação da densidade por sítio em relação 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.$$

O traçador de spin é o correlator conectado entre cada sítio e o sítio da vacância $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 são construídos a partir de operadores na base de Fock, então vêm dos mesmos bitstrings medidos e o segundo não custa nenhuma execução extra de hardware.

Ao construir os observáveis, o código registrou a posição de cada operador de carga (charge_indices), de spin (spin_indices) e de correlação de spin (correlator_indices) na lista. Como o simulador retorna os valores esperados na mesma ordem, você pode usar essas posições para extrair os resultados correspondentes. Para as correlações de spin, o código também usa correlator_sites para colocar cada resultado na posição correta da rede e deixa a vacância 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 as frentes de onda

Se você plotar cada correlator como um mapa de calor do espaço e do tempo, o cone de luz fica diretamente visível. O tempo flui para baixo, e o eixo horizontal é a posição ao longo da cadeia.

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()

Extrair as velocidades das frentes de onda

Transforme os mapas de calor em números seguindo o protocolo de detecção do estudo publicado [2]. Em cada fatia de tempo, simetrize o perfil espacial em torno da vacância, o que elimina por média a parte antissimétrica do ruído do hardware, normalize pela maior amplitude fora do centro e registre o sítio mais externo acima de um limiar. As detecções são então divididas em trechos de passos de tempo consecutivos, e apenas o trecho mais longo é ajustado, o que mantém o transiente de tempos iniciais fora do resultado.

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)

A velocidade ajustada depende do limiar de detecção e do quantil de corte, e essa sensibilidade domina o ruído de shots. Varra ambos e informe a mediana com o envelope completo, como faz a análise publicada.

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()

Desenhar as detecções de volta sobre os mapas de calor mostra se o algoritmo acompanhou a borda dianteira. Os marcadores preenchidos são as detecções que entraram no ajuste. Os marcadores vazios foram cortados pelas regras de máscara, seja porque estão muito perto da vacância, muito perto do fim da cadeia, ou eram tardios o bastante para que a frente já tenha se alargado e sua posição não seja mais nítida. A linha tracejada é o ajuste sobre os tempos que ele usou, e sua inclinação é a velocidade no título de cada painel.

A coluna da vacância está marcada abaixo do eixo. A carga é cortada ali porque seu valor é da ordem de um, e o spin não tem valor ali, por isso é desenhado em cinza.

O cone de carga se abre mais do que o cone de spin. Ambas as perturbações partem do mesmo defeito no mesmo instante e se propagam para fora 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()

A frente de carga ultrapassa a frente de spin. Um buraco em um fundo de Mott unidimensional não viaja como um elétron. Em vez disso, ele se fracionaliza em um holon que carrega a carga e um spinon que carrega a perturbação de spin, e em $U/t_h = 8$ o spinon é o mais lento dos dois.

A ordem depende da interação. Em $U = 0$, o modelo é de férmions livres com $v_c = 2$ e $v_s = 4$, então o spin é mais rápido. A repulsão suprime o transporte de spin enquanto deixa o transporte de carga quase intocado, e em $U/t_h = 8$ os valores publicados de hardware são $v_c = 1.78$ e $v_s = 1.46$ [2]. Executar novamente com vários valores de $U$ e plotar $v_s / v_c$ mapeia a transição.

Observe estas ressalvas sobre os números. A detecção baseada em limiar tem um viés ligeiramente para baixo, então é esperada uma velocidade de carga alguns por cento abaixo do valor publicado. O envelope de spin é largo porque o correlator conectado fica mais perto do piso de ruído, e por isso se cita a mediana da varredura em vez de um único ajuste.

Referências

[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.