Simulace 2D Fermi-Hubbardovy dynamiky pomocí množin toků

Důležité

Koncepty v tomto průvodci jsou aktuálně dostupné pouze v Python API. Ekvivalentní funkcionalita bude v C API zpřístupněna v budoucím vydání.

Viz také

Nejprve si přečti průvodce 1D množinami toků. Tento průvodce je jeho dvourozměrným pokračováním a neopakuje informace, které tam byly rozvinuty: transferové a vertexové operátory, co je množina toků, jak groups přežije do obvodu a jak se vlastní kódování zapojuje do transpilátoru přes MapperFnEvolutionSynthesis.

Fáze transpilátoru, na které se odkazujeme v celém textu, najdeš v průvodci transpilací a obecný recept na psaní vlastního mapperu v průvodci mappery.

1D průvodce skončil zmínkou, že při práci ve dvou dimenzích jsou potřeba lokální kódování. Tento průvodce tento krok provádí, pro dvourozměrný Fermi-Hubbardův model na čtvercové mřížce, s použitím Verstraete-Ciracova (VC) kódování [2].

Oproti 1D problému jsou tři rozdíly, které tvoří základ této práce:

• Množiny toků se množí. Řetězec má dvě orientace šipek; čtvercová mřížka má čtyři (východ, západ, sever, jih) a každá je sjednocením vrcholově disjunktních orientovaných cest (jedna na řádek či sloupec mřížky), a ne jedinou cestou.

• Kódování získává ancillu na každé místo, nikoli jednu pro celý systém, což dává

  poměr qubitů k módům $2$. Transferové operátory se stávají Pauliho operátory o váze 3 (horizontální) a váze 4 (vertikální) a každá ze čtyř množin toků se zobrazí na vlastní Pauliho tvar, takže každá potřebuje vlastní diagonalizující Cliffordův obvod.

• Zakódovaný hamiltonián je věrný jen na podprostoru: společném $+1$

  vlastním podprostoru plakettových omezení, která nese každé 2D lokální kódování [2] [3]. V 1D bylo kódování prostou izometrií a ověření bylo kontrolou spektra. Zde se srovnání musí promítnout, jak ukazuje krok 9.

Výsledek má stejný tvar jako v 1D: dvouqubitová hloubka jednoho Trotterova kroku se stává konstantní vzhledem k velikosti mřížky, a to na $24$ ($16$ pro přeskoky, shodně s Ref. [1], plus $8$ pro interakci), zatímco trotterizace téhož zakódovaného hamiltoniánu člen po členu roste lineárně. Počet hradel samozřejmě stále roste s počtem vazeb; na velikosti nezávislým se stává to, kolik vrstev tato hradla zabírají. Hloubky se všude počítají na abstraktním, plně propojeném registru; zobrazení kódování na omezenou mapu propojení zařízení je samostatný problém, o kterém se pojednává na konci.

Mřížka a její čtyři orientace

První změna je viditelná na orientovaném grafu interakcí. Místa jsou číslována po řádcích, site(r, c) = r * cols + c, a každá vazba nese obě orientace, jako řetězec v 1D, ale nyní vazby vedou podél dvou os:

Komutační pravidlo je 1D pravidlo beze změny: dva transferové operátory sdílející místo komutují, když šipky tímto místem protékají, a antikomutují, když se střetnou. Místo však má nyní až čtyři přilehlé vazby místo dvou, takže se šipky v něm mohou střetnout více způsoby a množina šipek, které se nikdy nestřetnou, už není jedinou cestou. Čtení horizontálních šipek zleva doprava dává jednu takovou množinu, vertikální šipky shora dolů další a obrácení každé dává další dvě. To jsou východní, západní, severní a jižní množiny toků, které organizují vše následující.

Verstraete-Ciracovo kódování

Verstraete-Cirac (VC) páruje každý fermionový mód $j$ s vlastním pomocným qubitem. Při $N$ místech mřížky mód $j$ žije na qubitu $j$ a jeho ancilla je na qubitu $N + j$, takže registr obsahuje $2N$ qubitů. Vertexové operátory zůstávají tak jednoduché, jak jen mohou být,

