Simule a dinâmica de Fermi-Hubbard 2D com conjuntos de fluxo

Importante

Os conceitos deste guia estão disponíveis atualmente apenas na API Python. Funcionalidade equivalente será disponibilizada na API C em uma versão futura.

Veja também

Leia primeiro o guia de conjuntos de fluxo 1D. Este guia é a sequência bidimensional dele e não repete as informações desenvolvidas lá: operadores de transferência e de vértice, o que é um conjunto de fluxo, como groups sobrevive no circuito e como uma codificação personalizada é integrada ao transpilador por meio de MapperFnEvolutionSynthesis.

Veja o guia de transpilação para as etapas do transpilador citadas ao longo do texto, e o guia de mapeadores para a receita geral de escrita de um mapeador personalizado.

O guia 1D terminou mencionando que codificações locais são necessárias ao trabalhar com duas dimensões. Este guia dá esse passo, para o modelo de Fermi-Hubbard bidimensional em uma rede quadrada, usando a codificação de Verstraete-Cirac (VC) [2].

Há três diferenças em relação ao problema 1D, que formam a base deste trabalho:

• Os conjuntos de fluxo se multiplicam. Uma cadeia tem duas orientações de seta; uma rede quadrada tem quatro (leste, oeste, norte, sul), e cada uma é uma união de caminhos direcionados disjuntos em vértices (um por linha ou coluna da rede), e não um único caminho.

• A codificação ganha um ancilla por sítio, e não um para o sistema inteiro, resultando em uma

  razão de $2$ Qubits por modo. Os operadores de transferência passam a ser Paulis de peso 3 (horizontais) e de peso 4 (verticais), e cada um dos quatro conjuntos de fluxo mapeia para seu próprio formato de Pauli, então cada um precisa de seu próprio Clifford diagonalizante.

• O Hamiltoniano codificado só é fiel em um subespaço: o $+1$ conjunto

  conjunto $+1$ das restrições de plaqueta que toda codificação local 2D carrega [2] [3]. Em 1D a codificação era uma isometria simples e a verificação era uma checagem de espectro. Aqui a comparação precisa ser projetada, como a etapa 9 ilustra.

O ganho tem o mesmo formato que em 1D: a profundidade de dois Qubits de um passo de Trotter se torna constante no tamanho da rede, em $24$ ($16$ para o salto, coincidindo com a Ref. [1], mais $8$ para a interação), enquanto uma Trotterização termo a termo do mesmo Hamiltoniano codificado cresce linearmente. A contagem de portas ainda cresce com o número de ligações, é claro; o que se torna independente do tamanho é quantas camadas essas portas ocupam. As profundidades em todo o texto são contadas em um registrador abstrato totalmente conectado; mapear a codificação no mapa de acoplamento restrito de um dispositivo é um problema à parte, discutido no final.

A rede e suas quatro orientações

A primeira mudança é visível no grafo de interação direcionado. Os sítios são numerados em ordem de linhas, site(r, c) = r * cols + c, e cada ligação carrega as duas orientações, como a cadeia em 1D, mas agora as ligações correm ao longo de dois eixos:

A regra de comutação é a do caso 1D, sem alteração: dois operadores de transferência que compartilham um sítio comutam quando as setas fluem através dele e anticomutam quando colidem. No entanto, um sítio agora tem até quatro ligações incidentes em vez de duas, então as setas que passam por ele podem colidir de mais maneiras, e um conjunto de setas que nunca colidem deixa de ser um único caminho. Ler as setas horizontais da esquerda para a direita dá um desses conjuntos, as setas verticais de cima para baixo dão outro, e invertê-las dá mais dois. Esses são os conjuntos de fluxo leste, oeste, norte e sul que organizam tudo o que vem a seguir.

A codificação de Verstraete-Cirac

