2D-Fermi-Hubbard-Dynamik mit Flow-Mengen simulieren

Wichtig

Die Konzepte in dieser Anleitung sind derzeit nur in der Python-API verfügbar. Entsprechende Funktionalität wird in einer zukünftigen Version in der C-API bereitgestellt.

Siehe auch

Lies zuerst die 1D-Flow-Mengen-Anleitung. Diese Anleitung ist ihre zweidimensionale Fortsetzung und wiederholt nicht, was dort entwickelt wurde: Transfer- und Vertex-Operatoren, was eine Flow-Menge ist, wie groups im Circuit erhalten bleibt und wie eine benutzerdefinierte Kodierung über MapperFnEvolutionSynthesis in den Transpiler eingebunden wird.

Siehe die Transpilationsanleitung für die Transpiler-Stufen, auf die durchgehend verwiesen wird, und die Mapper-Anleitung für das allgemeine Rezept zum Schreiben eines benutzerdefinierten Mappers.

Die 1D-Anleitung endete mit dem Hinweis, dass bei zwei Dimensionen lokale Kodierungen nötig sind. Diese Anleitung geht diesen Schritt, und zwar für das zweidimensionale Fermi-Hubbard-Modell auf einem quadratischen Gitter, unter Verwendung der Verstraete-Cirac-Kodierung (VC) [2].

Es gibt drei Unterschiede zum 1D-Problem, die die Grundlage dieser Arbeit bilden:

• Die Flow-Mengen vervielfachen sich. Eine Kette hat zwei Pfeilorientierungen; ein quadratisches Gitter hat vier (Ost, West, Nord, Süd), und jede ist eine Vereinigung knotendisjunkter gerichteter Pfade (einer pro Gitterzeile oder -spalte) statt eines einzelnen Pfades.

• Die Kodierung erhält ein Ancilla pro Platz, nicht eines für das ganze System, was ein

  Verhältnis von Qubits zu Moden von $2$ ergibt. Transferoperatoren werden zu Paulis mit Gewicht 3 (horizontal) und Gewicht 4 (vertikal), und jede der vier Flow-Mengen wird auf ihre eigene Pauli-Form abgebildet, sodass jede ihr eigenes diagonalisierendes Clifford braucht.

• Der kodierte Hamiltonoperator ist nur auf einem Unterraum treu: dem gemeinsamen $+1$-

  Eigenraum der Plaketten-Bedingungen, die jede lokale 2D-Kodierung mit sich bringt [2] [3]. In 1D war die Kodierung eine reine Isometrie und die Verifikation eine Spektrumsprüfung. Hier muss der Vergleich projiziert werden, wie Schritt 9 zeigt.

Der Gewinn hat dieselbe Form wie in 1D: Die Zwei-Qubit-Tiefe eines Trotter-Schritts wird konstant in der Gittergröße, nämlich $24$ ($16$ für das Hopping, passend zu Ref. [1], plus $8$ für die Wechselwirkung), während eine termweise Trotterisierung desselben kodierten Hamiltonoperators linear wächst. Die Gate-Anzahl wächst natürlich weiterhin mit der Anzahl der Bindungen; größenunabhängig wird, wie viele Schichten diese Gates belegen. Tiefen werden durchgehend auf einem abstrakten, vollständig verbundenen Register gezählt; die Kodierung auf die eingeschränkte Coupling Map eines Geräts abzubilden ist ein eigenes Problem, das am Ende besprochen wird.

Das Gitter und seine vier Orientierungen

Die erste Änderung ist am gerichteten Wechselwirkungsgraphen sichtbar. Die Plätze sind zeilenweise nummeriert, site(r, c) = r * cols + c, und jede Bindung trägt beide Orientierungen, wie die Kette in 1D, aber nun verlaufen die Bindungen entlang zweier Achsen:

Die Kommutationsregel ist die aus 1D, unverändert: Zwei Transferoperatoren, die einen Platz teilen, kommutieren, wenn die Pfeile durch ihn hindurchfließen, und antikommutieren, wenn sie kollidieren. Ein Platz hat jedoch nun bis zu vier anliegende Bindungen statt zwei, sodass die Pfeile durch ihn auf mehr Arten kollidieren können, und eine Menge von Pfeilen, die nie kollidiert, ist nicht mehr ein einzelner Pfad. Die horizontalen Pfeile von links nach rechts zu lesen ergibt eine solche Menge, die vertikalen Pfeile von oben nach unten eine weitere, und die Umkehrung jeder ergibt zwei weitere. Das sind die Ost-, West-, Nord- und Süd-Flow-Mengen, die alles Folgende ordnen.

Die Verstraete-Cirac-Kodierung

