Dinamica di spin e di carica nel modello di Fermi-Hubbard

Stima di utilizzo: 5 minuti su un processore Heron r3. (NOTA: si tratta solo di una stima. Il tuo tempo di esecuzione potrebbe variare.)

Contesto

Il modello di Fermi-Hubbard descrive elettroni che saltano su un reticolo con una penalità energetica per la doppia occupazione. Su una catena unidimensionale di $L$ siti con condizioni al contorno aperte, l'Hamiltoniana è definita come segue:

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

dove $c^\dagger_{i\sigma}$ crea un fermione di spin $\sigma$ sul sito $i$ e $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ li conta. I tre termini sono l'ampiezza di hopping $t_{h,i}$ sul legame tra i siti $i$ e $i+1$, l'interazione on-site $U_i$ e il potenziale chimico $\mu_i$. Il tutorial misura l'energia in unità di $t_h$ e il tempo in unità di $1/t_h$.

La dinamica in tempo reale di questo modello è difficile da calcolare in modo classico. La diagonalizzazione esatta si ferma intorno a $L \lesssim 14$ siti e, dopo un quench, l'entropia di entanglement di un solutore a rete di tensori cresce linearmente nel tempo, quindi la dimensione di legame necessaria per un'accuratezza fissa cresce in modo esponenziale.

In una dimensione, le eccitazioni elementari sono gli spinoni, che trasportano spin e non carica, e gli olonni, che trasportano carica e non spin. I due si propagano a velocità diverse, quindi una perturbazione localizzata si divide in un fronte veloce e uno lento. Questo effetto è chiamato separazione spin-carica, e questo tutorial lo misura [1].

Questo tutorial usa la Qiskit Function Dynamics Simulator, che prende l'Hamiltoniana, uno stato iniziale, uno schema di Trotter e un elenco di osservabili, e restituisce valori di aspettazione. Gestisce internamente la mappatura di Jordan-Wigner, la trotterizzazione, il layout dei qubit e la soppressione degli errori, quindi in questo tutorial non ci sono circuiti da costruire o traspilare.

Requisiti

Prima di iniziare questo tutorial, assicurati di avere installato i seguenti pacchetti:

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

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

• Qiskit Functions Catalog v0.21.0 o successivo (pip install qiskit-ibm-catalog)

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

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

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

Configurazione

Per prima cosa, autenticati usando la tua chiave API di IBM Quantum. Poi seleziona la Qiskit Function come segue. (Questo codice presuppone che tu abbia già salvato il tuo account nel tuo ambiente locale.)

# 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")

Passo 1: Mappare gli input classici su un problema quantistico

Definisci l'Hamiltoniana

La funzione accetta il modello come dizionario con tre chiavi. Ognuna accetta un singolo float, che si applica uniformemente lungo la catena, oppure una lista che specifica il valore sito per sito. Le liste hanno $L-1$ elementi per hopping_amplitude, perché una catena aperta ha un legame in meno rispetto ai siti, e $L$ elementi per on_site_interaction e chemical_potential.

Questo tutorial usa una catena uniforme di 31 siti con un'interazione repulsiva $U/t_h = 8$. Un potenziale chimico uniforme contribuisce solo con una fase globale quando il numero di particelle è fissato, quindi questo esempio lo pone a zero. Il potenziale chimico conta solo quando varia da sito a sito, ad esempio per aggiungere disordine o una singola impurità.

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

Imposta lo schema di Trotter

Il dizionario simulation contiene il passo temporale ($\Delta t$) e il numero di passi di Trotter, e il tempo di evoluzione totale è il loro prodotto. Passi più grandi raggiungono tempi più avanzati a profondità del circuito fissa ma comportano un errore di Trotter maggiore, che cresce con $U$ perché i termini dell'Hamiltoniana commutano meno. Lo studio hardware pubblicato ha usato $\Delta t = 0.15$ per $U/t_h = 8$ [2], e questo tutorial usa lo stesso valore. La profondità del circuito cresce linearmente con il numero di passi e non dipende affatto da $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 lo stato iniziale

