Simula la dinamica di Fermi-Hubbard 2D con i flow set

Importante

I concetti di questa guida sono attualmente disponibili solo nell'API Python. Una funzionalità equivalente sarà resa disponibile nell'API C in una versione futura.

Vedi anche

Leggi prima la guida ai flow set 1D. Questa guida ne è il seguito bidimensionale e non ripete le informazioni sviluppate lì: operatori di trasferimento e di vertice, che cos'è un flow set, come i groups sopravvivono nel circuito e come una codifica personalizzata viene collegata al Transpiler tramite MapperFnEvolutionSynthesis.

Consulta la guida al transpile per le fasi del Transpiler citate di seguito, e la guida ai mapper per la ricetta generale per scrivere un mapper personalizzato.

La guida 1D si concludeva ricordando che servono codifiche locali quando si lavora in due dimensioni. Questa guida compie quel passo, per il modello di Fermi-Hubbard bidimensionale su un reticolo quadrato, usando la codifica di Verstraete-Cirac (VC) [2].

Ci sono tre differenze rispetto al problema 1D, che formano la base di questo lavoro:

• I flow set si moltiplicano. Una catena ha due orientamenti delle frecce; un reticolo quadrato ne ha quattro (est, ovest, nord, sud), e ciascuno è un'unione di cammini orientati disgiunti nei vertici (uno per riga o colonna del reticolo) anziché un singolo cammino.

• La codifica guadagna un'ancilla per sito, non una per l'intero sistema, dando un

  rapporto qubit-modi di $2$. Gli operatori di trasferimento diventano Pauli di peso 3 (orizzontali) e di peso 4 (verticali), e ciascuno dei quattro flow set si mappa nella propria forma di Pauli, quindi ciascuno richiede il proprio Clifford diagonalizzante.

• L'Hamiltoniano codificato è fedele solo su un sottospazio: l'autospazio congiunto $+1$

  dei vincoli di plaquette che ogni codifica locale 2D porta con sé [2] [3]. In 1D la codifica era una semplice isometria e la verifica era un controllo dello spettro. Qui il confronto deve essere proiettato, come illustra il passaggio 9.

Il vantaggio ha la stessa forma che in 1D: la profondità a due qubit di un passo di Trotter diventa costante rispetto alla dimensione del reticolo, a $24$ ($16$ per l'hopping, in accordo con il rif. [1], più $8$ per l'interazione), mentre una trotterizzazione termine per termine dello stesso Hamiltoniano codificato cresce linearmente. Il numero di gate cresce comunque con il numero di legami, naturalmente; ciò che diventa indipendente dalla dimensione è quanti strati occupano quei gate. Le profondità sono sempre contate su un registro astratto completamente connesso; mappare la codifica sulla coupling map limitata di un dispositivo è un problema a sé, discusso alla fine.

Il reticolo e le sue quattro orientazioni

Il primo cambiamento è visibile sul grafo di interazione orientato. I siti sono numerati in ordine row-major, site(r, c) = r * cols + c, e ogni legame porta entrambi gli orientamenti, come la catena in 1D, ma ora i legami corrono lungo due assi:

La regola di commutazione è quella 1D, invariata: due operatori di trasferimento che condividono un sito commutano quando le frecce vi fluiscono attraverso e anticommutano quando si scontrano. Tuttavia, un sito ha ora fino a quattro legami incidenti anziché due, quindi le frecce che lo attraversano possono scontrarsi in più modi, e un insieme di frecce che non si scontrano mai non è più un singolo cammino. Leggere le frecce orizzontali da sinistra a destra dà un insieme di questo tipo, le frecce verticali dall'alto in basso ne danno un altro, e invertendo ciascuno se ne ottengono altri due. Questi sono i flow set est, ovest, nord e sud che organizzano tutto ciò che segue.

La codifica di Verstraete-Cirac

Verstraete-Cirac (VC) associa a ogni modo fermionico $j$ il proprio qubit ausiliario. Con $N$ siti del reticolo, il modo $j$ vive sul qubit $j$ e la sua ancilla sul qubit $N + j$, quindi il registro contiene $2N$ qubit. Gli operatori di vertice restano il più semplici possibile,

$$V_j = Z_j \, ,$$

il che rende le occupazioni direttamente leggibili più avanti. Gli operatori di trasferimento dipendono dall'orientamento dell'arco. Per un arco orizzontale $j \to k$ con $j$ il sito occidentale, e un arco verticale con $j$ il sito settentrionale,

$$T_{jk}^{\text{horizontal}} = \tfrac{1}{2} X_j X_k Z_{N+j} \, , \qquad T_{jk}^{\text{vertical}} = \tfrac{1}{2} X_j Y_k X_{N+j} Y_{N+k} \, ,$$

un “triangolo” di peso 3 e un “quadrato” di peso 4. Identificare ciascuna freccia è importante: con $j$ come sito occidentale, $j \to k$ punta a est, ma con $j$ come sito settentrionale punta a sud. Gli orientamenti inversi seguono dall'identità $T_{kj} = -V_j V_k T_{jk}$, che trasforma il triangolo XXZ in YYZ, e il quadrato XYXY in YXXY. Quelle quattro etichette distinte sono i quattro flow set, e il resto di questa guida li segue.

Nota

Queste non sono esattamente l'Eq. (A1) di Ref. [1], che afferma “fino a rotazioni a singolo qubit” e differisce in tre aspetti che vale la pena nominare, dato che chi controlla l'articolo se ne accorgerà:

• Il prefattore: l'Eq. (A1) non ne ha alcuno. Il $\tfrac12$ qui presente è ciò che rende $T_{jk}^2 = \tfrac14$, in accordo con la normalizzazione $T_{jk} = -\tfrac12 (a^\dagger_j - a_j)(a^\dagger_k + a_k)$ che usa TransferVertexOperator. L'articolo ha bisogno degli operatori di trasferimento solo all'interno di $\exp(-i\,dt\,H)$, dove un fattore di scala viene assorbito nell'angolo di rotazione; qui devono invece essere operatori che soddisfano un'identità.

• L'orientamento: l'Eq. (A1) seleziona le lettere in base alla parità di riga di $j$, che per un lato orizzontale è la stessa per entrambi gli estremi, quindi non distingue $T_{jk}$ da $T_{kj}$. L'articolo recupera l'orientamento separatamente, coniugando con uno strato di $R^Z(\pi/2)$. Questa guida invece lo incorpora nella codifica, ricavando l'inverso da $T_{kj} = -V_j V_k T_{jk}$, ed è questo che fa sì che i due orientamenti finiscano su Pauli diversi nel passo 3.

• La forma verticale: l'Eq. (A1) alterna XXXX e YYYY in base alla parità di riga, mentre questa guida usa un unico XYXY. Sono due modi per ottenere la stessa cosa: un'asimmetria di lettere che fa sì che un lato verticale entri in conflitto correttamente con i lati orizzontali che lo incontrano. Non interpretare la differenza come una libera scelta di base. Eliminare l'alternanza di parità e usare un unico XXXX rompe del tutto l'algebra, ed è per questo che le lettere distinte vengono verificate anziché date per scontate.

Nota

L'ancilla è associata a un sito, dando il rapporto $2$; è la codifica Derby-Klassen a collocare un qubit ausiliario per ogni plaquette (dispari), ottenendo il rapporto minore $3/2$. Qui si usa VC perché i suoi operatori di trasferimento sono quelli di cui i Clifford diagonalizzanti sono forniti in forma chiusa da Ref. [1].

1. L'Hamiltoniana su un reticolo quadrato

La costruzione è la stessa dell'1D, ma con una lista di lati bidimensionale. I siti mantengono la numerazione per righe usata per la figura precedente, site(r, c) = r * cols + c, così che i legami orizzontali collegano indici consecutivi e i legami verticali collegano indici che differiscono di cols:

