Dynamique de spin et de charge dans le modèle de Fermi-Hubbard

Estimation d'utilisation : 5 minutes sur un processeur Heron r3. (REMARQUE : il s'agit uniquement d'une estimation. Ton temps d'exécution peut varier.)

Contexte

Le modèle de Fermi-Hubbard décrit des électrons qui sautent sur un réseau avec une pénalité d'énergie en cas de double occupation. Sur une chaîne unidimensionnelle de $L$ sites avec des conditions aux limites ouvertes, le hamiltonien est défini comme suit :

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

où $c^\dagger_{i\sigma}$ crée un fermion de spin $\sigma$ sur le site $i$ et $n_{i\sigma} = c^\dagger_{i\sigma} c_{i\sigma}$ les compte. Les trois termes sont l'amplitude de saut $t_{h,i}$ sur la liaison entre les sites $i$ et $i+1$, l'interaction sur site $U_i$ et le potentiel chimique $\mu_i$. Le tutoriel mesure l'énergie en unités de $t_h$ et le temps en unités de $1/t_h$.

La dynamique en temps réel de ce modèle est difficile à calculer de manière classique. La diagonalisation exacte s'arrête autour de $L \lesssim 14$ sites, et après un quench, l'entropie d'intrication d'un solveur à réseau de tenseurs croît linéairement avec le temps, de sorte que la dimension de liaison nécessaire pour une précision fixe croît exponentiellement.

En une dimension, les excitations élémentaires sont les spinons, qui portent un spin et aucune charge, et les holons, qui portent une charge et aucun spin. Les deux se propagent à des vitesses différentes, de sorte qu'une perturbation localisée se scinde en un front rapide et un front lent. Cet effet s'appelle la séparation spin-charge, et ce tutoriel la mesure [1].

Ce tutoriel utilise la Qiskit Function Dynamics Simulator, qui prend le hamiltonien, un état initial, un calendrier de Trotter et une liste d'observables, et renvoie des valeurs moyennes. Elle gère en interne la correspondance de Jordan-Wigner, la trotterisation, la disposition des qubits et la suppression des erreurs, donc il n'y a aucun circuit à construire ni à transpiler dans ce tutoriel.

Exigences

Avant de commencer ce tutoriel, assure-toi d'avoir installé les paquets suivants :

• Qiskit SDK v2.3.0 ou ultérieur (pip install qiskit)

• Qiskit Runtime v0.50.0 ou ultérieur (pip install qiskit-ibm-runtime)

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

• NumPy v2.5.1 ou ultérieur (pip install numpy)

• Matplotlib v3.10.3 ou ultérieur (pip install matplotlib)

• Q-CTRL Visualizer v11.0.0 ou ultérieur (pip install qctrl-visualizer)

Configuration

Commence par t'authentifier avec ta clé API IBM Quantum. Sélectionne ensuite la Qiskit Function comme suit. (Ce code suppose que tu as déjà enregistré ton compte dans ton environnement 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")

Étape 1 : Associer les entrées classiques à un problème quantique

Définir le hamiltonien

La fonction prend le modèle sous forme de dictionnaire avec trois clés. Chacune accepte soit un seul flottant, qui s'applique uniformément le long de la chaîne, soit une liste qui précise la valeur site par site. Les listes ont $L-1$ entrées pour hopping_amplitude, car une chaîne ouverte a une liaison de moins que de sites, et $L$ entrées pour on_site_interaction et chemical_potential.

Ce tutoriel utilise une chaîne uniforme de 31 sites avec une interaction répulsive $U/t_h = 8$. Un potentiel chimique uniforme ne contribue qu'à une phase globale lorsque le nombre de particules est fixé, donc cet exemple le met à zéro. Le potentiel chimique n'a d'importance que lorsqu'il varie d'un site à l'autre, par exemple pour ajouter du désordre ou une impureté unique.

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

Définir le calendrier de Trotter

Le dictionnaire simulation contient le pas de temps ($\Delta t$) et le nombre de pas de Trotter, et le temps d'évolution total est leur produit. Des pas plus grands atteignent des temps plus tardifs à profondeur de circuit fixe, mais entraînent plus d'erreur de Trotter, qui croît avec $U$ parce que les termes du hamiltonien commutent moins. L'étude matérielle publiée a utilisé $\Delta t = 0.15$ pour $U/t_h = 8$ [2], et ce tutoriel utilise la même valeur. La profondeur du circuit croît linéairement avec le nombre de pas et ne dépend pas du tout 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)