Ogni sito ospita due spin-orbitali, quindi una catena di $L$ siti usa $2L$ qubit. Lo stato iniziale è una stringa di occupazione di Fock con un carattere per sito:

$$\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 stringa viene letta partendo dal sito meno significativo, quindi il carattere più a destra è il sito 0 e l'indice fisico del sito del carattere $k$ contato da sinistra è $i = L - 1 - k$. La stringa udud significa quindi che il sito 0 è d, il sito 1 è u, il sito 2 è d e il sito 3 è u.

Poi, prepara uno stato di Néel a riempimento metà, alternando u e d, con il sito centrale sostituito da una vacanza. La vacanza è la perturbazione di cui si traccia la diffusione.

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

Costruisci le osservabili

Ogni osservabile è un dizionario che rappresenta una somma pesata di prodotti di operatori di occupazione:

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

Ogni elemento di operators è un termine $T_{\alpha}$, scritto come lista di coppie [name, site] che si moltiplicano tra loro, e l'elemento corrispondente di coefficients è il suo peso. I due nomi sono "Nup" per $n_{j,\uparrow}$ e "Ndn" per $n_{j,\downarrow}$. Una singola osservabile a un punto $n_{3,\uparrow}$ si scrive come segue:

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

La cella successiva stampa il più piccolo insieme completo di input che la funzione accetta, utile come riferimento della struttura mentre leggi il codice che segue.

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

Servono tre gruppi di osservabili. La densità di carica $n_i = n_{i\uparrow} + n_{i\downarrow}$ e lo spin locale $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ sono entrambe somme di due termini sugli stessi operatori, che differiscono solo per i coefficienti. Il correlatore di spin $S^z_i S^z_{i_*}$ tra il sito $i$ e il sito della vacanza si espande in quattro termini a due punti con coefficienti $\pm 1/4$.

Gli operatori di occupazione commutano, quindi l'ordine dei fattori all'interno di un termine è irrilevante, ma l'ordine della lista observables fissa l'ordine dei valori restituiti. Il codice registra un indice per ogni gruppo mentre costruisce la lista. Il correlatore salta il sito della vacanza, dove il prodotto ripeterebbe un fattore.

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

Passo 2: Ottimizzare il problema per l'esecuzione su hardware quantistico

Questo tutorial non ha un passo di transpilazione. La funzione esegue internamente i seguenti passi [2], e non puoi configurarli dal notebook:

• Mappa i modi fermionici sui qubit secondo un ordinamento scelto per mantenere bassi i pesi di Pauli su un reticolo heavy-hex

• Usa una rete SWAP fermionica in modo che ogni termine di hopping diventi una rotazione a due qubit

• Incorpora la catena lungo un percorso che evita i qubit deboli

• Applica dynamical decoupling e Pauli twirling durante la compilazione

Scegli tu il backend e il numero di shot. La funzione accetta un nome di backend come stringa, quindi questo tutorial seleziona il dispositivo meno occupato abbastanza grande per la catena e ne passa il nome.

Il segnale più piccolo che devi risolvere determina il numero di shot. Il tracciante di carica raggiunge ampiezze vicine a $0.15$, ma il tracciante di spin connesso ha un picco intorno a $0.04$, circa quattro volte più piccolo. Per questo motivo lo studio pubblicato ha usato 20.000 shot per circuito [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

Passo 3: Eseguire usando le primitive di Qiskit

Ora invia la simulazione. run restituisce subito un handle del job, senza bloccarsi.

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() per seguire l'avanzamento e job.result() per bloccarti finché il job non termina. Il risultato è una mappatura la cui voce observable_expectation_values è indicizzata prima per osservabile, nell'ordine in cui sono state fornite, e poi per tempo. Ci sono step_count + 1 punti temporali perché è incluso il valore a $t=0$.

result = job.result()

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

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

Passo 4: Post-elaborare e restituire il risultato nel formato classico desiderato

Ricostruisci i correlatori

Dividi l'array piatto nei tre gruppi usando gli indici registrati in precedenza, poi forma i due correlatori traccianti. Il tracciante di carica è la variazione della densità per sito rispetto 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.$$

Il tracciante di spin è il correlatore connesso tra ogni sito e il sito della vacanza $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).$$