$$V_j = Z_j \, ,$$

což umožňuje později přímo číst obsazení. Transferové operátory závisí na orientaci hrany. Pro horizontální hranu $j \to k$, kde $j$ je západní místo, a vertikální hranu, kde $j$ je severní místo,

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

váhy 3 „trojúhelník“ a váhy 4 „čtverec“. Identifikace každé šipky je důležitá: když je $j$ západní místo, $j \to k$ míří na východ, ale když je $j$ severní místo, míří na jih. Opačné orientace vyplývají z identity $T_{kj} = -V_j V_k T_{jk}$, která mění trojúhelník XXZ na YYZ a čtverec XYXY na YXXY. Těmito čtyřmi odlišnými označeními jsou čtyři množiny toků a zbytek tohoto průvodce je sleduje.

Poznámka

Toto nejsou přesně rovnice (A1) z Ref. [1], která uvádí „až na jednoqubitové rotace“ a liší se ve třech ohledech, které stojí za pojmenování, protože si jich čtenář porovnávající text s článkem všimne:

• Prefaktor: Rovnice (A1) žádný nemá. Zde $\tfrac12$ zajišťuje, že $T_{jk}^2 = \tfrac14$, což odpovídá normalizaci $T_{jk} = -\tfrac12 (a^\dagger_j - a_j)(a^\dagger_k + a_k)$, kterou používá TransferVertexOperator. Článek své transferové operátory potřebuje pouze uvnitř $\exp(-i\,dt\,H)$, kde se škálovací faktor pohltí do úhlu rotace; zde to ale musí být operátory, které splňují určitou identitu.

• Orientace: Rovnice (A1) vybírá písmena podle parity řádku $j$, která je u vodorovné hrany pro oba koncové body stejná, takže nerozlišuje $T_{jk}$ od $T_{kj}$. Článek orientaci obnovuje zvlášť konjugací vrstvou $R^Z(\pi/2)$. Tento návod ji místo toho zapracovává přímo do kódování a opačnou orientaci odvozuje z $T_{kj} = -V_j V_k T_{jk}$, díky čemuž obě orientace v kroku 3 skončí na různých Paulih.

• Svislý tvar: Rovnice (A1) střídá XXXX a YYYY podle parity řádku, zatímco tento návod používá jediné XYXY. Jsou to dva způsoby, jak dosáhnout téhož: asymetrie písmen, díky níž svislá hrana správně koliduje s vodorovnými hranami, které se s ní setkávají. Nepovažuj ten rozdíl za volnou volbu báze. Vynechání střídání podle parity a použití jediného XXXX algebru přímo rozbije, proto se odlišná písmena kontrolují, a ne předpokládají.

Poznámka

Ancilla je připojena k uzlu, což dává poměr $2$; kódování Derby-Klassen naopak umisťuje pomocný qubit na každou (lichou) plaketu, což dává menší poměr $3/2$. VC se zde používá proto, že jeho transferové operátory jsou ty, jejichž diagonalizující Cliffordovy operace uvádí Ref. [1] v uzavřeném tvaru.

1. Hamiltonián na čtvercové mřížce

Konstrukce je stejná jako v 1D, jen s dvourozměrným seznamem hran. Uzly si zachovávají číslování po řádcích použité u obrázku výše, site(r, c) = r * cols + c, takže vodorovné vazby spojují po sobě jdoucí indexy a svislé vazby spojují indexy lišící se o 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. Čtyři toková rozdělení místo dvou

Každý hoppingový člen se klasifikuje podle směru, kterým jeho šipka na mřížce míří. Porovnání souřadnic řádku a sloupce (divmod) obou koncových bodů odliší vodorovné hrany od svislých a znaménko rozdílu určí orientaci. Vše diagonální (identita, členy $V_j$ a $V_j V_k$) patří do páté skupiny:

Poznámka