Verstraete-Cirac (VC) koppelt jeden fermionischen Mode $j$ an sein eigenes Hilfs-Qubit. Bei $N$ Gitterplätzen lebt Mode $j$ auf Qubit $j$ und sein Ancilla auf Qubit $N + j$, sodass das Register $2N$ Qubits enthält. Vertex-Operatoren bleiben so einfach wie möglich,

$$V_j = Z_j \, ,$$

was die Besetzungen später direkt ablesbar macht. Transferoperatoren hängen von der Kantenorientierung ab. Für eine horizontale Kante $j \to k$, wobei $j$ der westliche Platz ist, und eine vertikale Kante, wobei $j$ der nördliche Platz ist,

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

ein „Dreieck“ mit Gewicht 3 und ein „Quadrat“ mit Gewicht 4. Jeden Pfeil zu identifizieren ist wichtig: Mit $j$ als westlichem Platz zeigt $j \to k$ nach Osten, mit $j$ als nördlichem Platz jedoch nach Süden. Die umgekehrten Orientierungen folgen aus der Identität $T_{kj} = -V_j V_k T_{jk}$, die das XXZ-Dreieck zu YYZ und das XYXY-Quadrat zu YXXY macht. Diese vier verschiedenen Labels sind die vier Flow-Mengen, und der Rest dieser Anleitung verfolgt sie.

Hinweis

Dies sind nicht ganz Gl. (A1) aus Ref. [1], die „bis auf Einzel-Qubit-Rotationen“ formuliert ist und sich in drei Punkten unterscheidet, die man benennen sollte, da sie jedem auffallen, der das Paper nachprüft:

• Der Vorfaktor: Gl. (A1) enthält keinen. Der Faktor $\tfrac12$ hier sorgt dafür, dass $T_{jk}^2 = \tfrac14$ gilt, passend zur Normierung $T_{jk} = -\tfrac12 (a^\dagger_j - a_j)(a^\dagger_k + a_k)$, die TransferVertexOperator verwendet. Das Paper braucht seine Transferoperatoren immer nur innerhalb von $\exp(-i\,dt\,H)$, wo ein Skalierungsfaktor in den Rotationswinkel aufgenommen wird; hier müssen sie Operatoren sein, die eine Identität erfüllen.

• Die Orientierung: Gl. (A1) wählt die Buchstaben nach der Zeilenparität von $j$, die bei einer horizontalen Kante für beide Endpunkte gleich ist; deshalb unterscheidet sie $T_{jk}$ nicht von $T_{kj}$. Das Paper stellt die Orientierung separat durch Konjugation mit einer Lage $R^Z(\pi/2)$ wieder her. Dieser Leitfaden bindet sie stattdessen in die Kodierung ein und leitet die Gegenrichtung aus $T_{kj} = -V_j V_k T_{jk}$ ab, wodurch die beiden Orientierungen in Schritt 3 auf unterschiedlichen Paulis landen.

• Die vertikale Form: Gl. (A1) wechselt je nach Zeilenparität zwischen XXXX und YYYY, während dieser Leitfaden ein einzelnes XYXY verwendet. Das sind zwei Wege zum selben Ziel: eine Buchstabenasymmetrie, durch die eine vertikale Kante korrekt mit den horizontalen Kanten kollidiert, die auf sie treffen. Der Unterschied ist keine freie Basiswahl. Wer den Paritätswechsel weglässt und ein einzelnes XXXX verwendet, bricht die Algebra völlig, weshalb die unterschiedlichen Buchstaben geprüft statt angenommen werden.

Hinweis

Der Ancilla ist einem Site zugeordnet, was das Verhältnis $2$ ergibt; die Derby-Klassen-Kodierung dagegen platziert ein Hilfs-Qubit pro (ungerader) Plakette, für das kleinere Verhältnis $3/2$. VC wird hier verwendet, weil seine Transferoperatoren diejenigen sind, deren diagonalisierende Cliffords Ref. [1] in geschlossener Form angibt.

1. Der Hamiltonoperator auf einem Quadratgitter

Die Konstruktion ist dieselbe wie in 1D, nur mit einer zweidimensionalen Kantenliste. Die Sites behalten die zeilenweise Nummerierung der obigen Abbildung, site(r, c) = r * cols + c, sodass horizontale Bindungen aufeinanderfolgende Indizes verbinden und vertikale Bindungen Indizes, die sich um cols unterscheiden:

>>> 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. Vier Flow-Sets statt zwei

Jeder Hopping-Term wird danach klassifiziert, in welche Richtung sein Pfeil auf dem Gitter zeigt. Der Vergleich der divmod-Zeilen-/Spaltenkoordinaten der beiden Endpunkte trennt horizontale von vertikalen Kanten, und das Vorzeichen der Differenz bestimmt die Orientierung. Alles Diagonale (die Identität, die $V_j$- und die $V_j V_k$-Terme) kommt in eine fünfte Gruppe:

Hinweis