Préparer l'état initial

Chaque site contient deux orbitales de spin, donc une chaîne de $L$ sites utilise $2L$ qubits. L'état initial est une chaîne d'occupation de Fock avec un caractère par site :

$$\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 chaîne est lue en commençant par le site de poids le plus faible, donc le caractère le plus à droite est le site 0 et l'indice physique du site du caractère $k$ compté depuis la gauche est $i = L - 1 - k$. La chaîne udud signifie donc que le site 0 est d, le site 1 est u, le site 2 est d et le site 3 est u.

Prépare ensuite un état de Néel à demi-remplissage, en alternant u et d, avec le site central remplacé par une lacune. La lacune est la perturbation dont on suit la propagation.

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

Construire les observables

Chaque observable est un dictionnaire qui représente une somme pondérée de produits d'opérateurs d'occupation :

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

Chaque entrée de operators est un terme $T_{\alpha}$, écrit sous forme de liste de paires [name, site] qui se multiplient, et l'entrée correspondante de coefficients est son poids. Les deux noms sont "Nup" pour $n_{j,\uparrow}$ et "Ndn" pour $n_{j,\downarrow}$. Une seule observable à un point $n_{3,\uparrow}$ s'écrit comme suit :

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

La cellule suivante affiche le plus petit ensemble complet d'entrées que la fonction accepte, ce qui sert de référence de forme pendant la lecture du code qui suit.

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

Trois groupes d'observables sont nécessaires. La densité de charge $n_i = n_{i\uparrow} + n_{i\downarrow}$ et le spin local $S^z_i = (n_{i\uparrow} - n_{i\downarrow})/2$ sont deux sommes à deux termes sur les mêmes opérateurs, qui ne diffèrent que par leurs coefficients. Le corrélateur de spin $S^z_i S^z_{i_*}$ entre le site $i$ et le site de la lacune se développe en quatre termes à deux points avec des coefficients $\pm 1/4$.

Les opérateurs d'occupation commutent, donc l'ordre des facteurs dans un terme est sans importance, mais l'ordre de la liste observables fixe l'ordre des valeurs renvoyées. Le code enregistre un indice pour chaque groupe pendant qu'il construit la liste. Le corrélateur ignore le site de la lacune, où le produit répéterait un facteur.

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

Étape 2 : Optimiser le problème pour l'exécution sur matériel quantique

Ce tutoriel n'a pas d'étape de transpilation. La fonction effectue en interne les étapes suivantes [2], et tu ne peux pas les configurer depuis le notebook :

• Associe les modes fermioniques aux qubits selon un ordre choisi pour garder de faibles poids de Pauli sur un réseau heavy-hex

• Utilise un réseau de SWAP fermioniques pour que chaque terme de saut devienne une rotation à deux qubits

• Plonge la chaîne le long d'un chemin qui évite les qubits faibles

• Applique le découplage dynamique et le Pauli twirling pendant la compilation

Tu choisis le Backend et le nombre de shots. La fonction prend un nom de Backend sous forme de chaîne, donc ce tutoriel sélectionne l'appareil le moins occupé suffisamment grand pour la chaîne et passe son nom.

Le plus petit signal que tu dois résoudre fixe le nombre de shots. Le traceur de charge atteint des amplitudes proches de $0.15$, mais le traceur de spin connecté culmine autour de $0.04$, soit environ quatre fois moins. L'étude publiée a utilisé 20 000 shots par circuit pour cette raison [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

Étape 3 : exécuter avec les primitives Qiskit

Soumets ensuite la simulation. run renvoie immédiatement un handle de job au lieu de bloquer.

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

Utilise job.status() pour suivre la progression et job.result() pour attendre la fin du job. Le résultat est un mapping dont l'entrée observable_expectation_values est indexée d'abord par observable, dans l'ordre où ils ont été fournis, puis par temps. Il y a step_count + 1 points temporels car la valeur à $t=0$ est incluse.

result = job.result()

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

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