Konvence orientace. Názvy jsou jednoduše odečtené z obrázku výše: šipka míří na východ nebo západ, pokud ukazuje doprava nebo doleva podél řádku, a na sever nebo jih, pokud ukazuje nahoru nebo dolů podél sloupce. Protože číslování probíhá po řádcích, východ a jih jsou směry, ve kterých index roste. Nic dalšího na názvech nezávisí (jen označují skupiny), ale snadno se pletou, proto stojí za to je ověřit podle obrázku. Ref. [1] pojmenovává toková rozdělení na čtvercové mřížce stejně.

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

Uvnitř jednoho tokového rozdělení se každá šipka s jinou buď setkává hlavou na patu, nebo se s ní nesetkává vůbec, stejně jako v 1D, takže tokové rozdělení stále komutuje člen po členu a jeho evoluce nenese žádnou Trotterovu chybu. Novinkou je, že tokové rozdělení je nyní nesouvislé: východní rozdělení mřížky rows x cols se skládá z rows nezávislých cest zleva doprava. Toto rozpojení je spíš dar než komplikace, jak ukazuje krok 6.

3. Zakódování operátoru

Kódování je jedna funkce, která převádí jediný zobecněný transferový operátor na Pauliho řetězec, jenž se předá funkci map_transfer_vertex_generators(). Oproti verzi pro 1D má tělo více případů: diagonální, vodorovný, svislý a opačné orientace přes $T_{kj} = -V_j V_k T_{jk}$.

Je tu jeden problém. Funkce map_transfer_vertex_generators() volá zpět s generátorem a ničím jiným, zatímco rozhodnutí o Pauliho řetězci vazby potřebuje také cols, aby z indexu módu zjistilo souřadnice řádku a sloupce. Místo vnořování funkcí, které by ho uzavíraly, ho svazuje functools.partial(), takže obě funkce zůstávají na nejvyšší úrovni a lze je volat nezávisle.

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

Smysl seskupování před mapováním je vidět okamžitě: každé tokové rozdělení se zakóduje do jediného Pauliho tvaru, jednotného v celé množině.

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

Čtyři hoppingové tvary jsou navzájem různé, takže každé tokové rozdělení dostane vlastní Cliffordovu operaci v kroku 5. Opačné orientace jsou jiné Pauliho operátory, ne tytéž s přejmenovanými indexy; proto syntéza 2D tokových rozdělení potřebuje čtyři obvody, zatímco syntéza v 1D jen jeden.

4. Transpilace s vlastním kódováním

Tato fáze je stejná jako v 1D, jen s širším registrem qubitů: CustomF2QLayout mapuje registr $N$ módů na $2N$ qubitů, MapperFnEvolutionSynthesis nese kódování a jeho argument product_formula určuje, jak se výsledné PauliEvolutionGate dekomponuje. Porovnej s krokem 6 návodu pro 1D, který staví stejnou pipeline pro jedinou ancillu.

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

Pět hradel, nesoucích 6, 6, 6, 6 a 22 členů spočítaných v kroku 2. Každé se mapuje a syntetizuje nezávisle, takže rozdělení určené na fermionové úrovni přežije až do obvodu:

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

To je základní úroveň: kódování i seskupení jsou na místě, ale LieTrotter syntetizuje každé tokové rozdělení člen po členu a hloubka roste s velikostí mřížky.

5. Cliffordova operace hloubky 2 pro každé tokové rozdělení

Mechanismus je stejný jako v 1D: celé komutující tokové rozdělení se konjuguje Cliffordovou operací, která každý jeho člen zobrazí na odlišný Pauliho operátor váhy 1, takže evoluce celé množiny stojí jen jednoqubitové rotace. V 1D to zajišťovala jediná cihlová vrstva hradel CZ.

Ve 2D potřebuje každé tokové rozdělení vlastní obvod podle oddílu IV C a obr. 7 v Ref. [1]. Všechny čtyři působí vždy na jednu linii mřížky (řádek pro východ/západ, sloupec pro sever/jih) spolu s ancillami této linie a všechny mají stejný třídílný tvar:

• Vrstva jednoqubitových rotací, která uvede každý qubit do báze, kterou jeho tokové rozdělení potřebuje.

• První vrstva CX, vždy z uzlu na ancillu téhož uzlu.

• Druhá vrstva CX, která propojuje sousedy podél linie.