Orientierungskonvention. Die Namen sind einfach aus der obigen Abbildung abgelesen: Ein Pfeil ist east oder west, je nachdem, ob er entlang einer Zeile nach rechts oder links zeigt, und north oder south, je nachdem, ob er in einer Spalte nach oben oder unten zeigt. Da die Nummerierung zeilenweise läuft, sind east und south die Richtungen mit steigendem Index. Nichts im weiteren Verlauf hängt von den Namen ab (sie bezeichnen nur Gruppen), aber sie sind leicht zu verwechseln, daher lohnt es sich, sie anhand der Abbildung festzulegen. Ref. [1] benennt seine Flow-Sets auf dem Quadratgitter genauso.

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

Innerhalb eines Flow-Sets trifft jeder Pfeil entweder Kopf-an-Schwanz auf einen anderen oder gar nicht auf einen, wie in 1D; ein Flow-Set kommutiert also weiterhin Term für Term, und seine Evolution trägt keinen Trotter-Fehler. Neu ist, dass ein Flow-Set nun unzusammenhängend ist: Das east-Set eines rows x cols-Gitters besteht aus rows unabhängigen Pfaden von links nach rechts. Diese Unzusammenhängendheit ist ein Geschenk und keine Komplikation, wie Schritt 6 zeigt.

3. Den Operator kodieren

Die Kodierung ist eine einzelne Funktion, die einen einzelnen verallgemeinerten Transferoperator auf einen Pauli-String abbildet, der an map_transfer_vertex_generators() übergeben wird. Im Vergleich zur 1D-Version hat der Rumpf mehr Fälle: diagonal, horizontal, vertikal und die umgekehrten Orientierungen über $T_{kj} = -V_j V_k T_{jk}$.

Es gibt ein Problem. map_transfer_vertex_generators() ruft nur mit einem Generator zurück, während die Entscheidung über den Pauli-String einer Bindung auch cols benötigt, um Zeile und Spalte aus einem Modusindex zu gewinnen. Statt die Funktionen zu verschachteln, um darüber zu schließen, bindet functools.partial() den Wert, sodass beide Funktionen auf oberster Ebene bleiben und unabhängig aufrufbar sind.

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

Der Nutzen der Gruppierung vor der Abbildung ist sofort erkennbar: Jedes Flow-Set wird auf eine einzige Pauli-Form kodiert, einheitlich über das gesamte Set.

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

Die vier Hopping-Formen sind verschieden, daher erhält jedes Flow-Set sein eigenes Clifford in Schritt 5. Die umgekehrten Orientierungen sind andere Paulis, nicht dieselben umbenannt; deshalb braucht eine 2D-Flow-Set-Synthese vier Circuits, wo die 1D-Synthese einen braucht.

4. Transpilieren mit der benutzerdefinierten Kodierung

Diese Stufe ist die 1D-Stufe mit einem breiteren Qubit-Register: CustomF2QLayout bildet das $N$-Moden-Register auf $2N$ Qubits ab, MapperFnEvolutionSynthesis trägt die Kodierung, und sein Argument product_formula wählt, wie das resultierende PauliEvolutionGate zerlegt wird. Vergleiche Schritt 6 des 1D-Leitfadens, der dieselbe Pipeline für einen einzelnen Ancilla aufbaut.

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

Fünf Gates, die die in Schritt 2 gezählten 6, 6, 6, 6 und 22 Terme tragen. Jedes wird unabhängig abgebildet und synthetisiert, sodass die auf fermionischer Ebene getroffene Partition bis in den Circuit erhalten bleibt:

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

Das ist die Baseline: Kodierung und Gruppierung sind beide vorhanden, aber LieTrotter synthetisiert jedes Flow-Set Term für Term, und die Tiefe wächst mit dem Gitter.

5. Ein Clifford der Tiefe 2 pro Flow-Set

Der Mechanismus ist derselbe wie in 1D: Ein ganzes kommutierendes Flow-Set wird mit einem Clifford konjugiert, das jeden seiner Terme auf ein eigenes Pauli vom Gewicht 1 abbildet, sodass die Evolution des gesamten Sets nur Einzel-Qubit-Rotationen kostet. In 1D leistete das ein Brickwork aus CZ-Gates.

In 2D braucht jedes Flow-Set seinen eigenen Circuit, nach Abschnitt IV C und Abb. 7 von Ref. [1]. Alle vier wirken jeweils auf eine Linie des Gitters (eine Zeile für east/west, eine Spalte für north/south) zusammen mit den Ancillas dieser Linie, und alle vier haben dieselbe dreiteilige Form:

• Eine Lage von Einzel-Qubit-Rotationen, die jedes Qubit in die Basis bringt, die sein Flow-Set braucht.

• Eine erste CX-Lage, stets von Site zu Ancilla auf demselben Site.

• Eine zweite CX-Lage, die Nachbarn entlang der Linie koppelt.