Verstraete-Cirac (VC) emparelha cada modo fermiônico $j$ com seu próprio Qubit auxiliar. Com $N$ sítios na rede, o modo $j$ fica no Qubit $j$ e seu ancilla fica no Qubit $N + j$, então o registrador tem $2N$ Qubits. Os operadores de vértice permanecem tão simples quanto possível,

$$V_j = Z_j \, ,$$

o que torna as ocupações diretamente legíveis mais adiante. Os operadores de transferência dependem da orientação da aresta. Para uma aresta horizontal $j \to k$ com $j$ o sítio oeste, e uma aresta vertical com $j$ o sítio norte,

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

um “triângulo” de peso 3 e um “quadrado” de peso 4. Identificar cada seta importa: com $j$ como o sítio oeste, $j \to k$ aponta para leste, mas com $j$ como o sítio norte, aponta para o sul. As orientações inversas decorrem da identidade $T_{kj} = -V_j V_k T_{jk}$, que transforma o triângulo XXZ em YYZ, e o quadrado XYXY em YXXY. Esses quatro rótulos distintos são os quatro conjuntos de fluxo, e o restante deste guia os acompanha.

Nota

Estas não são exatamente a Eq. (A1) da Ref. [1], que afirma “até rotações de um único qubit” e difere em três pontos que vale nomear, já que quem conferir o artigo vai notar:

• O prefator: a Eq. (A1) não tem nenhum. O $\tfrac12$ aqui é o que faz $T_{jk}^2 = \tfrac14$, em linha com a normalização $T_{jk} = -\tfrac12 (a^\dagger_j - a_j)(a^\dagger_k + a_k)$ que TransferVertexOperator usa. O artigo só precisa dos seus operadores de transferência dentro de $\exp(-i\,dt\,H)$, onde um fator de escala é absorvido no ângulo de rotação; aqui eles precisam ser operadores que satisfazem uma identidade.

• A orientação: a Eq. (A1) escolhe as letras pela paridade da linha de $j$, que para uma aresta horizontal é a mesma nos dois extremos, então ela não distingue $T_{jk}$ de $T_{kj}$. O artigo recupera a orientação separadamente, conjugando com uma camada de $R^Z(\pi/2)$. Este guia, em vez disso, incorpora a orientação na codificação, derivando a inversa de $T_{kj} = -V_j V_k T_{jk}$, o que faz as duas orientações caírem em Paulis diferentes no passo 3.

• O formato vertical: a Eq. (A1) alterna XXXX e YYYY pela paridade da linha, enquanto este guia usa um único XYXY. São duas maneiras de obter a mesma coisa: uma assimetria de letras que faz uma aresta vertical colidir corretamente com as arestas horizontais que a encontram. Não leia a diferença como uma escolha livre de base. Abandonar a alternância de paridade e usar um único XXXX quebra a álgebra por completo, e é por isso que as letras distintas são verificadas em vez de presumidas.

Nota

O ancilla é associado a um sítio, dando a razão $2$; é a codificação de Derby-Klassen que coloca um qubit auxiliar por plaqueta (ímpar), para a razão menor $3/2$. O VC é usado aqui porque seus operadores de transferência são aqueles cujos Cliffords diagonalizadores a Ref. [1] apresenta em forma fechada.

1. O Hamiltoniano em uma rede quadrada

A construção é a mesma do caso 1D, mas com uma lista de arestas bidimensional. Os sítios mantêm a numeração por linhas usada na figura acima, site(r, c) = r * cols + c, de modo que as ligações horizontais conectam índices consecutivos, e as ligações verticais conectam índices que diferem por 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. Quatro conjuntos de fluxo em vez de dois

Cada termo de salto é classificado pela direção para a qual sua seta aponta na rede. Comparar as coordenadas de linha/coluna divmod dos dois extremos separa as arestas horizontais das verticais, e o sinal da diferença define a orientação. Tudo que é diagonal (a identidade, os termos $V_j$ e $V_j V_k$) vai para um quinto grupo:

Nota