Kroky 2 a 3 jsou jediné dvouqubitové vrstvy, takže dvouqubitová hloubka je $2$ bez ohledu na délku linie. Níže jsou rozepsány po jedné funkci na tokové rozdělení. Rotace jsou takové, jaké jsou, protože byly nalezeny konjugací členů každého tokového rozdělení s požadavkem na odlišné obrazy váhy 1, takže je ber jako dané a spoléhej na následující kontrolu.

Obě vodorovná rozdělení jsou nejjednodušší. Jejich členy mají váhu 3, XXZ a YYZ, a liší se jen dodatečným $S$, které na uzlech mění $X$ na $Y$:

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

Obě svislá rozdělení mají váhu 4, XYXY a YXXY, a jsou odpovídajícím způsobem složitější: rotace se střídají podle pozice podél linie a druhá vrstva CX se také střídá, spojuje uzly na sudých pozicích a ancilly na lichých (nebo naopak):

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

Nakreslený na jediné třiuzlové linii (nejdřív její uzly, pak jejich ancilly) se celý mechanismus vejde do jednoho obrázku. Tvar trojúhelníku je z těch dvou přehlednější: změna báze na ancillách, pak dvě vrstvy CX, první páruje každý uzel s jeho vlastní ancillou a druhá sahá k dalšímu uzlu podél linie.

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

Tvar čtverce potřebuje více jednoqubitové práce (rotace se střídají podél linie) a jeho druhá vrstva CX se střídá mezi propojováním uzlů a propojováním ancill:

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

Obě mají hloubku dvou vrstev CX a přidání uzlů do linie je rozšiřuje, ale neprohlubuje; to je celé tvrzení o konstantní hloubce, viditelné pouhým pohledem. Západní a severní Cliffordovy operace se liší jen jednoqubitovými vrstvami, proto nejsou kresleny zvlášť.

Syntéza stále potřebuje zjistit, kam který člen skončí. Linie o $L$ uzlech nese $L - 1$ vazeb, z nichž každá přispívá jedním členem, takže pomocná funkce níže vrací jeden záznam na vazbu: qubit, na který člen dopadl, a jednoqubitový Pauliho operátor, kterým se stal. Podstatné je, že se tyto hodnoty odečítají z konjugace, a nepředpovídají. Na rozdíl od 1D, kde jedna cihlová vrstva CZ zobrazuje $X_{j+1} \mapsto Z_j X_{j+1} Z_{j+2}$ v uzavřeném tvaru, zde žádný takový vzorec neexistuje a ancilly sloupcové linie nejsou souvislé. Jediné, co se mezi tokovými rozděleními liší, je vstupní Pauliho řetězec, který už vypsal krok 3:

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

Obrazy musí vyjít odlišné a váhy 1, aby každý člen šlo rotovat nezávisle na jeho vlastním qubitu; dvouqubitová hloubka $2$ je již vidět ve čtyřech konstrukcích výše, z nichž každá vysílá dvě vrstvy CX bez ohledu na délku linie. Níže se obě vlastnosti ověřují pro každé tokové rozdělení, každou linii a několik velikostí mřížky. Krok 6 pak volá clifford_supports, takže je to funkční součást syntézy, a ne jen ukázka.

6. Syntéza respektující tokové rozdělení

Syntéza má užitečný náskok. Protože se Evolution dělí sama po skupinách, synthesize se volá jednou na tokové rozdělení (zde pětkrát, čtyři hoppingová rozdělení plus diagonální skupina), takže každý hoppingový člen, který vidí, sdílí jediný Pauliho tvar a jedna Cliffordova operace diagonalizuje všechny. Stačí přečíst tvar z prvního členu a zbylé seřadit podle linie, na které leží. Pak přichází pozorování, díky němuž jsou dvě dimenze levnější, než by se dalo obávat:

Důležité

Komponenty jednoho tokového rozdělení jsou vrcholově disjunktní: řádky západního rozdělení nesdílejí žádné qubity. Všechny Cliffordovy operace tokového rozdělení proto lze vyslat společně a sdílejí stejné dvě dvouqubitové vrstvy, místo aby čekaly v řadě jeden řádek po druhém. Tokové rozdělení tedy stojí pevný počet vrstev (dvě pro Cliffordovu operaci a dvě pro její inverzi) bez ohledu na to, kolik řádků nebo sloupců mřížka má, a proto je celková hloubka konstantní.

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