Die Schritte 2 und 3 sind die einzigen Zwei-Qubit-Lagen, daher beträgt die Zwei-Qubit-Tiefe $2$, unabhängig davon, wie lang die Linie ist. Sie sind unten pro Flow-Set als eine Funktion ausgeschrieben. Die Rotationen sind so, wie sie sind, weil sie gefunden wurden, indem die Terme jedes Flow-Sets konjugiert wurden und verschiedene Bilder vom Gewicht 1 verlangt wurden; betrachte sie also als gegeben und verlasse dich auf die anschließende Prüfung.

Die beiden horizontalen Sets sind die einfachsten. Ihre Terme haben Gewicht 3, XXZ und YYZ, und unterscheiden sich nur durch das zusätzliche $S$, das auf den Sites $X$ in $Y$ verwandelt:

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

Die beiden vertikalen Sets haben Gewicht 4, XYXY und YXXY, und sind entsprechend aufwendiger: Die Rotationen wechseln mit der Position entlang der Linie, und auch die zweite CX-Lage wechselt, indem sie Sites auf geraden Positionen und Ancillas auf ungeraden koppelt (oder umgekehrt):

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

Auf einer einzelnen Linie mit drei Sites gezeichnet (zuerst ihre Sites, dann ihre Ancillas), passt der gesamte Mechanismus in ein Bild. Die Dreiecksform ist die klarere der beiden: ein Basiswechsel auf den Ancillas, dann zwei CX-Lagen, wobei die erste jeden Site mit seinem eigenen Ancilla paart und die zweite zum nächsten Site entlang der Linie reicht.

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

Die Quadratform braucht mehr Einzel-Qubit-Arbeit (die Rotationen wechseln entlang der Linie), und ihre zweite CX-Lage wechselt zwischen der Kopplung von Sites und der Kopplung von Ancillas:

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

Beide sind zwei CX-Lagen tief, und das Hinzufügen von Sites zur Linie macht sie breiter, aber nicht tiefer; das ist die gesamte Behauptung konstanter Tiefe, sichtbar auf einen Blick. Die west- und north-Cliffords unterscheiden sich nur in ihren Einzel-Qubit-Lagen und sind daher nicht separat gezeichnet.

Was die Synthese noch braucht, ist die Bestimmung, wo jeder Term landet. Eine Linie mit $L$ Sites trägt $L - 1$ Bindungen, die je einen Term beitragen; daher liefert der folgende Helfer einen Eintrag pro Bindung: das Qubit, auf dem der Term gelandet ist, und das Einzel-Qubit-Pauli, zu dem er wurde. Entscheidend ist, dass diese von der Konjugation abgelesen und nicht vorhergesagt werden. Anders als in 1D, wo ein CZ-Brickwork $X_{j+1} \mapsto Z_j X_{j+1} Z_{j+2}$ in geschlossener Form abbildet, gibt es hier keine solche Formel, und die Ancillas einer Spaltenlinie sind nicht zusammenhängend. Das Einzige, was sich zwischen den Flow-Sets unterscheidet, ist der eingehende Pauli-String, den Schritt 3 bereits ausgegeben hat:

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

Die Bilder müssen verschieden sein und Gewicht 1 haben, damit jeder Term unabhängig auf seinem eigenen Qubit rotiert werden kann; die Zwei-Qubit-Tiefe von $2$ ist in den vier obigen Konstruktionen bereits sichtbar, da jede zwei CX-Lagen ausgibt, unabhängig von der Linienlänge. Unten werden beide Eigenschaften über jedes Flow-Set, jede Linie und mehrere Gittergrößen geprüft. Schritt 6 ruft dann clifford_supports auf, es ist also ein funktionierender Teil der Synthese und nicht nur eine Demonstration.

6. Flow-Set-bewusste Synthese

Die Synthese bekommt einen hilfreichen Vorsprung. Da sich Evolution gruppenweise selbst aufteilt, wird synthesize einmal pro Flow-Set aufgerufen (hier fünfmal, vier Hopping-Sets plus die diagonale Gruppe), sodass jeder Hopping-Term, den es sieht, dieselbe Pauli-Form hat und ein Clifford alle diagonalisiert. Es muss diese Form nur am ersten Term ablesen und die übrigen danach sortieren, auf welcher Linie sie liegen. Dann folgt die Beobachtung, die zwei Dimensionen billiger macht, als man befürchten könnte:

Wichtig

Die Komponenten eines einzelnen Flow-Sets sind knotendisjunkt: Die Zeilen des west-Sets teilen keine Qubits. Daher können alle Cliffords eines Flow-Sets gemeinsam ausgegeben werden und teilen sich dieselben zwei Zwei-Qubit-Lagen, statt sich Zeile für Zeile hintereinander einzureihen. Ein Flow-Set kostet deshalb eine feste Anzahl von Lagen (zwei für das Clifford und zwei für sein Inverses), egal wie viele Zeilen oder Spalten das Gitter hat; deshalb bleibt die Gesamttiefe konstant.

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