Convenção de orientação. Os nomes são simplesmente lidos da figura acima: uma seta é leste ou oeste conforme aponta para a direita ou para a esquerda ao longo de uma linha, e norte ou sul conforme aponta para cima ou para baixo ao longo de uma coluna. Como a numeração segue as linhas, leste e sul são as direções de índice crescente. Nada adiante depende dos nomes (eles apenas rotulam grupos), mas é fácil confundi-los, então vale fixá-los em relação à imagem. A Ref. [1] nomeia seus conjuntos de fluxo da rede quadrada da mesma forma.

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

Dentro de um conjunto de fluxo, cada seta ou encontra outra cabeça com cauda ou não encontra nenhuma, como em 1D, então um conjunto de fluxo ainda comuta termo a termo e sua evolução não carrega erro de Trotter. O que é novo é que um conjunto de fluxo agora é desconexo: o conjunto leste de uma rede rows x cols consiste em rows caminhos independentes da esquerda para a direita. Essa desconexão é uma vantagem, e não uma complicação, como mostra o passo 6.

3. Codifique o operador

A codificação é uma função que mapeia um único operador de transferência generalizado para uma string de Pauli, que é passada para map_transfer_vertex_generators(). Em comparação com a versão 1D, o corpo tem mais casos: diagonal, horizontal, vertical e as orientações inversas por meio de $T_{kj} = -V_j V_k T_{jk}$.

Há um problema. map_transfer_vertex_generators() chama de volta com um gerador e nada mais, enquanto decidir a string de Pauli de uma ligação também exige cols para recuperar as coordenadas de linha e coluna a partir de um índice de modo. Em vez de aninhar as funções para capturá-lo, functools.partial() o vincula, o que mantém as duas funções no nível superior e chamáveis de forma independente.

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

O ponto de agrupar antes de mapear pode ser visto de imediato: cada conjunto de fluxo é codificado em um único formato de Pauli, uniforme em todo o conjunto.

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

Os quatro formatos de salto são distintos, então cada conjunto de fluxo recebe seu próprio Clifford no passo 5. As orientações inversas são Paulis diferentes, não os mesmos renomeados; é por isso que uma síntese 2D por conjuntos de fluxo precisa de quatro circuitos, enquanto a síntese 1D precisa de um.

4. Transpile com a codificação personalizada

Esta etapa é a versão 1D com um registrador de qubits mais largo: CustomF2QLayout mapeia o registrador de $N$ modos em $2N$ qubits, MapperFnEvolutionSynthesis carrega a codificação, e seu argumento product_formula seleciona como a PauliEvolutionGate resultante é decomposta. Compare com o passo 6 do guia 1D, que constrói o mesmo pipeline para um único ancilla.

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

Cinco portas, carregando os 6, 6, 6, 6 e 22 termos contados no passo 2. Cada uma é mapeada e sintetizada de forma independente, então a partição decidida no nível fermiônico sobrevive até o circuito:

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

Essa é a linha de base: a codificação e o agrupamento estão ambos no lugar, mas LieTrotter sintetiza cada conjunto de fluxo termo a termo, e a profundidade cresce com a rede.

5. Um Clifford de profundidade 2 por conjunto de fluxo

O mecanismo é o mesmo do caso 1D: conjugar um conjunto de fluxo inteiro de termos comutantes por um Clifford que leva cada um de seus termos a um Pauli distinto de peso 1, de modo que a evolução do conjunto inteiro custe apenas rotações de um único qubit. Em 1D, um brickwork de portas CZ fazia isso.

Em 2D, cada conjunto de fluxo precisa do seu próprio circuito, seguindo a Seção IV C e a Fig. 7 da Ref. [1]. Os quatro atuam em uma linha da rede por vez (uma linha para leste/oeste, uma coluna para norte/sul) junto com os ancillas dessa linha, e os quatro têm o mesmo formato em três partes:

• Uma camada de rotações de um único qubit, colocando cada qubit na base de que seu conjunto de fluxo precisa.