>>> from collections import defaultdict
>>> from qiskit_fermions.operators import TransferVertexOperator
>>>
>>> def fermi_hubbard_2d(rows, cols, tunneling, interaction):
...     site = lambda r, c: r * cols + c
...     edges = [(site(r, c), site(r, c + 1)) for r in range(rows) for c in range(cols - 1)]
...     edges += [(site(r, c), site(r + 1, c)) for r in range(rows - 1) for c in range(cols)]
...     data: defaultdict[tuple[tuple[int, int], ...], complex] = defaultdict(complex)
...     for j, k in edges:
...         # hopping: -t (a^dag_j a_k + a^dag_k a_j) = t (T_jk + T_kj)
...         data[((j, k),)] += tunneling
...         data[((k, j),)] += tunneling
...         # interaction: U n_j n_k with n_j = (1 - V_j) / 2
...         data[()] += interaction / 4
...         data[((j, j),)] -= interaction / 4
...         data[((k, k),)] -= interaction / 4
...         data[((j, j), (k, k))] += interaction / 4
...     return TransferVertexOperator.from_dict(data)
>>>
>>> rows, cols = 3, 3
>>> num_sites = rows * cols
>>> hamiltonian = fermi_hubbard_2d(rows, cols, tunneling=1.0, interaction=2.0)
>>> print(f"Number of terms: {len(hamiltonian)}")
Number of terms: 46

2. Quattro insiemi di flusso invece di due

Ogni termine di hopping viene classificato in base alla direzione in cui punta la sua freccia sul reticolo. Confrontando le coordinate riga/colonna divmod dei due estremi si distinguono i lati orizzontali da quelli verticali, e il segno della differenza determina l'orientamento. Tutto ciò che è diagonale (l'identità, i termini $V_j$ e $V_j V_k$) confluisce in un quinto gruppo:

Nota

Convenzione di orientamento. I nomi sono semplicemente ricavati dalla figura precedente: una freccia è est o ovest a seconda che punti a destra o a sinistra lungo una riga, e nord o sud a seconda che punti in alto o in basso lungo una colonna. Poiché la numerazione procede per righe, est e sud sono le direzioni in cui l'indice aumenta. Nulla a valle dipende dai nomi (servono solo a etichettare i gruppi), ma è facile confonderli, quindi conviene fissarli rispetto alla figura. Ref. [1] assegna gli stessi nomi ai suoi insiemi di flusso su reticolo quadrato.

>>> EAST, WEST, NORTH, SOUTH, DIAGONAL = range(5)
>>>
>>> def flow_set_groups(operator, cols):
...     groups = []
...     for terms, _ in operator.iter_terms():
...         if len(terms) == 1 and terms[0][0] != terms[0][1]:
...             j, k = terms[0]
...             row_j, col_j = divmod(j, cols)
...             row_k, col_k = divmod(k, cols)
...             if row_j == row_k:
...                 groups.append(EAST if col_k > col_j else WEST)
...             else:
...                 groups.append(SOUTH if row_k > row_j else NORTH)
...         else:
...             groups.append(DIAGONAL)
...     return groups
>>>
>>> hamiltonian.groups = flow_set_groups(hamiltonian, cols)
>>>
>>> names = ["east", "west", "north", "south", "diagonal"]
>>> for index, name in enumerate(names):
...     count = sum(1 for group in hamiltonian.groups if group == index)
...     print(f"{name:9s} {count:2d}")
east       6
west       6
north      6
south      6
diagonal  22

All'interno di un insieme di flusso ogni freccia o ne incontra un'altra testa-coda oppure non la incontra affatto, come in 1D, quindi un insieme di flusso continua a commutare termine per termine e la sua evoluzione non comporta errore di Trotter. La novità è che un insieme di flusso ora è sconnesso: l'insieme est di un reticolo rows x cols è formato da rows cammini indipendenti da sinistra a destra. Questa sconnessione è un vantaggio più che una complicazione, come mostra il passo 6.

3. Codificare l'operatore

La codifica è una funzione che associa un singolo operatore di trasferimento generalizzato a una stringa di Pauli, che viene passata a map_transfer_vertex_generators(). Rispetto alla versione 1D, il corpo ha più casi: diagonale, orizzontale, verticale e gli orientamenti inversi tramite $T_{kj} = -V_j V_k T_{jk}$.

C'è un problema. map_transfer_vertex_generators() richiama la callback con un generatore e nient'altro, mentre per decidere la stringa di Pauli di un legame serve anche cols, per ricavare riga e colonna a partire dall'indice di un modo. Invece di annidare le funzioni per chiuderle su di esso, lo vincola functools.partial(), il che mantiene entrambe le funzioni a livello superiore e richiamabili in modo indipendente.

>>> from functools import partial
>>>
>>> from qiskit.quantum_info import SparseObservable
>>> from qiskit_fermions.mappers import map_transfer_vertex_generators
>>>
>>> def ancilla_qubit(mode, num_sites):
...     """The auxiliary qubit paired with fermionic mode `mode`."""
...     return num_sites + mode
>>>
>>> def vc_action(generator, width, cols):
...     """Map one generalized transfer operator to its Verstraete-Cirac Pauli."""
...     num_sites = width // 2
...     ancilla = partial(ancilla_qubit, num_sites=num_sites)
...     j, k = generator
...     if j == k:  # vertex operator V_j = Z_j
...         return SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     row_j, col_j = divmod(j, cols)
...     row_k, col_k = divmod(k, cols)
...     if row_j == row_k:  # horizontal: weight-3 triangle
...         left, right = (j, k) if col_j < col_k else (k, j)
...         operator = SparseObservable.from_sparse_list(
...             [("XXZ", [left, right, ancilla(left)], 0.5)], num_qubits=width
...         )
...         forward = (j, k) == (left, right)
...     else:  # vertical: weight-4 square
...         up, down = (j, k) if row_j < row_k else (k, j)
...         operator = SparseObservable.from_sparse_list(
...             [("XYXY", [up, down, ancilla(up), ancilla(down)], 0.5)], num_qubits=width
...         )
...         forward = (j, k) == (up, down)
...     if forward:
...         return operator
...     # reverse orientation: T_kj = -V_j V_k T_jk
...     v_j = SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     v_k = SparseObservable.from_sparse_list([("Z", [k], 1.0)], num_qubits=width)
...     return (-1.0 * v_j.compose(v_k).compose(operator)).simplify()
>>>
>>> def vc_encoding(operator, width, cols):
...     """Map a TransferVertexOperator to a qubit observable, Verstraete-Cirac style."""
...     return map_transfer_vertex_generators(
...         operator,
...         partial(vc_action, width=width, cols=cols),
...         identity=lambda: SparseObservable.identity(width),
...         compose=SparseObservable.compose,
...     ).simplify()
>>>
>>> # bind the lattice width once; the transpiler only ever passes (operator, width)
>>> encoding = partial(vc_encoding, cols=cols)
>>> num_qubits = 2 * num_sites
>>> num_qubits
18

Il senso di raggruppare prima di mappare si vede subito: ogni insieme di flusso viene codificato in un'unica forma di Pauli, uniforme su tutto l'insieme.

>>> for index, name in enumerate(names):
...     observable = encoding(hamiltonian.split_out_groups()[index], num_qubits)
...     shapes = sorted({label for label, _, _ in observable.to_sparse_list()})
...     print(f"{name:9s} {shapes}")
east      ['XXZ']
west      ['YYZ']
north     ['YXXY']
south     ['XYXY']
diagonal  ['', 'Z', 'ZZ']