Wie in 1D ist dies bewusst für diese Kodierung geschrieben: Es erkennt die vier Pauli-Formen und löst bei allem anderen eine Ausnahme aus, statt stillschweigend auszuweichen. Es setzt außerdem die oben beschriebene Aufteilung ein Flow-Set pro Gate voraus und sagt das auch; mehr als eine Hopping-Form in einem einzelnen Gate bedeutet, dass die Gruppierung nicht angewendet wurde, was ein Fehler ist und nichts, was man umgehen sollte.

Hinweis

Wie sein 1D-Gegenstück nimmt diese Synthese keine Anzahl von Trotter-Schritten entgegen. Ein Trotter-Schritt muss zwischen den Flow-Sets wechseln, aber transpile_trotter_step unten zerlegt den fermionischen Circuit in ein Evolution-Gate pro Set, sodass diese Methode immer nur ein Set auf einmal sieht und nur innerhalb davon wiederholen könnte, was nichts ändert: Die Wiederholung eines Operators teilt lediglich seine eigene Evolution in gleiche Stücke.

Mehrschritt-Trotterisierung braucht daher eine Produktformel, die die groups des Operators innerhalb eines einzelnen Gates berücksichtigt, vor der Zerlegung. FermionicLieTrotter leistet das, und Schritt 10 nutzt es, um die Konvergenz zu messen.

Beachte die beiden ähnlich benannten Klassen im Spiel. Qiskits LieTrotter, oben als Baseline verwendet, wirkt auf den abgebildeten Hamiltonoperator, in dem die Gruppierung bereits verloren ist. FermionicLieTrotter ist das fermionische Gegenstück, und es ist die Klasse, die die Flow-Sets noch sehen kann.

7. Das Ergebnis: konstante Tiefe

Alles bisher lief auf dem einzelnen $3 \times 3$-Gitter aus Schritt 1. Der Vergleich der beiden Synthesepfade über wachsende Gitter erfordert, dass die gesamte Abfolge (Operator, Gruppierung, Kodierung, Circuit, Pass Manager) für jede Größe wiederholt wird; daher lohnt es sich, sie einmal zu kapseln:

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

Beide Metriken über einen Bereich von Gittergrößen zu sammeln und gegen die Gitterkante aufzutragen, macht den Unterschied deutlich:

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

Die Zwei-Qubit-Tiefe ist 24 und konstant von $3 \times 3$ bis $8 \times 8$ (das sind 18 bis 128 Qubits), während LieTrotter, angewandt auf den identisch kodierten Hamiltonoperator, linear mit der Gitterkante wächst, von 47 auf 92. Die Gate-Anzahl wächst bei beiden, da ohnehin $O(N)$ Terme anzuwenden sind; was die Flow-Sets beseitigen, ist das Wachstum der Zahl der Lagen, die diese Gates belegen.

Das Budget teilt sich auf in 16 für das Hopping und acht für die Wechselwirkung:

• Jedes der vier Flow-Sets kostet $2 + 2$ Zwei-Qubit-Lagen: sein Clifford beim Hineingehen und das Inverse beim Herauskommen, beide konstruktionsbedingt von Tiefe zwei, mit nur Einzel-Qubit-Rotationen dazwischen. Das ergibt $4 \times 4 = 16$, die Zahl, die Ref. [1] nennt, da sie nur die Hopping-Terme zählt. Setzt man im obigen Code $U = 0$, erhält man sie reproduziert.

• Die diagonale Gruppe trägt die restlichen acht bei. Ihre Rzz-Gates packen sich in vier Lagen, und jedes Rzz zerfällt in zwei CX-Gates um eine Einzel-Qubit-Rotation, also $4 \times 2 = 8$. Vier Lagen sind optimal: Auf einem Quadratgitter reicht die Wechselwirkung nur bis zu nächsten Nachbarn, sodass jeder innere Site vier Rzz-Gates trägt (eines pro Bindung) und immer nur an einem gleichzeitig teilnehmen kann.

Nichts in einer der beiden Summen bezieht sich auf die Gittergröße, und genau das ist der Punkt.

8. Den Circuit verifizieren

Die Flow-Set-Synthese und LieTrotter beschreiben denselben Trotter-Schritt desselben kodierten Hamiltonoperators, daher müssen ihre Unitäre bis auf eine globale Phase übereinstimmen. Zufällige Statevektoren erkennen jede Abweichung weit günstiger, als die volle Unitäre aufzubauen:

Hinweis