• Uma primeira camada de CX, sempre de sítio para ancilla no mesmo sítio.

• Uma segunda camada de CX, acoplando vizinhos ao longo da linha.

Os passos 2 e 3 são as únicas camadas de dois qubits, então a profundidade de dois qubits é $2$, independentemente do comprimento da linha. Elas são escritas abaixo, uma função por conjunto de fluxo. As rotações são como são porque foram encontradas, conjugando os termos de cada conjunto de fluxo e exigindo imagens distintas de peso 1, então trate-as como dadas e confie na verificação que se segue.

Os dois conjuntos horizontais são os mais simples. Seus termos têm peso 3, XXZ e YYZ, e diferem apenas pelo $S$ extra que transforma $X$ em $Y$ nos sítios:

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

Os dois conjuntos verticais têm peso 4, XYXY e YXXY, e são correspondentemente mais complexos: as rotações alternam com a posição ao longo da linha, e a segunda camada de CX também alterna, acoplando sítios em posições pares e ancillas em posições ímpares (ou o inverso):

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

Desenhado em uma única linha de três sítios (primeiro seus sítios, depois seus ancillas), todo o mecanismo cabe em uma só figura. O formato triangular é o mais claro dos dois: uma mudança de base nos ancillas, depois duas camadas de CX, a primeira emparelha cada sítio com seu próprio ancilla e a segunda alcança o próximo sítio ao longo da linha.

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

O formato quadrado precisa de mais trabalho de um único qubit (as rotações alternam ao longo da linha) e sua segunda camada de CX alterna entre acoplar sítios e acoplar ancillas:

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

Ambos têm duas camadas de CX de profundidade, e adicionar sítios à linha os alarga sem aprofundá-los; essa é toda a afirmação de profundidade constante, visível por inspeção. Os Cliffords oeste e norte diferem apenas em suas camadas de um único qubit, então não são desenhados separadamente.

O que a síntese ainda precisa é determinar onde cada termo vai parar. Uma linha de $L$ sítios carrega $L - 1$ ligações, cada uma contribuindo com um termo, então o auxiliar abaixo retorna uma entrada por ligação: o qubit em que o termo caiu e o Pauli de um único qubit em que ele se transformou. É crucial que isso seja lido da conjugação, e não previsto. Diferentemente do 1D, em que um brickwork de CZ mapeia $X_{j+1} \mapsto Z_j X_{j+1} Z_{j+2}$ em forma fechada, aqui não existe tal fórmula, e os ancillas de uma linha de coluna não são contíguos. A única coisa que difere entre os conjuntos de fluxo é a string de Pauli de entrada, que o passo 3 já imprimiu:

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

As imagens devem sair distintas e de peso 1, para que cada termo possa ser rotacionado de forma independente em seu próprio qubit; a profundidade de dois qubits de $2$ já é visível nas quatro construções acima, cada uma emitindo duas camadas de CX independentemente do comprimento da linha. Abaixo, ambas as propriedades são verificadas em todos os conjuntos de fluxo, todas as linhas e vários tamanhos de rede. O passo 6 então chama clifford_supports, de modo que ela é uma parte funcional da síntese, e não apenas uma demonstração.

6. Síntese ciente dos conjuntos de fluxo

A síntese ganha uma vantagem útil. Como Evolution se divide por grupo, synthesize é chamada uma vez por conjunto de fluxo (cinco vezes aqui, quatro conjuntos de salto mais o grupo diagonal), então todo termo de salto que ela vê compartilha um formato de Pauli, e um Clifford diagonaliza todos eles. Ela só precisa ler esse formato do primeiro termo e separar o resto pela linha em que estão. Então vem a observação que torna duas dimensões mais baratas do que se poderia temer:

Importante