Le quattro forme di hopping sono distinte, quindi ogni insieme di flusso ha il proprio Clifford nel passo 5. Gli orientamenti inversi sono Pauli diversi, non gli stessi rietichettati; per questo la sintesi di un insieme di flusso 2D richiede quattro circuiti, mentre la sintesi 1D ne richiede uno.

4. Transpilare con la codifica personalizzata

Questa fase è quella 1D con un registro di qubit più ampio: CustomF2QLayout mappa il registro di $N$ modi su $2N$ qubit, MapperFnEvolutionSynthesis porta con sé la codifica, e il suo argomento product_formula seleziona come viene scomposto il PauliEvolutionGate risultante. Confronta il passo 6 della guida 1D, che costruisce la stessa pipeline per una singola ancilla.

>>> from qiskit.circuit import QuantumRegister
>>> from qiskit.passmanager import MultiStagePassManager
>>> from qiskit.synthesis import LieTrotter
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution
>>> from qiskit_fermions.circuit.library.synthesis import FermionicLieTrotter
>>> from qiskit_fermions.transpiler import FermionicCircuitToDAG, QuantumDAGToCircuit
>>> from qiskit_fermions.transpiler.passes import (
...     CustomF2QLayout,
...     F2QSynthesis,
...     MapperFnEvolutionSynthesis,
... )
>>>
>>> total_time = 1.0
>>>
>>> def flow_set_pass_manager(mode_register, encoding, product_formula=None):
...     """Build the pipeline mapping `mode_register` with the VC `encoding`."""
...     synthesis = F2QSynthesis()
...     synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(
...         encoding, product_formula=product_formula
...     )
...     # this encoding spends one ancilla per fermionic mode
...     qubit_register = QuantumRegister(2 * mode_register.size, "q")
...     return MultiStagePassManager(
...         input=FermionicCircuitToDAG(),
...         layout=CustomF2QLayout({mode_register: qubit_register}),
...         synthesis=synthesis,
...         output=QuantumDAGToCircuit(),
...     )
>>>
>>> circuit = FermionicCircuit(num_sites)
>>> circuit.append(
...     Evolution(
...         num_sites, hamiltonian, time=total_time, synthesis=FermionicLieTrotter()
...     ),
...     circuit.modes,
... )
>>>
>>> # as in 1D, one decompose() splits the Evolution group by group (five gates here,
>>> # the four hopping flow sets in group order, then the diagonal interaction)
>>> flow_set_circuit = circuit.decompose()
>>> flow_set_circuit.count_ops()["Evolution"]
5

Cinque Gate, che portano i 6, 6, 6, 6 e 22 termini contati nel passo 2. Ciascuno viene mappato e sintetizzato in modo indipendente, quindi la partizione decisa a livello fermionico sopravvive fino al circuito:

>>> pass_manager = flow_set_pass_manager(flow_set_circuit.register, encoding, LieTrotter())
>>> qubit_circuit = pass_manager.run(flow_set_circuit)
>>> qubit_circuit.num_qubits
18
>>>
>>> def two_qubit_stats(circuit):
...     decomposed = circuit.decompose(reps=6)
...     two_qubit = lambda instruction: len(instruction.qubits) == 2
...     return decomposed.count_ops().get("cx", 0), decomposed.depth(two_qubit)
>>>
>>> two_qubit_stats(qubit_circuit)
(144, 47)

Questa è la baseline: la codifica e il raggruppamento sono entrambi presenti, ma LieTrotter sintetizza ogni insieme di flusso termine per termine, e la profondità cresce con il reticolo.

5. Un Clifford di profondità 2 per insieme di flusso

Il meccanismo è lo stesso dell'1D: coniugare un intero insieme di flusso commutante con un Clifford che mappa ciascuno dei suoi termini su un Pauli distinto di peso 1, in modo che l'evoluzione dell'intero insieme costi solo rotazioni a singolo qubit. In 1D lo faceva un'unica struttura a mattoni di Gate CZ.

In 2D ogni insieme di flusso richiede un proprio circuito, seguendo la Sezione IV C e la Fig. 7 di Ref. [1]. Tutti e quattro agiscono su una linea del reticolo alla volta (una riga per est/ovest, una colonna per nord/sud) insieme alle ancille di quella linea, e tutti e quattro hanno la stessa forma in tre parti:

• Uno strato di rotazioni a singolo qubit, che pone ciascun qubit nella base richiesta dal suo insieme di flusso.

• Un primo strato di CX, sempre da sito ad ancilla sullo stesso sito.

• Un secondo strato di CX, che accoppia i vicini lungo la linea.

I passi 2 e 3 sono gli unici due strati a due qubit, quindi la profondità a due qubit è $2$, indipendentemente dalla lunghezza della linea. Sono scritti qui sotto, una funzione per insieme di flusso. Le rotazioni sono quelle che sono perché sono state trovate coniugando i termini di ciascun insieme di flusso e richiedendo immagini distinte di peso 1; quindi considerale date e affidati alla verifica che segue.

I due insiemi orizzontali sono i più semplici. I loro termini hanno peso 3, XXZ e YYZ, e differiscono solo per l'$S$ aggiuntivo che trasforma $X$ in $Y$ sui siti:

>>> from qiskit.circuit import QuantumCircuit
>>>
>>> def east_clifford(circuit, line, ancilla):
...     """Diagonalize the east flow set (XXZ) on one lattice row."""
...     for j in line:
...         circuit.h(ancilla(j))
...     for j in line[:-1]:  # layer 1: each site with its own ancilla
...         circuit.cx(ancilla(j), j)
...     for j, k in zip(line, line[1:]):  # layer 2: ancilla with the next site along
...         circuit.cx(k, ancilla(j))
>>>
>>> def west_clifford(circuit, line, ancilla):
...     """Diagonalize the west flow set (YYZ) on one lattice row."""
...     for j in line:
...         circuit.s(j)  # the only difference from east: X -> Y on the sites
...         circuit.h(ancilla(j))
...     for j in line[:-1]:
...         circuit.cx(ancilla(j), j)
...     for j, k in zip(line, line[1:]):
...         circuit.cx(k, ancilla(j))

I due insiemi verticali hanno peso 4, XYXY e YXXY, e sono corrispondentemente più elaborati: le rotazioni si alternano con la posizione lungo la linea, e anche il secondo strato di CX si alterna, accoppiando i siti nelle posizioni pari e le ancille in quelle dispari (o viceversa):

>>> def south_clifford(circuit, line, ancilla):
...     """Diagonalize the south flow set (XYXY) on one lattice column."""
...     for position, j in enumerate(line):
...         if position % 2 == 0:
...             circuit.sdg(ancilla(j))
...             circuit.h(ancilla(j))
...         else:
...             circuit.h(j)
...             circuit.s(j)
...             circuit.h(ancilla(j))
...             circuit.s(ancilla(j))
...     for j in line:  # layer 1
...         circuit.cx(ancilla(j), j)
...     for position, (j, k) in enumerate(zip(line, line[1:])):  # layer 2
...         if position % 2 == 0:
...             circuit.cx(ancilla(j), ancilla(k))  # ancilla to ancilla
...         else:
...             circuit.cx(j, k)  # site to site
>>>
>>> def north_clifford(circuit, line, ancilla):
...     """Diagonalize the north flow set (YXXY) on one lattice column."""
...     for position, j in enumerate(line):
...         if position % 2 == 0:
...             circuit.h(ancilla(j))
...             circuit.s(ancilla(j))
...         else:
...             circuit.h(j)
...             circuit.s(j)
...             circuit.sdg(ancilla(j))
...             circuit.h(ancilla(j))
...     for j in line:  # layer 1
...         circuit.cx(ancilla(j), j)
...     for position, (j, k) in enumerate(zip(line, line[1:])):  # layer 2
...         if position % 2 == 0:
...             circuit.cx(k, j)  # site to site, and the other way round from south
...         else:
...             circuit.cx(ancilla(j), ancilla(k))