Beide Circuits werden vor der Anwendung auf einen Zustand in Ein- und Zwei-Qubit-Gates zerlegt. Ohne Zerlegung zwingt jede High-Level-Evolution-Anweisung Statevector, ihren vollen dichten $2^{2N} \times 2^{2N}$-Operator aufzubauen; wird zuerst zerlegt, kann es stattdessen die Gates einzeln anwenden, was hier mehrere hundert Mal schneller ist.

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

Hinweis

Beide Hälften dieser Prüfung sind wichtig. Ein Überlapp mit Betrag eins auf einem einzelnen Zustand ist notwendig, aber nicht hinreichend: Eine Synthese, die eine zustandsabhängige Phase falsch hat, würde ihn trotzdem zeigen. Dieselbe Phase über mehrere Zufallszustände zu verlangen, schließt diese Lücke.

9. Die Physik auf dem Stabilisator-Unterraum verifizieren

Dies ist die zweidimensionale Besonderheit. In 1D war die Kodierung eine Isometrie, und der Vergleich der Spektren mit einer Jordan-Wigner-Referenz genügte. VC ist keine Isometrie: Es verwendet $2N$ Qubits für $N$ Moden, sodass ein $2^{2N}$-dimensionaler Raum einen $2^N$-dimensionalen physikalischen enthält. Das ist typisch für lokale Kodierungen in zwei Dimensionen (die überzähligen Qubits werden durch lokale kommutierende Symmetrien eingeschränkt, eine pro Plakette [2] [3]), und es bedeutet, dass der kodierte Hamiltonoperator den fermionischen nur auf dem gemeinsamen $+1$-Eigenraum dieser Nebenbedingungen getreu darstellt.

Warnung

Das Wort Stabilisator hat in diesem Gebiet eine Doppelrolle, und die beiden Bedeutungen sind unabhängig voneinander. Ref. [1] verwendet es für die kommutierenden Transferoperatoren eines Flow-Sets, die eine Stabilisatorgruppe bilden, die die Cliffords aus Schritt 5 diagonalisieren (ein Kompilierungskniff und die Bedeutung, die oben überall gemeint ist). Dieser Abschnitt verwendet es im Sinne der Fehlerkorrektur: Operatoren, deren $+1$-Eigenraum der physikalische Unterraum ist. Ref. [1] weist ausdrücklich auf diese Kollision hin und legt die zweite Bedeutung beiseite; die Plakettenbedingungen unten sind eine Eigenschaft der VC-Kodierung [2] und kein Ergebnis des Flow-Set-Rahmens.

Die Generatoren werden hier konstruiert statt zitiert: einer pro Quadratfläche, und jeder ist einfach die kodierte Kanten-Operator-Schleife um diese Fläche.

Kantenoperatoren sind dafür das natürliche Objekt. Eine geschlossene Schleife aus ihnen ist eine Symmetrie jedes teilchenzahlerhaltenden Hamiltonoperators, wohingegen die entsprechende Schleife aus Transferoperatoren es nicht ist: Die beiden unterscheiden sich durch einen Vertex-Operator auf dem Quell-Site jeder Bindung,

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

Um eine geschlossene Schleife ist jeder Site Quelle von genau einer Bindung, sodass die Schleife ein $V_j$ pro Site einsammelt, nämlich die lokale Parität der vier Sites der Plakette. Diese lokale Parität macht den Unterschied: Der Hamiltonoperator lässt Teilchen über den Rand der Plakette hüpfen und erhält daher nur die globale Parität. Die reine Transferschleife kommutiert deshalb nicht mit ihr, während die Kantenschleife (die den Paritätsfaktor trägt) kommutiert.

Da die Kodierung auf Transferoperatoren definiert ist, braucht die Abbildung eines Kantenoperators die obige Beziehung. Das ist ein Callback von fünf Zeilen, und er passt in map_edge_vertex_generators() genauso, wie vc_action in map_transfer_vertex_generators() passt:

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

Jeder Generator hat Gewicht $6$, verteilt auf zwei physikalische Qubits seiner Plakette und vier Ancillas. Sie sind involutorisch und kommutieren mit dem kodierten Hamiltonoperator, sodass der physikalische Unterraum von der Dynamik erhalten bleibt. Beide Eigenschaften sind Aussagen über Operatoren, daher werden sie symbolisch auf der SparseObservable-Algebra geprüft, ganz ohne $2^{18} \times 2^{18}$-Matrizen:

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

Mit den Stabilisatoren in der Hand braucht der Vergleich keine Basiszuordnung, denn $V_j = Z_j$ bedeutet, dass die Besetzungen in beiden Bildern auf denselben physikalischen Qubits liegen. Also: einen Produktzustand direkt auf den physikalischen Qubits schreiben, mit den Ancillas in $|0\rangle$, ihn auf den Stabilisator-Unterraum projizieren, mit dem kodierten Hamiltonoperator evolvieren und die Site-Dichten $\langle n_j \rangle = (1 - \langle Z_j \rangle)/2$ auf beiden Seiten ablesen.