Os componentes de um único conjunto de fluxo são disjuntos em vértices: as linhas do conjunto oeste não compartilham qubits. Assim, todos os Cliffords de um conjunto de fluxo podem ser emitidos juntos e compartilham as mesmas duas camadas de dois qubits, em vez de ficarem na fila uma linha após a outra. Um conjunto de fluxo custa, portanto, um número fixo de camadas (duas para o Clifford e duas para o seu inverso), não importa quantas linhas ou colunas a rede tenha, e é por isso que a profundidade total acaba constante.

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

Como em 1D, isto é escrito deliberadamente para esta codificação: ela reconhece os quatro formatos de Pauli e levanta um erro diante de qualquer outro, em vez de recorrer silenciosamente a uma alternativa. Ela também assume a divisão de um conjunto de fluxo por porta descrita acima, e diz isso; mais de um formato de salto em uma única porta significa que o agrupamento não foi aplicado, o que é um erro, e não algo a contornar.

Nota

Como sua equivalente 1D, esta síntese não recebe um número de passos de Trotter. Um passo de Trotter precisa alternar entre os conjuntos de fluxo, mas transpile_trotter_step abaixo decompõe o circuito fermiônico em uma porta Evolution por conjunto, então este método só vê um conjunto de cada vez e só poderia repetir dentro dele, o que não muda nada: repetir um operador apenas divide sua própria evolução em partes iguais.

A Trotterização de vários passos precisa, portanto, de uma fórmula de produto que respeite os groups do operador dentro de uma única porta, antes da decomposição. FermionicLieTrotter faz isso, e o passo 10 o usa para medir a convergência.

Observe as duas classes de nomes parecidos em jogo. O LieTrotter do Qiskit, usado acima como linha de base, atua sobre o Hamiltoniano mapeado, onde o agrupamento já se perdeu. FermionicLieTrotter é a contrapartida fermiônica, e é a que ainda consegue enxergar os conjuntos de fluxo.

7. O resultado: profundidade constante

Tudo até aqui rodou na única rede $3 \times 3$ construída no passo 1. Comparar os dois caminhos de síntese em redes crescentes exige repetir toda essa sequência (operador, agrupamento, codificação, circuito, gerenciador de passes) para cada tamanho, então vale encapsulá-la uma vez:

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

Coletar as duas métricas em uma faixa de tamanhos de rede e plotá-las contra a aresta da rede deixa a diferença evidente:

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

A profundidade de dois qubits é 24 e constante de $3 \times 3$ a $8 \times 8$ (ou seja, de 18 qubits até 128), enquanto LieTrotter aplicado ao Hamiltoniano codificado de forma idêntica cresce linearmente com a aresta da rede, de 47 a 92. A contagem de portas cresce nos dois casos, já que há $O(N)$ termos a aplicar de qualquer forma; o que os conjuntos de fluxo removem é o crescimento no número de camadas que essas portas ocupam.

O orçamento se divide em 16 para o salto e oito para a interação:

• Cada um dos quatro conjuntos de fluxo custa $2 + 2$ camadas de dois qubits: seu Clifford na entrada e o inverso na saída, ambos de profundidade dois por construção, com apenas rotações de um único qubit no meio. Isso dá $4 \times 4 = 16$, que é o valor que a Ref. [1] cita, já que ela conta apenas os termos de salto. Definir $U = 0$ no código acima o reproduz.

• O grupo diagonal adiciona os oito restantes. Suas portas Rzz se empacotam em quatro camadas, e cada Rzz se decompõe em duas portas CX em torno de uma rotação de um único qubit, então $4 \times 2 = 8$. Quatro camadas é o ótimo: em uma rede quadrada, a interação alcança apenas vizinhos mais próximos, então todo sítio interior carrega quatro portas Rzz (uma por ligação) e só pode participar de uma por vez.

Nada em nenhuma das duas somas se refere ao tamanho da rede, e esse é todo o ponto.

8. Verifique o circuito

A síntese por conjuntos de fluxo e LieTrotter descrevem o mesmo passo de Trotter do mesmo Hamiltoniano codificado, então seus unitários devem concordar até uma fase global. Vetores de estado aleatórios detectam qualquer discrepância de forma muito mais barata do que construir o unitário completo:

Nota

Ambos os circuitos são decompostos até portas de um e dois qubits antes de serem aplicados a um estado. Se não decompostas, cada instrução de alto nível Evolution força Statevector a construir seu operador denso completo $2^{2N} \times 2^{2N}$; decompor primeiro permite aplicar as portas uma de cada vez, o que aqui é várias centenas de vezes mais rápido.

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

Nota

As duas metades dessa verificação importam. Uma sobreposição de módulo unitário em um único estado é necessária, mas não suficiente: uma síntese que errasse uma fase dependente do estado ainda passaria. Exigir a mesma fase em vários estados aleatórios é o que fecha a lacuna.

9. Verifique a física no subespaço estabilizador

Esta é a complicação bidimensional. Em 1D, a codificação era uma isometria, e comparar espectros com uma referência de Jordan-Wigner bastava. O VC não é uma isometria: ele usa $2N$ qubits para $N$ modos, então um espaço de dimensão $2^{2N}$ abriga um físico de dimensão $2^N$. Isso é genérico para codificações locais em duas dimensões (os qubits excedentes são restringidos por simetrias locais comutantes, uma por plaqueta [2] [3]) e significa que o Hamiltoniano codificado representa o fermiônico fielmente apenas no autoespaço conjunto $+1$ dessas restrições.

Aviso

A palavra estabilizador tem duplo uso neste assunto, e os dois significados não têm relação. A Ref. [1] a usa para os operadores de transferência comutantes de um conjunto de fluxo, que formam um grupo estabilizador que os Cliffords do passo 5 diagonalizam (um artifício de compilação, e o sentido usado em tudo acima). Esta seção a usa no sentido de correção de erros: operadores cujo autoespaço $+1$ é o subespaço físico. A Ref. [1] assinala explicitamente essa colisão e deixa o segundo significado de lado; as restrições de plaqueta abaixo são uma propriedade da codificação VC [2], e não um resultado do arcabouço de conjuntos de fluxo.

Os geradores são construídos aqui em vez de citados: um por face quadrada, e cada um é simplesmente o laço de operadores de aresta codificados em torno dessa face.

Os operadores de aresta são o objeto natural para isso. Um laço fechado deles é uma simetria de qualquer Hamiltoniano que conserva o número de partículas, enquanto o laço correspondente de operadores de transferência não é: os dois diferem por um operador de vértice no sítio de origem de cada ligação,

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

Ao longo de um laço fechado, todo sítio é a origem de exatamente uma ligação, então o laço coleta um $V_j$ por sítio, que é a paridade local dos quatro sítios da plaqueta. Essa paridade local é o que faz a diferença: o Hamiltoniano move partículas através da fronteira da plaqueta, então conserva apenas a paridade global. O laço de transferência puro, portanto, não comuta com ela, enquanto o laço de aresta (que carrega o fator de paridade) comuta.

Como a codificação é definida sobre operadores de transferência, mapear um operador de aresta exige a relação acima. Isso é um callback de cinco linhas, e ele se encaixa em map_edge_vertex_generators() exatamente como vc_action se encaixa em 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))]

Cada gerador tem peso $6$, distribuído por dois dos qubits físicos de sua plaqueta e quatro ancillas. Eles são involutórios e comutam com o Hamiltoniano codificado, então o subespaço físico é preservado pela dinâmica. Ambas as propriedades são afirmações sobre operadores, então são verificadas simbolicamente na álgebra de SparseObservable, sem exigir matrizes $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

Com os estabilizadores em mãos, a comparação não precisa de correspondência de bases, porque $V_j = Z_j$ significa que as ocupações vivem nos mesmos qubits físicos nas duas visões. Então: escreva um estado produto diretamente nos qubits físicos com os ancillas em $|0\rangle$, projete-o no subespaço estabilizador, evolua com o Hamiltoniano codificado e leia as densidades dos sítios $\langle n_j \rangle = (1 - \langle Z_j \rangle)/2$ nos dois lados.