Entrambi sono costruiti da operatori nella base di Fock, quindi provengono dalle stesse bitstring misurate e il secondo non richiede alcuna esecuzione hardware aggiuntiva.

Quando hai costruito le osservabili, il codice ha registrato nella lista la posizione di ogni operatore di carica (charge_indices), di spin (spin_indices) e di correlazione di spin (correlator_indices). Poiché il simulatore restituisce i valori di aspettazione nello stesso ordine, puoi usare queste posizioni per estrarre i risultati corrispondenti. Per le correlazioni di spin, il codice usa inoltre correlator_sites per collocare ogni risultato nella posizione corretta del reticolo, e lascia la vacanza come 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

Visualizza i fronti d'onda

Se rappresenti ogni correlatore come una mappa di calore di spazio e tempo, il cono di luce è direttamente visibile. Il tempo scorre verso il basso e l'asse orizzontale è la posizione lungo la catena.

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

Estrai le velocità dei fronti d'onda

Trasforma le heatmap in numeri seguendo il protocollo di rilevamento dello studio pubblicato [2]. Ad ogni intervallo di tempo, simmetrizza il profilo spaziale rispetto alla vacanza, il che elimina per media la parte antisimmetrica del rumore dell'hardware, normalizza per la maggiore ampiezza fuori dal centro e registra il sito più esterno al di sopra di una soglia. I rilevamenti vengono poi suddivisi in sequenze di passi temporali consecutivi e si esegue il fit solo sulla sequenza più lunga, il che esclude dal risultato il transiente dei tempi iniziali.

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 velocità ottenuta dal fit dipende dalla soglia di rilevamento e dal quantile di troncamento, e questa sensibilità domina sul rumore di shot. Esegui la scansione di entrambi e riporta la mediana con l'intero inviluppo, come fa l'analisi pubblicata.

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

Ridisegnare i rilevamenti sulle heatmap mostra se l'algoritmo ha seguito il fronte di propagazione. I marcatori pieni sono i rilevamenti entrati nel fit. I marcatori vuoti sono stati scartati dalle regole di mascheramento, perché si trovano troppo vicino alla vacanza, troppo vicino all'estremità della catena, oppure erano abbastanza tardivi da far sì che il fronte si sia allargato e la sua posizione non sia più netta. La linea tratteggiata è il fit sui tempi utilizzati e la sua pendenza è la velocità riportata nel titolo di ogni pannello.

La colonna della vacanza è indicata sotto l'asse. La carica viene troncata lì perché il suo valore è dell'ordine di uno, mentre lo spin non ha alcun valore in quel punto, quindi è disegnato in grigio.

Il cono della carica si apre più ampiamente di quello dello spin. Entrambe le perturbazioni partono dallo stesso difetto nello stesso istante e si propagano verso l'esterno a velocità diverse.

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

Il fronte della carica supera quello dello spin. Una lacuna in uno sfondo di Mott unidimensionale non si propaga come un elettrone. Si frazionalizza invece in un olone che porta la carica e uno spinone che porta la perturbazione dello spin, e a $U/t_h = 8$ lo spinone è il più lento dei due.

L'ordine dipende dall'interazione. A $U = 0$, il modello è di fermioni liberi con $v_c = 2$ e $v_s = 4$, quindi lo spin è più veloce. La repulsione sopprime il trasporto dello spin lasciando quasi intatto il trasporto della carica, e a $U/t_h = 8$ i valori hardware pubblicati sono $v_c = 1.78$ e $v_s = 1.46$ [2]. Ripetere l'esecuzione con diversi valori di $U$ e tracciare $v_s / v_c$ mostra il crossover.

Tieni presenti queste avvertenze sui numeri. Il rilevamento basato su soglia è leggermente distorto verso il basso, quindi è atteso un valore della velocità della carica inferiore di qualche punto percentuale a quello pubblicato. L'inviluppo dello spin è ampio perché il correlatore connesso è più vicino al livello di rumore di fondo, ed è per questo che si riporta la mediana della scansione anziché un singolo fit.

Riferimenti

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