Stejně jako v 1D je to záměrně napsáno pro toto kódování: rozpoznává čtyři Pauliho tvary a u čehokoli jiného vyvolá chybu, místo aby tiše přešlo na záložní řešení. Také předpokládá výše uvedené rozdělení jedno tokové rozdělení na hradlo a říká to; více než jeden hoppingový tvar v jediném hradle znamená, že se seskupení neuplatnilo, což je chyba, ne něco, co se má obcházet.

Poznámka

Stejně jako její protějšek v 1D, ani tato syntéza nebere počet Trotterových kroků. Trotterův krok musí cyklit mezi tokovými rozděleními, ale transpile_trotter_step níže dekomponuje fermionový obvod na jedno hradlo Evolution na každé rozdělení, takže tato metoda vidí vždy jen jedno rozdělení a mohla by se opakovat jen uvnitř něj, což nic nemění: opakování jednoho operátoru jen rozdělí jeho vlastní evoluci na stejné díly.

Vícekroková Trotterizace proto potřebuje součinový vzorec, který respektuje groups operátoru uvnitř jednoho hradla, před dekompozicí. FermionicLieTrotter to dělá a krok 10 ho používá k měření konvergence.

Všimni si dvou podobně pojmenovaných tříd. LieTrotter z Qiskitu, použitý výše jako základní úroveň, působí na namapovaný hamiltonián, kde už je seskupení pryč. FermionicLieTrotter je jeho fermionový protějšek a je to ten, který tokové rozdělení ještě vidí.

7. Výsledek: konstantní hloubka

Vše dosud běželo na jediné mřížce $3 \times 3$ vytvořené v kroku 1. Porovnání obou syntetických cest na rostoucích mřížkách vyžaduje opakovat celou tu posloupnost (operátor, seskupení, kódování, obvod, pass manager) pro každou velikost, proto ji stojí za to jednou zabalit:

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

Sesbírání obou metrik pro řadu velikostí mřížky a jejich vynesení proti hraně mřížky ten rozdíl dobře ukáže:

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

Dvouqubitová hloubka je 24 a konstantní od $3 \times 3$ do $8 \times 8$ (tedy od 18 qubitů do 128), zatímco LieTrotter aplikovaný na stejně zakódovaný hamiltonián roste lineárně s hranou mřížky, od 47 do 92. Počet hradel roste u obou, protože je třeba aplikovat $O(N)$ členů tak či onak; tokové rozdělení odstraňuje růst počtu vrstev, které tato hradla zabírají.

Rozpočet se dělí na 16 pro hopping a osm pro interakci:

• Každé ze čtyř tokových rozdělení stojí $2 + 2$ dvouqubitové vrstvy: jeho Cliffordova operace na vstupu a inverze na výstupu, obě podle konstrukce hloubky dva, s pouze jednoqubitovými rotacemi mezi nimi. To je $4 \times 4 = 16$, což je údaj, který uvádí Ref. [1], protože počítá jen hoppingové členy. Nastavení $U = 0$ ve výše uvedeném kódu ho reprodukuje.

• Diagonální skupina přidává zbývajících osm. Její hradla Rzz se vejdou do čtyř vrstev a každé Rzz se dekomponuje na dvě hradla CX kolem jedné jednoqubitové rotace, takže $4 \times 2 = 8$. Čtyři vrstvy jsou optimum: na čtvercové mřížce interakce sahá jen k nejbližším sousedům, takže každý vnitřní uzel nese čtyři hradla Rzz (jedno na vazbu) a může se účastnit vždy jen jednoho.

Nic v žádném z těch součtů se nevztahuje k velikosti mřížky, což je celý smysl.

8. Ověření obvodu

Syntéza respektující tokové rozdělení a LieTrotter popisují tentýž Trotterův krok téhož zakódovaného hamiltoniánu, takže jejich unitární operátory se musí shodovat až na globální fázi. Náhodné stavové vektory odhalí jakoukoli odchylku mnohem levněji než sestavení celého unitárního operátoru:

Poznámka

Oba obvody se dekomponují na jedno- a dvouqubitová hradla před aplikací na stav. Pokud by zůstaly nedekomponované, každá vysokoúrovňová instrukce Evolution by nutila Statevector sestavit její plný hustý operátor $2^{2N} \times 2^{2N}$; dekompozice napřed mu umožní aplikovat hradla po jednom, což je zde několik set krát rychlejší.

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

Poznámka

Obě poloviny té kontroly jsou důležité. Překryv s jednotkovým modulem na jediném stavu je nutný, ale ne postačující: syntéza, která se u fáze závislé na stavu zmýlila, by to stále ukázala. Požadavek na stejnou fázi napříč několika náhodnými stavy tu mezeru uzavírá.

9. Ověření fyziky na podprostoru stabilizátorů

Toto je dvourozměrný zádrhel. V 1D bylo kódování izometrií a stačilo porovnat spektra s referencí Jordan-Wigner. VC izometrií není: používá $2N$ qubitů pro $N$ módů, takže $2^{2N}$-rozměrný prostor obsahuje $2^N$-rozměrný fyzikální. To je typické pro lokální kódování ve dvou dimenzích (přebytečné qubity jsou vázány lokálními komutujícími symetriemi, jednou na plaketu [2] [3]) a znamená to, že zakódovaný hamiltonián reprezentuje fermionový věrně jen na společném vlastním podprostoru $+1$ těchto vazeb.

Varování

Slovo stabilizátor má v tomto oboru dvojí význam a oba spolu nesouvisejí. Ref. [1] ho používá pro komutující transferové operátory tokového rozdělení, které tvoří stabilizátorovou grupu diagonalizovanou Cliffordovými operacemi z kroku 5 (kompilační pomůcka a význam použitý všude výše). Tato sekce ho používá ve smyslu kvantové korekce chyb: operátory, jejichž vlastní podprostor $+1$ je fyzikálním podprostorem. Ref. [1] tuto kolizi výslovně zmiňuje a druhý význam odkládá stranou; plaketové vazby níže jsou vlastností kódování VC [2], ne výsledkem rámce tokových rozdělení.

Generátory se zde konstruují, a necitují: jeden na čtvercovou stěnu a každý je jednoduše zakódovaná smyčka hranových operátorů kolem této stěny.

Hranové operátory jsou pro to přirozeným objektem. Uzavřená smyčka z nich je symetrií každého hamiltoniánu zachovávajícího počet částic, zatímco odpovídající smyčka transferových operátorů ne: liší se vrcholovým operátorem na zdrojovém uzlu každé vazby,

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

Kolem uzavřené smyčky je každý uzel zdrojem právě jedné vazby, takže smyčka sbírá jedno $V_j$ na uzel, což je lokální parita čtyř uzlů plakety. Právě ta lokální parita dělá rozdíl: hamiltonián přenáší částice přes hranici plakety, takže zachovává jen globální paritu. Holá transferová smyčka s ním proto nekomutuje, zatímco hranová smyčka (nesoucí paritní faktor) ano.

Protože je kódování definováno na transferových operátorech, mapování hranového operátoru potřebuje výše uvedený vztah. To je pětiřádkový callback a zapadá do map_edge_vertex_generators() přesně tak, jako vc_action zapadá do 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))]

Každý generátor má váhu $6$, rozloženou na dva fyzikální qubity své plakety a čtyři ancilly. Jsou involutorní a komutují se zakódovaným hamiltoniánem, takže dynamika zachovává fyzikální podprostor. Obě vlastnosti jsou tvrzení o operátorech, takže se kontrolují symbolicky na algebře SparseObservable, bez matic $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

Se stabilizátory v ruce porovnání nepotřebuje žádnou bázovou korespondenci, protože $V_j = Z_j$ znamená, že obsazení žijí na týchž fyzikálních qubitech v obou obrazech. Takže: zapiš součinový stav přímo na fyzikální qubity s ancillami v $|0\rangle$, promítni ho na podprostor stabilizátorů, vyvíjej se zakódovaným hamiltoniánem a z obou stran odečti hustoty uzlů $\langle n_j \rangle = (1 - \langle Z_j \rangle)/2$.