Cada passo aqui permanece no nível de um vetor. O projetor é aplicado um fator $(1 + S)/2$ de cada vez em vez de ser multiplicado por inteiro, e a evolução temporal usa expm_multiply(), que calcula $e^{-iHt} |\psi\rangle$ sem nunca formar $e^{-iHt}$. Como o Hamiltoniano codificado é um operador esparso de dimensão $2^{2N}$, essa distinção é o que mantém esta seção barata.

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

Essas são as densidades codificadas; elas coincidem com a referência de Jordan-Wigner até a precisão de máquina, em uma rede cuja codificação só é fiel após a projeção.

10. Meça a convergência de Trotter

O passo 8 estabeleceu que as duas sínteses concordam entre si, o que é uma afirmação sobre equivalência, e não sobre precisão: ambas descrevem o mesmo passo de Trotter, então ambas carregam o mesmo erro de Trotter. Medir esse erro exige uma referência exata, e o passo 9 acabou de construir uma.

O número de passos vem de FermionicLieTrotter, a contrapartida fermiônica da fórmula em nível de qubits usada como linha de base antes. Aplicada ao operador agrupado inteiro antes da decomposição, sua repetição alterna entre os quatro conjuntos de fluxo e o grupo diagonal, em vez de dentro de qualquer um deles:

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

A fidelidade sobe monotonicamente em direção a um à medida que os passos encurtam, que é o comportamento que uma fórmula de produto deve mostrar. Encurtar os passos é apenas um dos dois botões; FermionicSuzukiTrotter os troca um pelo outro em uma ordem mais alta, como mede Sintetize uma porta Evolution no espaço fermiônico.

Um único passo deste Hamiltoniano em $t = 1$ é uma aproximação ruim (os quatro conjuntos de salto e a interação não comutam), então a primeira linha fica longe de um, e esse é o ponto: a convergência é o que confirma a construção, não uma linha isolada.

Nota

Cada fator do cronograma abrange todos os modos desta rede, então restringir um fator ao seu suporte não tem efeito aqui e os índices absolutos de qubit com que FlowSetSynthesis trabalha permanecem inalterados. Uma síntese escrita para um operador cujos grupos tocam apenas parte do registrador teria de considerar essa renumeração.

O que levar daqui

• Duas dimensões mudam a contabilidade, não a receita. Todo ingrediente do guia 1D se mantém sem modificação: groups para o agrupamento, map_transfer_vertex_generators() para a codificação, CustomF2QLayout para o registrador incompatível e MapperFnEvolutionSynthesis para ligá-lo ao transpiler. O que cresceu foi o número de casos dentro de cada um.

• A disjunção é o que torna o 2D de profundidade constante. Um conjunto de fluxo se decompõe em linhas disjuntas em vértices, então seus Cliffords são emitidos em paralelo e compartilham duas camadas de dois qubits, não importa o tamanho da rede. Quatro conjuntos de fluxo e um grupo diagonal dão então um orçamento fixo: aqui $16$ para os quatro conjuntos de fluxo mais $8$ para o grupo de interação diagonal.

• Cada orientação é seu próprio formato de Pauli. $T_{kj} = -V_j V_k T_{jk}$ transforma XXZ em YYZ e XYXY em YXXY, então os quatro conjuntos de fluxo precisam de quatro Cliffords. Supor que um Clifford cobre as duas orientações de um eixo é a maneira mais fácil de errar uma síntese 2D por conjuntos de fluxo.

• Leia os suportes a partir da conjugação. Em 1D, um brickwork de CZ dá as imagens em forma fechada. Em 2D não há tal fórmula, e as linhas de coluna têm ancillas não contíguos, então um palpite derivado à mão sobre onde cada termo cai é frágil. Conjugar e inspecionar o resultado custa uma construção de Clifford e não pode estar errado; isso também significa que um Clifford igualmente válido que leve os termos para outro lugar continua funcionando.