Disegnato su una singola linea di tre siti (prima i siti, poi le loro ancille) l'intero meccanismo sta in un'unica figura. La forma a triangolo è la più chiara delle due: un cambio di base sulle ancille, poi due strati di CX, il primo accoppia ogni sito con la propria ancilla e il secondo raggiunge il sito successivo lungo la linea.

>>> line = [0, 1, 2]  # three sites, with ancillas on qubits 3, 4, 5
>>> triangle = QuantumCircuit(6)
>>> east_clifford(triangle, line, lambda mode: 3 + mode)
>>> triangle.draw("mpl")
<Figure size ... with 1 Axes>

La forma a quadrato richiede più lavoro a singolo qubit (le rotazioni si alternano lungo la linea) e il suo secondo strato di CX si alterna tra l'accoppiare siti e l'accoppiare ancille:

>>> square = QuantumCircuit(6)
>>> south_clifford(square, line, lambda mode: 3 + mode)
>>> square.draw("mpl")
<Figure size ... with 1 Axes>

Entrambe hanno profondità di due strati di CX, e aggiungere siti alla linea le allarga senza approfondirle; questa è l'intera affermazione di profondità costante, visibile a occhio. I Clifford ovest e nord differiscono solo per gli strati a singolo qubit, quindi non vengono disegnati separatamente.

Ciò di cui la sintesi ha ancora bisogno è determinare dove finisce ogni termine. Una linea di $L$ siti ha $L - 1$ legami, ciascuno dei quali contribuisce con un termine, quindi la funzione di supporto qui sotto restituisce una voce per legame: il qubit su cui il termine è finito e il Pauli a singolo qubit in cui si è trasformato. Fondamentale: questi valori sono letti dalla coniugazione, non previsti. A differenza dell'1D, dove una struttura a mattoni di CZ mappa $X_{j+1} \mapsto Z_j X_{j+1} Z_{j+2}$ in forma chiusa, qui non esiste una formula simile, e le ancille di una linea di colonna non sono contigue. L'unica cosa che differisce tra gli insiemi di flusso è la stringa di Pauli in ingresso, che il passo 3 ha già stampato:

>>> from qiskit.quantum_info import Clifford, SparsePauliOp
>>>
>>> def clifford_supports(clifford, pauli, line, ancilla, num_qubits):
...     """Conjugate a flow set's terms to find which qubit each one lands on.
...
...     `pauli` is the flow set's Pauli string, as printed in step 3. Its letters apply
...     to the i-th term's qubits in the order below: the two sites of the bond, then
...     their ancillas.
...
...     This is the one check that needs SparsePauliOp: Clifford conjugation wants a
...     `Pauli` and its symplectic bits, and SparseObservable exposes neither -- its own
...     `evolve` conjugates by a Pauli, not by a Clifford.
...     """
...     supports = []
...     for i in range(len(line) - 1):
...         qubits = [line[i], line[i + 1], ancilla(line[i]), ancilla(line[i + 1])]
...         term = SparsePauliOp.from_sparse_list(
...             [(pauli, qubits[: len(pauli)], 1.0)], num_qubits=num_qubits
...         )
...         image = term.paulis[0].evolve(clifford, frame="s")
...         # a weight-1 image has exactly one qubit with an X and/or Z bit set
...         (support,) = np.nonzero(image.x | image.z)
...         assert len(support) == 1, f"term {i} is not weight-1"
...         supports.append((int(support[0]), image[support].to_label()))
...     return supports

Le immagini devono risultare distinte e di peso 1, in modo che ogni termine possa essere ruotato indipendentemente sul proprio qubit; la profondità a due qubit di $2$ è già visibile nelle quattro costruzioni precedenti, ciascuna delle quali emette due strati di CX a prescindere dalla lunghezza della linea. Di seguito, entrambe le proprietà vengono verificate su ogni insieme di flusso, ogni linea e diverse dimensioni del reticolo. Il passo 6 chiama poi clifford_supports, quindi è una parte operativa della sintesi e non solo una dimostrazione.

6. Sintesi consapevole degli insiemi di flusso

La sintesi parte avvantaggiata. Poiché Evolution si divide da sé gruppo per gruppo, synthesize viene chiamata una volta per insieme di flusso (cinque volte qui, quattro insiemi di hopping più il gruppo diagonale), così ogni termine di hopping che vede condivide un'unica forma di Pauli, e un solo Clifford li diagonalizza tutti. Deve solo leggere quella forma dal primo termine e ordinare gli altri in base alla linea su cui si trovano. Segue poi l'osservazione che rende le due dimensioni più economiche di quanto si potrebbe temere:

Importante

Le componenti di un singolo insieme di flusso sono disgiunte sui vertici: le righe dell'insieme ovest non condividono alcun qubit. Quindi tutti i Clifford di un insieme di flusso possono essere emessi insieme e condividono gli stessi due strati a due qubit, anziché mettersi in coda una riga dopo l'altra. Un insieme di flusso costa quindi un numero fisso di strati (due per il Clifford e due per il suo inverso) indipendentemente dal numero di righe o colonne del reticolo, ed è per questo che la profondità totale risulta costante.