Každý krok zde zůstává na úrovni vektoru. Projektor se aplikuje po jednom faktoru $(1 + S)/2$, místo aby se roznásobil, a časový vývoj používá expm_multiply(), který počítá $e^{-iHt} |\psi\rangle$, aniž by kdy sestavil $e^{-iHt}$. Protože je zakódovaný hamiltonián řídký $2^{2N}$-rozměrný operátor, právě to činí tuto sekci levnou.

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

To jsou zakódované hustoty; shodují se s referencí Jordan-Wigner až na strojovou přesnost, na mřížce, jejíž kódování je věrné až po projekci.

10. Měření Trotterovy konvergence

Krok 8 ukázal, že se obě syntézy navzájem shodují, což je tvrzení o ekvivalenci, ne o přesnosti: obě popisují tentýž jediný Trotterův krok, takže obě nesou stejnou Trotterovu chybu. Měření té chyby vyžaduje přesnou referenci a krok 9 ji právě vytvořil.

Počet kroků pochází z FermionicLieTrotter, fermionového protějšku vzorce na úrovni qubitů použitého výše jako základní úroveň. Aplikovaný na celý seskupený operátor před dekompozicí cykluje své opakování mezi čtyřmi tokovými rozděleními a diagonální skupinou, a ne uvnitř kteréhokoli z nich:

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

Věrnost monotónně roste k jedné, jak se kroky zkracují, což je chování, které součinový vzorec musí vykazovat. Zkracování kroků je jen jedním ze dvou ovladačů; FermionicSuzukiTrotter je mezi sebou vyvažuje při vyšším řádu, jak měří Syntéza hradla Evolution ve fermionovém prostoru.

Jediný krok tohoto hamiltoniánu při $t = 1$ je slabá aproximace (čtyři hoppingová rozdělení a interakce nekomutují), takže první řádek je daleko od jedné, a o to jde: konstrukci potvrzuje konvergence, ne jakýkoli jednotlivý řádek.

Poznámka

Každý faktor rozvrhu zasahuje všechny módy této mřížky, takže zúžení faktoru na jeho nosič je zde bezvýznamné a absolutní indexy qubitů, podle kterých FlowSetSynthesis páruje, se nemění. Syntéza napsaná pro operátor, jehož skupiny se dotýkají jen části registru, by musela s tímto přeznačením počítat.

Co si odnést

• Dvě dimenze mění účetnictví, ne recept. Každá ingredience z návodu pro 1D se přenáší beze změny: groups pro seskupení, map_transfer_vertex_generators() pro kódování, CustomF2QLayout pro nesouladný registr a MapperFnEvolutionSynthesis pro zapojení do transpilátoru. Vyrostl počet případů uvnitř každého z nich.

• Disjunktnost dělá 2D konstantně hlubokým. Tokové rozdělení se rozpadá na vrcholově disjunktní linie, takže jeho Cliffordovy operace se vysílají paralelně a sdílejí dvě dvouqubitové vrstvy bez ohledu na velikost mřížky. Čtyři tokové rozdělení a diagonální skupina pak dávají pevný rozpočet: zde $16$ pro čtyři tokové rozdělení plus $8$ pro diagonální interakční skupinu.

• Každá orientace je vlastní Pauliho tvar. $T_{kj} = -V_j V_k T_{jk}$ mění XXZ na YYZ a XYXY na YXXY, takže čtyři tokové rozdělení potřebují čtyři Cliffordovy operace. Předpoklad, že jedna Cliffordova operace pokryje obě orientace osy, je nejsnazší způsob, jak syntézu 2D tokových rozdělení pokazit.

• Nosiče odečítej z konjugace. V 1D dává jedna cihlová vrstva CZ obrazy v uzavřeném tvaru. Ve 2D takový vzorec není a sloupcové linie mají nesouvislé ancilly, takže ručně odvozený odhad, kam který člen dopadne, je křehký. Konjugace a prohlídka výsledku stojí jednu konstrukci Clifford a nemůže se mýlit; znamená to také, že stejně platná Cliffordova operace, která členy umístí jinam, dál funguje.