• Deixe o agrupamento fazer a separação. Como Evolution se divide por grupo, a síntese recebe um conjunto de fluxo de cada vez e nunca precisa rederivar a partição: ela lê um formato de Pauli do primeiro termo e escolhe o Clifford correspondente. Assumir isso, e levantar um erro quando não vale, é o que mantém synthesize curta.

• A contagem de passos fica acima da divisão. A síntese aqui vê um conjunto de fluxo por porta, então não consegue subdividir a evolução entre os conjuntos. FermionicLieTrotter faz isso no espaço fermiônico, antes da decomposição, o que torna possível um teste de convergência.

• Verifique no subespaço. Uma codificação de $N$ modos em $2N$ qubits só é fiel onde os estabilizadores dizem que é. Como o VC mantém $V_j = Z_j$, a verificação confiável mais barata é projetar um estado produto e comparar densidades: sem intertwiner e sem correspondência espectral, que não é confiável quando o Hamiltoniano restrito é degenerado.

Nota

O Hamiltoniano codificado comuta com os estabilizadores, então a projeção só precisa ser aplicada ao estado inicial. A evolução preserva o subespaço por si só. Em hardware real, isso também é o que torna a codificação útil como recurso de detecção de erros: os estabilizadores podem ser medidos.

Próximos passos

Duas direções estendem o que está acima, e ambas são mudanças de ingrediente, não de método. A primeira é a codificação: Verstraete-Cirac foi escolhida aqui porque a Ref. [1] dá seus Cliffords diagonalizadores em forma fechada, mas Derby-Klassen gasta um ancilla por plaqueta ímpar em vez de um por sítio, para uma razão qubit-modo de $3/2$ em vez de $2$. Seus operadores de transferência têm formatos de Pauli diferentes, então uma portagem equivale a reescrever a função de codificação do passo 3 e encontrar o Clifford para cada novo formato; o agrupamento, o layout e o arcabouço de síntese se mantêm intactos.

A segunda é a própria rede. Tudo aqui assume fronteiras abertas e uma única espécie fermiônica. Fronteiras periódicas fecham cada linha em um ciclo: as setas ainda comutam todas, já que cada uma encontra a próxima cabeça com cauda, mas os dois formatos de Clifford do passo 5 são construídos para uma linha com extremos. Um ciclo de $L$ sítios carrega $L$ ligações por orientação onde uma linha carrega $L - 1$, então a ligação de fechamento não tem lugar na construção e precisa de tratamento próprio. Um modelo com spin é mais fácil: dobra o registrador e adiciona termos de interação no sítio ao grupo diagonal, sem tocar os conjuntos de fluxo. Nenhum dos dois perturba o argumento de profundidade constante, mas ambos mudam a contabilidade, que, como diz o primeiro ponto, é a parte que duas dimensões tornam maior.

Importante

A conectividade do hardware não é tratada aqui. Todo circuito acima é construído sobre um registrador abstrato em que qualquer par de qubits pode interagir, e as profundidades relatadas são contadas nesse cenário. Dispositivos reais têm um mapa de acoplamento restrito, e mapear esta codificação em um deles é um problema substancial por si só: o registrador VC é um grafo de $2N$ qubits cujas portas ocorrem entre um sítio e seu próprio ancilla, entre sítios vizinhos e entre ancillas vizinhos, então um layout precisa manter os três tipos de pares adjacentes simultaneamente. Se não consegue, o transpiler insere portas SWAP, e as camadas que elas adicionam são o que a construção de profundidade constante vinha compensando. A Ref. [1] também cita seu valor de profundidade sob um grafo de conectividade quadrada assumido, com CX como único emaranhador nativo. Encontrar bons layouts e estratégias de roteamento para codificações locais em uma dada topologia de dispositivo está além do escopo deste guia.

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