>>> from qiskit.synthesis.evolution import EvolutionSynthesis
>>>
>>> # The interaction couples nearest neighbors only, so every interior site carries four
>>> # Rzz gates -- one per bond -- and can only do one at a time: four layers is the best
>>> # any schedule can achieve. Horizontal bonds split by column parity, vertical ones by
>>> # row parity, which achieves it.
>>> def rzz_layers(circuit, angles, cols):
...     """Emit the diagonal Rzz gates in four layers of disjoint qubit pairs."""
...     horizontal = [(j, k) for j, k in angles if k - j == 1]
...     vertical = [(j, k) for j, k in angles if k - j == cols]
...     schedule = [
...         [(j, k) for j, k in horizontal if (j % cols) % 2 == 0],
...         [(j, k) for j, k in horizontal if (j % cols) % 2 == 1],
...         [(j, k) for j, k in vertical if (j // cols) % 2 == 0],
...         [(j, k) for j, k in vertical if (j // cols) % 2 == 1],
...     ]
...     # emit layer by layer, since Qiskit's depth follows the order gates were added
...     for layer in schedule:
...         for j, k in layer:
...             circuit.rzz(angles[j, k], j, k)
>>>
>>> class FlowSetSynthesis(EvolutionSynthesis):
...     """Synthesize one 2D line flow set at a time, via its depth-2 Clifford."""
...
...     preserve_order = True
...
...     def __init__(self, rows, cols):
...         self.rows, self.cols = rows, cols
...
...     def rows_of_lattice(self):
...         return [[r * self.cols + c for c in range(self.cols)] for r in range(self.rows)]
...
...     def columns_of_lattice(self):
...         return [[r * self.cols + c for r in range(self.rows)] for c in range(self.cols)]
...
...     def synthesize(self, evolution):
...         num_qubits = 2 * self.rows * self.cols
...         ancilla = partial(ancilla_qubit, num_sites=self.rows * self.cols)
...         circuit = QuantumCircuit(evolution.operator.num_qubits)
...         terms = [
...             (label, indices, 2 * evolution.time * coeff.real)
...             for label, indices, coeff in evolution.operator.to_sparse_list()
...         ]
...
...         # `Evolution.decompose()` hands us one gate per flow set, so every hopping term
...         # here shares a single Pauli shape and one Clifford diagonalizes all of them
...         shapes = tuple({label for label, _, _ in terms} - {"", "Z", "ZZ"})
...         match shapes:
...             case ("XXZ",):
...                 clifford_function, lines = east_clifford, self.rows_of_lattice()
...                 along_rows = True
...             case ("YYZ",):
...                 clifford_function, lines = west_clifford, self.rows_of_lattice()
...                 along_rows = True
...             case ("XYXY",):
...                 clifford_function, lines = south_clifford, self.columns_of_lattice()
...                 along_rows = False
...             case ("YXXY",):
...                 clifford_function, lines = north_clifford, self.columns_of_lattice()
...                 along_rows = False
...             case ():  # the diagonal group: no conjugation needed
...                 self.diagonal(circuit, terms)
...                 return circuit
...             case _:
...                 raise ValueError(f"expected one known flow set per gate, got {shapes}")
...
...         # a term's line is the row it sits in and its slot the column, or vice versa
...         rotations: dict[int, list] = defaultdict(list)
...         for _, indices, angle in terms:
...             row, col = divmod(indices[0], self.cols)
...             line, slot = (row, col) if along_rows else (col, row)
...             rotations[line].append((slot, angle))
...
...         active = sorted(rotations)
...         # the lines of one flow set are vertex-disjoint, so all of their Cliffords
...         # are emitted together and share the same two two-qubit layers
...         conjugation = QuantumCircuit(num_qubits)
...         for line in active:
...             clifford_function(conjugation, lines[line], ancilla)
...         circuit.compose(conjugation, inplace=True)
...         for line in active:
...             supports = clifford_supports(
...                 Clifford(conjugation), shapes[0], lines[line], ancilla, num_qubits
...             )
...             for slot, angle in rotations[line]:
...                 match supports[slot]:
...                     case qubit, "X":
...                         circuit.rx(angle, qubit)
...                     case qubit, "Y":
...                         circuit.ry(angle, qubit)
...                     case qubit, "Z":
...                         circuit.rz(angle, qubit)
...         circuit.compose(conjugation.inverse(), inplace=True)
...         return circuit
...
...     def diagonal(self, circuit, terms):
...         """Synthesize the diagonal group: a phase, single-qubit Z and paired Rzz."""
...         rzz = {}
...         for label, indices, angle in terms:
...             match label, indices:
...                 case "", _:  # the identity contributes only a global phase
...                     circuit.global_phase -= angle / 2
...                 case "Z", [qubit]:  # a lone vertex operator
...                     circuit.rz(angle, qubit)
...                 case "ZZ", [j, k]:  # interaction, see rzz_layers above
...                     rzz[min(j, k), max(j, k)] = angle
...                 case _:
...                     raise ValueError(f"unexpected Pauli term: {label} on {indices}")
...         rzz_layers(circuit, rzz, self.cols)

Come in 1D, questa è scritta deliberatamente per questa codifica: riconosce le quattro forme di Pauli e solleva un errore per qualsiasi altra cosa, anziché ripiegare silenziosamente su altro. Presuppone inoltre la suddivisione un-insieme-di-flusso-per-Gate vista sopra, e lo dichiara; più di una forma di hopping in un singolo Gate significa che il raggruppamento non è stato applicato, il che è un errore e non qualcosa da aggirare.

Nota

Come la sua controparte 1D, questa sintesi non accetta un numero di passi di Trotter. Un passo di Trotter deve ciclare tra gli insiemi di flusso, ma transpile_trotter_step qui sotto scompone il circuito fermionico in un Gate Evolution per insieme, quindi questo metodo vede sempre un solo insieme alla volta e potrebbe ripetere solo al suo interno, il che non cambia nulla: ripetere un solo operatore divide semplicemente la sua evoluzione in parti uguali.

La Trotterizzazione a più passi richiede quindi una formula di prodotto che rispetti i groups dell'operatore all'interno di un singolo Gate, prima della scomposizione. FermionicLieTrotter fa proprio questo, e il passo 10 lo usa per misurare la convergenza.

Nota le due classi dal nome simile in gioco. LieTrotter di Qiskit, usata come baseline sopra, agisce sull'Hamiltoniana mappata, dove il raggruppamento è già perso. FermionicLieTrotter è la controparte fermionica, ed è quella che può ancora vedere gli insiemi di flusso.

7. Il risultato: profondità costante

Tutto fin qui è stato eseguito sul singolo reticolo $3 \times 3$ costruito nel passo 1. Confrontare i due percorsi di sintesi su reticoli crescenti richiede che l'intera sequenza (operatore, raggruppamento, codifica, circuito, pass manager) sia ripetuta per ogni dimensione, quindi conviene racchiuderla una volta per tutte:

>>> def transpile_trotter_step(rows, cols, product_formula, steps=1):
...     """Transpile `steps` Trotter steps on a `rows` x `cols` lattice."""
...     num_sites = rows * cols
...     operator = fermi_hubbard_2d(rows, cols, 1.0, 2.0)
...     operator.groups = flow_set_groups(operator, cols)
...     encoding = partial(vc_encoding, cols=cols)
...
...     circuit = FermionicCircuit(num_sites)
...     circuit.append(
...         Evolution(
...             num_sites,
...             operator,
...             time=total_time,
...             synthesis=FermionicLieTrotter(reps=steps),
...         ),
...         circuit.modes,
...     )
...     circuit = circuit.decompose()  # split into one Evolution per flow set
...
...     pass_manager = flow_set_pass_manager(circuit.register, encoding, product_formula)
...     return pass_manager.run(circuit)

Raccogliere entrambe le metriche su un intervallo di dimensioni del reticolo e tracciarle rispetto al lato del reticolo rende chiara la differenza:

>>> sizes = [3, 4, 5, 6, 7, 8]
>>> lie_trotter = [
...     two_qubit_stats(transpile_trotter_step(n, n, LieTrotter())) for n in sizes
... ]
>>> flow = [
...     two_qubit_stats(transpile_trotter_step(n, n, FlowSetSynthesis(n, n))) for n in sizes
... ]
>>>
>>> figure, axes = plt.subplots(1, 2, figsize=(9, 3.5), layout="constrained")
>>> for panel, index, title in zip(axes, (1, 0), ("two-qubit depth", "CX count")):
...     _ = panel.plot(sizes, [stat[index] for stat in lie_trotter], "o-", label="qubit-level")
...     _ = panel.plot(sizes, [stat[index] for stat in flow], "s-", label="flow set")
...     _ = panel.set(xlabel="lattice edge $L$ (an $L \\times L$ lattice)", title=title)
...     _ = panel.legend()
>>> figure
<Figure size ... with 2 Axes>

La profondità a due qubit è 24 e costante da $3 \times 3$ a $8 \times 8$ (cioè da 18 qubit fino a 128), mentre LieTrotter applicato all'Hamiltoniana codificata in modo identico cresce linearmente con il lato del reticolo, da 47 a 92. Il numero di Gate cresce in entrambi i casi, dato che ci sono comunque $O(N)$ termini da applicare; ciò che gli insiemi di flusso eliminano è la crescita del numero di strati che quei Gate occupano.

Il budget si divide in 16 per l'hopping e otto per l'interazione:

• Ciascuno dei quattro insiemi di flusso costa $2 + 2$ strati a due qubit: il suo Clifford in ingresso e l'inverso in uscita, entrambi di profondità due per costruzione, con solo rotazioni a singolo qubit nel mezzo. Fa $4 \times 4 = 16$, che è la cifra citata da Ref. [1], poiché conta solo i termini di hopping. Impostando $U = 0$ nel codice precedente la si riproduce.

• Il gruppo diagonale aggiunge i restanti otto. I suoi Gate Rzz si impacchettano in quattro strati, e ogni Rzz si scompone in due Gate CX attorno a una rotazione a singolo qubit, quindi $4 \times 2 = 8$. Quattro strati è l'ottimo: su un reticolo quadrato l'interazione raggiunge solo i primi vicini, quindi ogni sito interno porta quattro Gate Rzz (uno per legame) e può partecipare a uno solo alla volta.

Nessuna delle due somme fa riferimento alla dimensione del reticolo, che è proprio il punto.

8. Verificare il circuito

La sintesi per insiemi di flusso e LieTrotter descrivono lo stesso passo di Trotter della stessa Hamiltoniana codificata, quindi le loro unitarie devono coincidere a meno di una fase globale. I vettori di stato casuali rilevano qualsiasi discrepanza molto più economicamente che costruire l'unitaria completa:

Nota

Entrambi i circuiti vengono scomposti in Gate a uno e due qubit prima di essere applicati a uno stato. Se non scomposta, ogni istruzione Evolution di alto livello costringe Statevector a costruire il suo intero operatore denso $2^{2N} \times 2^{2N}$; scomporre prima gli permette invece di applicare i Gate uno alla volta, il che qui è centinaia di volte più veloce.

>>> import numpy as np
>>> from qiskit.quantum_info import Statevector
>>>
>>> rng = np.random.default_rng(1234)
>>>
>>> def agrees_up_to_one_phase(rows, cols, samples=3):
...     """Compare both syntheses on random states.
...
...     Returns whether they agree up to a phase, and whether it is the *same* phase
...     for every input state.
...     """
...     lie_trotter = transpile_trotter_step(
...         rows, cols, LieTrotter()
...     ).decompose(reps=6)
...     flow = transpile_trotter_step(
...         rows, cols, FlowSetSynthesis(rows, cols)
...     ).decompose(reps=6)
...     dimension = 2**lie_trotter.num_qubits
...     residuals, phases = [], []
...     for _ in range(samples):
...         vector = rng.normal(size=dimension) + 1j * rng.normal(size=dimension)
...         state = Statevector(vector / np.linalg.norm(vector))
...         a = state.evolve(lie_trotter).data
...         b = state.evolve(flow).data
...         overlap = np.vdot(a, b)
...         phases.append(overlap / abs(overlap))
...         residuals.append(np.linalg.norm(b - a * overlap / abs(overlap)))
...     # the SAME global phase must work for every input state
...     spread = max(abs(phase - phases[0]) for phase in phases)
...     return bool(max(residuals) < 1e-9), bool(spread < 1e-9)
>>>
>>> agrees_up_to_one_phase(rows, cols)
(True, True)

Nota

Entrambe le metà di quella verifica contano. Una sovrapposizione di modulo unitario su un singolo stato è necessaria ma non sufficiente: una sintesi che sbagliasse una fase dipendente dallo stato lo mostrerebbe comunque. Richiedere la stessa fase su diversi stati casuali è ciò che colma la lacuna.

9. Verificare la fisica sul sottospazio degli stabilizzatori

Questa è la complicazione della dimensione due. In 1D la codifica era un'isometria, e bastava confrontare gli spettri con un riferimento di Jordan-Wigner. VC non è un'isometria: usa $2N$ qubit per $N$ modi, quindi uno spazio di dimensione $2^{2N}$ ospita uno fisico di dimensione $2^N$. Questo è generico per le codifiche locali in due dimensioni (i qubit in eccesso sono vincolati da simmetrie locali commutanti, una per plaquette [2] [3]) e significa che l'Hamiltoniana codificata rappresenta fedelmente quella fermionica solo sull'autospazio congiunto $+1$ di quei vincoli.

Avviso

La parola stabilizzatore ha un doppio uso in questo ambito, e i due significati non hanno relazione tra loro. Ref. [1] la usa per gli operatori di trasferimento commutanti di un insieme di flusso, che formano un gruppo stabilizzatore che i Clifford del passo 5 diagonalizzano (un espediente di compilazione, e il senso inteso in tutto quanto sopra). Questa sezione la usa nel senso della correzione degli errori: operatori il cui autospazio $+1$ è il sottospazio fisico. Ref. [1] segnala esplicitamente questa collisione e mette da parte il secondo significato; i vincoli di plaquette qui sotto sono una proprietà della codifica VC [2], non un risultato del quadro degli insiemi di flusso.

I generatori sono costruiti qui anziché citati: uno per faccia quadrata, e ciascuno è semplicemente il ciclo di operatori di lato codificati attorno a quella faccia.

Gli operatori di lato sono l'oggetto naturale per questo. Un loop chiuso di essi è una simmetria di qualsiasi Hamiltoniana che conserva il numero di particelle, mentre il corrispondente loop di operatori di trasferimento non lo è: i due differiscono per un operatore di vertice sul sito sorgente di ciascun legame,

$$T_{jk} = \tfrac{i}{2} \, V_j E_{jk} \qquad \Longleftrightarrow \qquad E_{jk} = -2i \, V_j T_{jk} \, ,$$

Lungo un ciclo chiuso ogni sito è la sorgente di esattamente un legame, quindi il ciclo raccoglie un $V_j$ per sito, che è la parità locale dei quattro siti della plaquette. È quella parità locale a fare la differenza: l'Hamiltoniana fa saltare le particelle attraverso il bordo della plaquette, quindi conserva solo la parità globale. Il ciclo di trasferimento nudo quindi non commuta con essa, mentre il ciclo di lato (che porta il fattore di parità) commuta.

Poiché la codifica è definita sugli operatori di trasferimento, mappare un operatore di lato richiede la relazione precedente. È una callback di cinque righe, e si inserisce in map_edge_vertex_generators() esattamente come vc_action si inserisce in map_transfer_vertex_generators():

>>> from qiskit_fermions.mappers import map_edge_vertex_generators
>>> from qiskit_fermions.operators import EdgeVertexOperator
>>>
>>> def vc_edge_action(generator, width, cols):
...     """Map one generalized *edge* operator to its Verstraete-Cirac Pauli."""
...     j, k = generator
...     if j == k:  # a vertex operator is the same object in both algebras
...         return SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     # The library decomposes T_jk = (i/2) V_j E_jk, so E_jk = -2i V_j T_jk. Reuse the
...     # transfer encoding for T_jk rather than repeating its four cases here. Note that
...     # `compose` left-multiplies its argument, so this composes to V_j T_jk; the other
...     # order carries the opposite sign, since a vertex and a transfer operator sharing a
...     # mode anticommute.
...     vertex = SparseObservable.from_sparse_list([("Z", [j], 1.0)], num_qubits=width)
...     return (-2.0j * vc_action(generator, width, cols).compose(vertex)).simplify()
>>>
>>> def vc_edge_encoding(operator, width, cols):
...     """Map an EdgeVertexOperator to a qubit observable, Verstraete-Cirac style."""
...     return map_edge_vertex_generators(
...         operator,
...         partial(vc_edge_action, width=width, cols=cols),
...         identity=lambda: SparseObservable.identity(width),
...         compose=SparseObservable.compose,
...     ).simplify()
>>>
>>> def plaquette_loops(rows, cols):
...     """The directed edge loop around each square face."""
...     at = lambda r, c: r * cols + c
...     return [
...         [(at(r, c), at(r, c + 1)), (at(r, c + 1), at(r + 1, c + 1)),
...          (at(r + 1, c + 1), at(r + 1, c)), (at(r + 1, c), at(r, c))]
...         for r in range(rows - 1)
...         for c in range(cols - 1)
...     ]
>>>
>>> def stabilizers(rows, cols):
...     """One plaquette stabilizer per square face of the lattice."""
...     num_qubits = 2 * rows * cols
...     out = []
...     for loop in plaquette_loops(rows, cols):
...         loop_operator = EdgeVertexOperator.from_dict({tuple(loop): 1.0})
...         encoded_loop = vc_edge_encoding(loop_operator, num_qubits, cols)
...         out.append(encoded_loop / complex(encoded_loop.coeffs[0]))
...     return out
>>>
>>> plaquette_stabilizers = stabilizers(rows, cols)
>>> for stabilizer in plaquette_stabilizers:
...     print(stabilizer.to_sparse_list())
[('ZZYXXY', [1, 4, 9, 10, 12, 13], (1+0j))]
[('ZZYXXY', [2, 5, 10, 11, 13, 14], (1+0j))]
[('ZZYXXY', [4, 7, 12, 13, 15, 16], (1+0j))]
[('ZZYXXY', [5, 8, 13, 14, 16, 17], (1+0j))]

Ogni generatore ha peso $6$, distribuito su due qubit fisici della sua plaquette e quattro ancille. Sono involutori e commutano con l'Hamiltoniana codificata, quindi il sottospazio fisico è preservato dalla dinamica. Entrambe le proprietà sono affermazioni su operatori, quindi vengono verificate simbolicamente sull'algebra di SparseObservable, senza bisogno di matrici $2^{18} \times 2^{18}$:

>>> encoded = encoding(hamiltonian, num_qubits)
>>> identity = SparseObservable.identity(num_qubits)
>>> for stabilizer in plaquette_stabilizers:
...     square = stabilizer.compose(stabilizer).simplify()
...     print(square == identity, stabilizer.commutes(encoded))
True True
True True
True True
True True

Con gli stabilizzatori a disposizione, il confronto non richiede alcuna corrispondenza di base, perché $V_j = Z_j$ significa che le occupazioni vivono sugli stessi qubit fisici in entrambe le rappresentazioni. Quindi: si scrive uno stato prodotto direttamente sui qubit fisici con le ancille in $|0\rangle$, lo si proietta sul sottospazio degli stabilizzatori, lo si evolve con l'Hamiltoniana codificata e si leggono le densità dei siti $\langle n_j \rangle = (1 - \langle Z_j \rangle)/2$ da entrambi i lati.

Ogni passo qui rimane a livello di vettore. Il proiettore viene applicato un fattore $(1 + S)/2$ alla volta anziché essere moltiplicato per intero, e l'evoluzione temporale usa expm_multiply(), che calcola $e^{-iHt} |\psi\rangle$ senza mai formare $e^{-iHt}$. Poiché l'Hamiltoniana codificata è un operatore sparso di dimensione $2^{2N}$, questa distinzione è ciò che rende economica questa sezione.

>>> from scipy.sparse.linalg import expm_multiply
>>> from qiskit_fermions.mappers.library import transfer_vertex_jordan_wigner
>>>
>>> # an exact evolution time, unrelated to the Trotter step above
>>> evolution_time = 0.6
>>> reference = transfer_vertex_jordan_wigner(hamiltonian, num_sites).simplify()
>>>
>>> def to_sparse_matrix(observable):
...     # SparseObservable carries no dense- or sparse-matrix method, so route through
...     # SparsePauliOp. This is only ever needed for the numerical checks below --
...     # the encoding itself never leaves SparseObservable.
...     return SparsePauliOp.from_sparse_observable(observable).to_matrix(sparse=True)
>>>
>>> def basis_state(occupation, num_qubits):
...     """The computational basis state with `occupation` on the first qubits.
...
...     Qubit `q` carries bit `2**q`, and the ancillas are the *high* qubits `N..2N-1`,
...     so the same integer indexes the same site occupations on either register and
...     leaves every ancilla in |0>.
...     """
...     index = int("".join(str(bit) for bit in reversed(occupation)), 2)
...     state = np.zeros(2**num_qubits, dtype=complex)
...     state[index] = 1.0
...     return state
>>>
>>> def project_onto_codespace(state, stabilizers):
...     """Apply the projector one (1 + S) factor at a time, then renormalize.
...
...     The factors commute, so their product is itself a projector and the order does
...     not matter. Multiplying it out would build a dense 2^2N x 2^2N operator; applied
...     factor by factor this stays a handful of sparse matrix-vector products.
...
...     Each factor is really (1 + S)/2, but every one contributes the same scalar, so the
...     halvings only rescale the result and the final normalization absorbs them.
...     """
...     for stabilizer in stabilizers:
...         state += to_sparse_matrix(stabilizer) @ state
...     return state / np.linalg.norm(state)
>>>
>>> occupation = [1, 0, 1, 0, 1, 0, 0, 1, 0]  # site occupations of the initial state
>>> initial_reference = basis_state(occupation, num_sites)
>>> initial_encoded = project_onto_codespace(
...     basis_state(occupation, num_qubits), plaquette_stabilizers
... )
>>>
>>> evolved_encoded = expm_multiply(
...     -1j * evolution_time * to_sparse_matrix(encoded), initial_encoded
... )
>>> evolved_reference = expm_multiply(
...     -1j * evolution_time * to_sparse_matrix(reference), initial_reference
... )
>>>
>>> def density(state, qubit, width):
...     """Site density <n_j> = (1 - <Z_j>) / 2 of a statevector."""
...     z = to_sparse_matrix(
...         SparseObservable.from_sparse_list([("Z", [qubit], 1.0)], num_qubits=width)
...     )
...     return (1 - (state.conj() @ (z @ state)).real) / 2
>>>
>>> encoded_densities = [density(evolved_encoded, j, num_qubits) for j in range(num_sites)]
>>> reference_densities = [density(evolved_reference, j, num_sites) for j in range(num_sites)]
>>>
>>> for j, density_value in enumerate(encoded_densities):
...     print(f"site {j}: {density_value:.10f}")
site 0: 0.6750993356
site 1: 0.3689791258
site 2: 0.6750993356
site 3: 0.3818130396
site 4: 0.4868585996
site 5: 0.3818130396
site 6: 0.2693544696
site 7: 0.4916285849
site 8: 0.2693544696
>>> np.allclose(encoded_densities, reference_densities, atol=1e-12)
True

Queste sono le densità codificate; coincidono con il riferimento di Jordan-Wigner fino alla precisione di macchina, su un reticolo la cui codifica è fedele solo dopo la proiezione.

10. Misurare la convergenza di Trotter

Il passo 8 ha stabilito che le due sintesi concordano tra loro, il che è un'affermazione di equivalenza e non di accuratezza: entrambe descrivono lo stesso singolo passo di Trotter, quindi entrambe portano lo stesso errore di Trotter. Misurare quell'errore richiede un riferimento esatto, e il passo 9 ne ha appena costruito uno.

Il numero di passi proviene da FermionicLieTrotter, la controparte fermionica della formula a livello di qubit usata come baseline in precedenza. Applicata all'intero operatore raggruppato prima della scomposizione, la sua ripetizione cicla tra i quattro insiemi di flusso e il gruppo diagonale anziché all'interno di uno qualsiasi di essi:

>>> initial_state = project_onto_codespace(
...     basis_state(occupation, num_qubits), plaquette_stabilizers
... )
>>> exact = expm_multiply(
...     -1j * total_time * to_sparse_matrix(encoded), initial_state
... )
>>>
>>> def trotter_fidelity(steps):
...     # `steps` goes to the fermionic formula inside transpile_trotter_step, so the
...     # repetition cycles between the flow sets rather than within any one of them
...     circuit = transpile_trotter_step(
...         rows, cols, FlowSetSynthesis(rows, cols), steps=steps
...     )
...     evolved = Statevector(initial_state).evolve(circuit.decompose(reps=6))
...     return abs(np.vdot(evolved.data, exact)) ** 2
>>>
>>> for steps in (1, 2, 4, 8):
...     print(f"{steps} step(s): fidelity {trotter_fidelity(steps):.6f}")
1 step(s): fidelity 0.053820
2 step(s): fidelity 0.704014
4 step(s): fidelity 0.932937
8 step(s): fidelity 0.982687

La fedeltà sale monotonicamente verso uno al ridursi dei passi, che è il comportamento che una formula di prodotto deve mostrare. Accorciare i passi è solo una delle due manopole; FermionicSuzukiTrotter le bilancia l'una contro l'altra a un ordine superiore, come misura Sintetizzare un Gate Evolution nello spazio fermionico.

Un singolo passo di questa Hamiltoniana a $t = 1$ è un'approssimazione scarsa (i quattro insiemi di hopping e l'interazione non commutano), quindi la prima riga è lontana da uno, ed è proprio questo il punto: è la convergenza a confermare la costruzione, non una singola riga.

Nota

Ogni fattore della schedule copre ogni modo di questo reticolo, quindi restringere un fattore al suo supporto è un'operazione nulla qui e gli indici assoluti dei qubit su cui FlowSetSynthesis fa il confronto restano invariati. Una sintesi scritta per un operatore i cui gruppi toccano solo parte del registro dovrebbe tenere conto di quella rietichettatura.

Cosa portare a casa

• Due dimensioni cambiano la contabilità, non la ricetta. Ogni ingrediente della guida 1D si riporta invariato: groups per il raggruppamento, map_transfer_vertex_generators() per la codifica, CustomF2QLayout per il registro non corrispondente, e MapperFnEvolutionSynthesis per collegarlo al Transpiler. Ciò che è cresciuto è il numero di casi all'interno di ciascuno.

• La disgiunzione è ciò che rende costante la profondità in 2D. Un insieme di flusso si scompone in linee disgiunte sui vertici, quindi i suoi Clifford vengono emessi in parallelo e condividono due strati a due qubit indipendentemente dalle dimensioni del reticolo. Quattro insiemi di flusso e un gruppo diagonale danno quindi un budget fisso: qui $16$ per i quattro insiemi di flusso più $8$ per il gruppo di interazione diagonale.

• Ogni orientamento è una propria forma di Pauli. $T_{kj} = -V_j V_k T_{jk}$ trasforma XXZ in YYZ e XYXY in YXXY, quindi i quattro insiemi di flusso richiedono quattro Clifford. Supporre che un solo Clifford copra entrambi gli orientamenti di un asse è il modo più facile per sbagliare una sintesi di insiemi di flusso 2D.

• Leggi i supporti dalla coniugazione. In 1D una struttura a mattoni di CZ fornisce le immagini in forma chiusa. In 2D non esiste una formula simile, e le linee di colonna hanno ancille non contigue, quindi una congettura ricavata a mano su dove finisce ogni termine è fragile. Coniugare e ispezionare il risultato costa una costruzione di Clifford e non può essere sbagliato; significa anche che un Clifford ugualmente valido che porta i termini altrove continua a funzionare.

• Lascia che sia il raggruppamento a ordinare. Poiché Evolution si divide da sé per gruppo, alla sintesi viene consegnato un insieme di flusso alla volta e non deve mai ricavare di nuovo la partizione: legge una forma di Pauli dal primo termine e sceglie il Clifford corrispondente. Assumere questo, e sollevare un errore quando non vale, è ciò che mantiene breve synthesize.

• Il numero di passi sta sopra la divisione. La sintesi qui vede un insieme di flusso per Gate, quindi non può suddividere l'evoluzione tra gli insiemi. FermionicLieTrotter lo fa nello spazio fermionico, prima della scomposizione, ed è ciò che rende possibile un test di convergenza.

• Verifica sul sottospazio. Una codifica di $N$ modi con $2N$ qubit è fedele solo dove lo dicono gli stabilizzatori. Poiché VC mantiene $V_j = Z_j$, il controllo affidabile più economico è proiettare uno stato prodotto e confrontare le densità: niente intertwiner e niente confronto spettrale, che è inaffidabile quando l'Hamiltoniana ristretta è degenere.

Nota

L'Hamiltoniana codificata commuta con gli stabilizzatori, quindi la proiezione deve essere applicata solo allo stato iniziale. L'evoluzione preserva il sottospazio da sola. Su hardware reale questo è anche ciò che rende utile la codifica come risorsa di rilevamento degli errori: gli stabilizzatori possono essere misurati.

Prossimi passi

Due direzioni estendono quanto sopra, ed entrambe sono cambi di ingrediente e non di metodo. La prima è la codifica: qui è stata scelta Verstraete-Cirac perché Ref. [1] ne fornisce i Clifford diagonalizzanti in forma chiusa, ma Derby-Klassen impiega un'ancilla per plaquette dispari anziché una per sito, per un rapporto qubit-modi di $3/2$ anziché $2$. I suoi operatori di trasferimento hanno forme di Pauli diverse, quindi un porting equivale a riscrivere la funzione di codifica del passo 3 e a trovare il Clifford per ogni nuova forma; il raggruppamento, il layout e l'impalcatura della sintesi si riportano intatti.

La seconda è il reticolo stesso. Tutto qui presuppone condizioni al contorno aperte e una sola specie fermionica. Le condizioni al contorno periodiche chiudono ogni linea in un ciclo: le frecce commutano ancora tutte, dato che ciascuna incontra la successiva testa-coda, ma le due forme di Clifford del passo 5 sono costruite per una linea con estremi. Un ciclo di $L$ siti ha $L$ legami per orientamento dove una linea ne ha $L - 1$, quindi il legame di chiusura non trova posto nella costruzione e richiede una gestione a sé. Un modello con spin è più semplice: raddoppia il registro e aggiunge termini di interazione on-site al gruppo diagonale, senza toccare gli insiemi di flusso. Nessuno dei due altera l'argomento della profondità costante, ma entrambi cambiano la contabilità, che, come dice il primo punto, è la parte che due dimensioni ingrandiscono.

Importante

La connettività dell'hardware non è trattata qui. Ogni circuito precedente è costruito su un registro astratto in cui qualsiasi coppia di qubit può interagire, e le profondità riportate sono contate in quel contesto. I dispositivi reali hanno una mappa di accoppiamento ristretta, e mappare questa codifica su uno di essi è un problema sostanziale a sé stante: il registro VC è un grafo di $2N$ qubit i cui Gate agiscono tra un sito e la propria ancilla, tra siti vicini e tra ancille vicine, quindi un layout deve mantenere adiacenti contemporaneamente tutti e tre i tipi di coppie. Se non ci riesce, il Transpiler inserisce Gate SWAP, e gli strati che questi aggiungono sono ciò che la costruzione a profondità costante stava recuperando. Anche Ref. [1] cita la sua cifra di profondità sotto l'ipotesi di un grafo di connettività quadrato con CX come unico entangler nativo. Trovare buoni layout e strategie di routing per codifiche locali su una data topologia di dispositivo esula dallo scopo di questa guida.

[1] (1,2,3,4,5,6,7,8,9,10)

A. Gandon, S. Piccinelli, M. Rossmannek, F. Tacchino, A. Baiardi, J. Nys, and I. Tavernelli, Stabilizer-based quantum simulation of fermion dynamics with local qubit encodings (2026), arXiv:2512.11418v2.

[2] (1,2,3,4)

F. Verstraete and J. I. Cirac, Mapping local Hamiltonians of fermions to local Hamiltonians of spins, J. Stat. Mech. 2005, P09012 (2005).

[3] (1,2)

Y.-A. Chen and Y. Xu, Equivalence between fermion-to-qubit mappings in two spatial dimensions, PRX Quantum 4, 010326 (2023).