Étape 4 : post-traiter et renvoyer le résultat dans le format classique souhaité

Reconstruire les corrélateurs

Redécoupe le tableau plat en trois groupes à l'aide des indices enregistrés plus tôt, puis forme les deux corrélateurs traceurs. Le traceur de charge est la variation de la densité par site par rapport à $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.$$

Le traceur de spin est le corrélateur connecté entre chaque site et le site de la lacune $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).$$

Les deux sont construits à partir d'opérateurs en base de Fock, ils proviennent donc des mêmes bitstrings mesurés et le second ne coûte aucune exécution matérielle supplémentaire.

Lors de la construction des observables, le code a enregistré la position de chaque opérateur de charge (charge_indices), de spin (spin_indices) et de corrélation de spin (correlator_indices) dans la liste. Comme le simulateur renvoie les valeurs d'espérance dans le même ordre, tu peux utiliser ces positions pour extraire les résultats correspondants. Pour les corrélations de spin, le code utilise aussi correlator_sites pour placer chaque résultat à la bonne position du réseau, et laisse la lacune à 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

Visualiser les fronts d'onde

Si tu traces chaque corrélateur sous forme de carte thermique de l'espace et du temps, le cône de lumière est directement visible. Le temps s'écoule vers le bas, et l'axe horizontal est la position le long de la chaîne.

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

Extraire les vitesses des fronts d'onde

Transforme les cartes thermiques en nombres en suivant le protocole de détection de l'étude publiée [2]. À chaque tranche de temps, symétrise le profil spatial autour de la lacune, ce qui élimine en moyenne la partie antisymétrique du bruit matériel, normalise par la plus grande amplitude hors centre, et note le site le plus éloigné au-dessus d'un seuil. Les détections sont ensuite découpées en séquences de pas de temps consécutifs, et seule la plus longue séquence est ajustée, ce qui écarte du résultat le transitoire des premiers instants.

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 vitesse ajustée dépend du seuil de détection et du quantile de rognage, et cette sensibilité domine le bruit de shots. Balaie les deux et indique la médiane avec l'enveloppe complète, comme le fait l'analyse publiée.

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

Redessiner les détections sur les cartes thermiques montre si l'algorithme a suivi le front avant. Les marqueurs pleins sont les détections qui sont entrées dans l'ajustement. Les marqueurs vides ont été écartés par les règles de masquage, soit parce qu'ils sont trop près de la lacune, trop près de l'extrémité de la chaîne, soit parce qu'ils sont assez tardifs pour que le front soit élargi et que sa position ne soit plus nette. La ligne en pointillés est l'ajustement sur les temps utilisés, et sa pente est la vitesse indiquée dans le titre de chaque panneau.

La colonne de la lacune est marquée sous l'axe. La charge y est tronquée car sa valeur est de l'ordre de un, et le spin n'a pas de valeur à cet endroit, il est donc dessiné en gris.

Le cône de charge s'ouvre plus largement que le cône de spin. Les deux perturbations partent du même défaut au même instant et se propagent vers l'extérieur à des vitesses différentes.

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

Le front de charge dépasse le front de spin. Un trou dans un fond de Mott unidimensionnel ne se propage pas comme un électron. Il se fractionne plutôt en un holon qui porte la charge et un spinon qui porte la perturbation de spin, et à $U/t_h = 8$ le spinon est le plus lent des deux.

L'ordre dépend de l'interaction. À $U = 0$, le modèle est celui de fermions libres avec $v_c = 2$ et $v_s = 4$, donc le spin est plus rapide. La répulsion supprime le transport de spin tout en laissant le transport de charge presque intact, et à $U/t_h = 8$ les valeurs matérielles publiées sont $v_c = 1.78$ et $v_s = 1.46$ [2]. Relancer avec plusieurs valeurs de $U$ et tracer $v_s / v_c$ cartographie la transition.

Note ces réserves sur les chiffres. La détection par seuil est légèrement biaisée vers le bas, donc une vitesse de charge inférieure de quelques pour cent à la valeur publiée est attendue. L'enveloppe du spin est large car le corrélateur connecté est plus proche du plancher de bruit, c'est pourquoi on indique la médiane du balayage plutôt qu'un seul ajustement.

Références

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