Jeder Schritt hier bleibt auf der Ebene eines Vektors. Der Projektor wird Faktor für Faktor als $(1 + S)/2$ angewandt, statt ausmultipliziert zu werden, und die Zeitentwicklung verwendet expm_multiply(), das $e^{-iHt} |\psi\rangle$ berechnet, ohne $e^{-iHt}$ je zu bilden. Da der kodierte Hamiltonoperator ein dünnbesetzter $2^{2N}$-dimensionaler Operator ist, hält diese Unterscheidung diesen Abschnitt günstig.

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

Das sind die kodierten Dichten; sie stimmen mit der Jordan-Wigner-Referenz bis auf Maschinengenauigkeit überein, auf einem Gitter, dessen Kodierung erst nach der Projektion getreu ist.

10. Die Trotter-Konvergenz messen

Schritt 8 hat gezeigt, dass die beiden Synthesen miteinander übereinstimmen, was eine Aussage über Äquivalenz und nicht über Genauigkeit ist: Beide beschreiben denselben einzelnen Trotter-Schritt, also tragen beide denselben Trotter-Fehler. Diesen Fehler zu messen, braucht eine exakte Referenz, und Schritt 9 hat gerade eine gebaut.

Die Schrittzahl stammt von FermionicLieTrotter, dem fermionischen Gegenstück der zuvor als Baseline verwendeten Formel auf Qubit-Ebene. Auf den gesamten gruppierten Operator vor der Zerlegung angewandt, wechselt seine Wiederholung zwischen den vier Flow-Sets und der diagonalen Gruppe, statt innerhalb einer von ihnen zu bleiben:

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

Die Fidelity steigt monoton gegen eins, je kürzer die Schritte werden, was das Verhalten ist, das eine Produktformel zeigen muss. Die Verkürzung der Schritte ist nur einer von zwei Reglern; FermionicSuzukiTrotter wägt sie bei höherer Ordnung gegeneinander ab, wie Einen Evolution-Gate im fermionischen Raum synthetisieren misst.

Ein einzelner Schritt dieses Hamiltonoperators bei $t = 1$ ist eine schlechte Näherung (die vier Hopping-Sets und die Wechselwirkung kommutieren nicht), daher ist die erste Zeile weit von eins entfernt, und genau das ist der Punkt: Die Konvergenz bestätigt die Konstruktion, nicht irgendeine einzelne Zeile.

Hinweis

Jeder Faktor des Schedules umfasst jede Mode dieses Gitters, daher ist das Einengen eines Faktors auf seinen Träger hier ein No-op, und die absoluten Qubit-Indizes, auf die FlowSetSynthesis passt, bleiben unverändert. Eine Synthese, die für einen Operator geschrieben ist, dessen Gruppen nur einen Teil des Registers berühren, müsste diese Umbenennung berücksichtigen.

Was man mitnehmen sollte

• Zwei Dimensionen ändern die Buchhaltung, nicht das Rezept. Jede Zutat aus dem 1D-Leitfaden bleibt unverändert erhalten: groups für die Gruppierung, map_transfer_vertex_generators() für die Kodierung, CustomF2QLayout für das nicht passende Register und MapperFnEvolutionSynthesis, um sie in den Transpiler einzubinden. Gewachsen ist die Zahl der Fälle in jeder dieser Zutaten.

• Disjunktheit macht 2D zu konstanter Tiefe. Ein Flow-Set zerfällt in knotendisjunkte Linien, daher werden seine Cliffords parallel ausgegeben und teilen sich zwei Zwei-Qubit-Lagen, egal wie groß das Gitter ist. Vier Flow-Sets und eine diagonale Gruppe ergeben dann ein festes Budget: hier $16$ für die vier Flow-Sets plus $8$ für die diagonale Wechselwirkungsgruppe.

• Jede Orientierung ist eine eigene Pauli-Form. $T_{kj} = -V_j V_k T_{jk}$ verwandelt XXZ in YYZ und XYXY in YXXY, daher brauchen die vier Flow-Sets vier Cliffords. Anzunehmen, dass ein Clifford beide Orientierungen einer Achse abdeckt, ist der einfachste Weg, eine 2D-Flow-Set-Synthese falsch zu machen.

• Träger von der Konjugation ablesen. In 1D liefert ein CZ-Brickwork die Bilder in geschlossener Form. In 2D gibt es keine solche Formel, und Spaltenlinien haben nicht zusammenhängende Ancillas, sodass eine von Hand abgeleitete Vermutung, wo jeder Term landet, fragil ist. Das Konjugieren und Untersuchen des Ergebnisses kostet eine Clifford-Konstruktion und kann nicht falsch sein; es bedeutet auch, dass ein ebenso gültiges Clifford, das die Terme woanders landen lässt, weiter funktioniert.