• Nech seskupení třídit. Protože se Evolution dělí sama podle skupin, syntéza dostává vždy jedno tokové rozdělení a nikdy nemusí rozdělení znovu odvozovat: přečte jeden Pauliho tvar z prvního členu a vybere odpovídající Cliffordovu operaci. Tento předpoklad a vyvolání chyby, když neplatí, drží synthesize krátké.

• Počet kroků leží nad rozdělením. Syntéza zde vidí jedno tokové rozdělení na hradlo, takže nemůže evoluci rozdělit napříč množinami. FermionicLieTrotter to dělá ve fermionovém prostoru, před dekompozicí, což vůbec umožňuje test konvergence.

• Ověřuj na podprostoru. $2N$-qubitové kódování $N$ módů je věrné jen tam, kde to stabilizátory říkají. Protože VC zachovává $V_j = Z_j$, nejlevnější spolehlivou kontrolou je promítnout součinový stav a porovnat hustoty: žádný intertwiner a žádné párování spekter, které je nespolehlivé, když je omezený hamiltonián degenerovaný.

Poznámka

Zakódovaný hamiltonián komutuje se stabilizátory, takže projekci stačí aplikovat jen na počáteční stav. Evoluce zachovává podprostor sama. Na skutečném hardwaru právě to činí kódování užitečným jako zdroj detekce chyb: stabilizátory lze měřit.

Další

Dva směry rozšiřují to, co je výše, a oba jsou změnou ingredience, ne metody. Prvním je kódování: Verstraete-Cirac byl zvolen zde, protože Ref. [1] uvádí jeho diagonalizující Cliffordovy operace v uzavřeném tvaru, ale Derby-Klassen spotřebuje jednu ancillu na lichou plaketu místo jedné na uzel, tedy poměr qubitů k módům $3/2$ místo $2$. Jeho transferové operátory mají jiné Pauliho tvary, takže přenos spočívá v přepsání kódovací funkce z kroku 3 a nalezení Cliffordovy operace pro každý nový tvar; seskupení, layout a lešení syntézy se přenášejí beze změny.

Druhým je samotná mřížka. Vše zde předpokládá otevřené okrajové podmínky a jediný fermionový druh. Periodické okrajové podmínky uzavírají každou linii do cyklu: šipky stále všechny komutují, protože každá se s další setkává hlavou na patu, ale dva Cliffordovy tvary z kroku 5 jsou vytvořeny pro linii s koncovými body. Cyklus o $L$ uzlech nese $L$ vazeb na orientaci tam, kde linie nese $L - 1$, takže vazba přes okraj nemá v konstrukci místo a vyžaduje vlastní zacházení. Spinový model je snazší: zdvojnásobí registr a přidá do diagonální skupiny členy interakce na uzlu, aniž by se dotkl tokových rozdělení. Ani jedno neruší argument konstantní hloubky, ale obojí mění účetnictví, což je, jak říká první poznatek, část, kterou dvě dimenze zvětšují.

Důležité

Konektivita hardwaru se zde neřeší. Každý obvod výše je postaven na abstraktním registru, kde mohou interagovat libovolné dva qubity, a uváděné hloubky se počítají v tomto nastavení. Skutečná zařízení mají omezenou mapu propojení a mapování tohoto kódování na ni je samo o sobě podstatný problém: registr VC je graf o $2N$ qubitech, jehož hradla běží mezi uzlem a jeho vlastní ancillou, mezi sousedními uzly a mezi sousedními ancillami, takže layout musí udržet všechny tři druhy dvojic současně sousedící. Pokud to nejde, transpilátor vloží hradla SWAP a vrstvy, které přidají, jsou přesně to, co konstrukce s konstantní hloubkou ušetřila. Ref. [1] rovněž uvádí svůj údaj o hloubce za předpokladu čtvercového grafu konektivity s CX jako jediným nativním entanglerem. Hledání dobrých layoutů a směrovacích strategií pro lokální kódování na dané topologii zařízení přesahuje rozsah tohoto návodu.

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