• Lass die Gruppierung das Sortieren übernehmen. Da sich Evolution pro Gruppe selbst aufteilt, bekommt die Synthese jeweils ein Flow-Set und muss die Partition nie neu herleiten: Sie liest eine Pauli-Form am ersten Term ab und wählt das passende Clifford. Das anzunehmen und eine Ausnahme auszulösen, wenn es nicht zutrifft, hält synthesize kurz.

• Die Schrittzahl liegt oberhalb der Aufteilung. Die Synthese hier sieht ein Flow-Set pro Gate und kann die Evolution daher nicht über die Sets hinweg unterteilen. FermionicLieTrotter tut das im fermionischen Raum, vor der Zerlegung, und das macht einen Konvergenztest überhaupt erst möglich.

• Auf dem Unterraum verifizieren. Eine $2N$-Qubit-Kodierung von $N$ Moden ist nur dort getreu, wo die Stabilisatoren es sagen. Da VC $V_j = Z_j$ beibehält, ist die günstigste zuverlässige Prüfung, einen Produktzustand zu projizieren und Dichten zu vergleichen: kein Intertwiner und kein Spektralabgleich, der unzuverlässig ist, wenn der eingeschränkte Hamiltonoperator entartet ist.

Hinweis

Der kodierte Hamiltonoperator kommutiert mit den Stabilisatoren, daher muss die Projektion nur auf den Anfangszustand angewandt werden. Die Evolution erhält den Unterraum von selbst. Auf echter Hardware macht das die Kodierung auch als Ressource zur Fehlererkennung nützlich: Die Stabilisatoren können gemessen werden.

Als Nächstes

Zwei Richtungen erweitern das Obige, und beide sind Änderungen einer Zutat und nicht der Methode. Die erste ist die Kodierung: Verstraete-Cirac wurde hier gewählt, weil Ref. [1] seine diagonalisierenden Cliffords in geschlossener Form angibt, aber Derby-Klassen verwendet einen Ancilla pro ungerader Plakette statt einen pro Site, für ein Qubit-zu-Moden-Verhältnis von $3/2$ statt $2$. Seine Transferoperatoren haben andere Pauli-Formen, sodass eine Portierung darauf hinausläuft, die Kodierungsfunktion aus Schritt 3 umzuschreiben und für jede neue Form das Clifford zu finden; die Gruppierung, das Layout und das Synthesegerüst bleiben unberührt.

Die zweite ist das Gitter selbst. Alles hier setzt offene Ränder und eine einzelne fermionische Spezies voraus. Periodische Ränder schließen jede Linie zu einem Zyklus: Die Pfeile kommutieren weiterhin alle, da jeder den nächsten Kopf-an-Schwanz trifft, aber die beiden Clifford-Formen aus Schritt 5 sind für eine Linie mit Endpunkten gebaut. Ein Zyklus aus $L$ Sites trägt $L$ Bindungen pro Orientierung, wo eine Linie $L - 1$ trägt, sodass die Wrap-around-Bindung in der Konstruktion keinen Platz hat und eine eigene Behandlung braucht. Ein spinbehaftetes Modell ist einfacher: Es verdoppelt das Register und fügt der diagonalen Gruppe On-Site-Wechselwirkungsterme hinzu, ohne die Flow-Sets zu berühren. Keines von beiden stört das Argument der konstanten Tiefe, aber beide ändern die Buchhaltung, die, wie die erste Erkenntnis sagt, der Teil ist, den zwei Dimensionen größer machen.

Wichtig

Die Hardware-Konnektivität wird hier nicht behandelt. Jeder obige Circuit ist auf einem abstrakten Register aufgebaut, in dem jedes Qubit-Paar wechselwirken kann, und die angegebenen Tiefen werden in diesem Rahmen gezählt. Echte Geräte haben eine eingeschränkte Coupling Map, und die Abbildung dieser Kodierung darauf ist ein eigenes, umfangreiches Problem: Das VC-Register ist ein $2N$-Qubit-Graph, dessen Gates zwischen einem Site und seinem eigenen Ancilla, zwischen benachbarten Sites und zwischen benachbarten Ancillas verlaufen, sodass ein Layout alle drei Arten von Paaren gleichzeitig benachbart halten muss. Gelingt das nicht, fügt der Transpiler SWAP-Gates ein, und die Lagen, die sie hinzufügen, sind genau das, was die Konstruktion konstanter Tiefe zurückgekauft hat. Ref. [1] nennt ebenso seine Tiefenzahl unter einem angenommenen quadratischen Konnektivitätsgraphen mit CX als einzigem nativen Verschränker. Gute Layouts und Routing-Strategien für lokale Kodierungen auf einer gegebenen Gerätetopologie zu finden, liegt jenseits des Umfangs dieses Leitfadens.

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