Simula la dinámica de Fermi-Hubbard 2D con conjuntos de flujo

Importante

Los conceptos de esta guía están disponibles actualmente solo en la API de Python. Se ofrecerá una funcionalidad equivalente en la API de C en una versión futura.

Véase también

Lee primero la guía de conjuntos de flujo 1D. Esta guía es su secuela bidimensional y no repite la información desarrollada allí: operadores de transferencia y de vértice, qué es un conjunto de flujo, cómo groups sobrevive en el circuito y cómo se integra una codificación personalizada en el transpilador mediante MapperFnEvolutionSynthesis.

Consulta la guía de transpilación para las etapas del transpilador a las que se hace referencia a lo largo del texto, y la guía de mapeadores para la receta general de escritura de un mapeador personalizado.

La guía 1D terminaba mencionando que hacen falta codificaciones locales al trabajar con dos dimensiones. Esta guía da ese paso, para el modelo de Fermi-Hubbard bidimensional en una red cuadrada, usando la codificación de Verstraete-Cirac (VC) [2].

Hay tres diferencias con el problema 1D, que forman la base de este trabajo:

• Los conjuntos de flujo se multiplican. Una cadena tiene dos orientaciones de flecha; una red cuadrada tiene cuatro (este, oeste, norte, sur), y cada una es una unión de caminos dirigidos disjuntos en vértices (uno por fila o columna de la red) y no un único camino.

• La codificación gana un auxiliar por sitio, no uno para todo el sistema, lo que da una

  razón de qubits por modo de $2$. Los operadores de transferencia pasan a ser Paulis de peso 3 (horizontales) y de peso 4 (verticales), y cada uno de los cuatro conjuntos de flujo se mapea a su propia forma de Pauli, por lo que cada uno necesita su propio Clifford diagonalizador.

• El hamiltoniano codificado solo es fiel en un subespacio: el autoespacio conjunto $+1$

  de las restricciones de plaqueta que lleva toda codificación local 2D [2] [3]. En 1D la codificación era una isometría simple y la verificación era una comprobación del espectro. Aquí, la comparación debe proyectarse, como ilustra el paso 9.

El beneficio tiene la misma forma que en 1D: la profundidad de dos qubits de un paso de Trotter pasa a ser constante en el tamaño de la red, en $24$ ($16$ para el salto, coincidiendo con la Ref. [1], más $8$ para la interacción), mientras que una trotterización término a término del mismo hamiltoniano codificado crece linealmente. El número de compuertas sigue creciendo con el número de enlaces, por supuesto; lo que pasa a ser independiente del tamaño es cuántas capas ocupan esas compuertas. Las profundidades se cuentan siempre sobre un registro abstracto totalmente conectado; mapear la codificación sobre el mapa de acoplamiento restringido de un dispositivo es un problema aparte, que se trata al final.

La red y sus cuatro orientaciones

El primer cambio es visible en el grafo de interacción dirigido. Los sitios se numeran en orden por filas, site(r, c) = r * cols + c, y cada enlace lleva ambas orientaciones, como la cadena en 1D, pero ahora los enlaces recorren dos ejes:

La regla de conmutación es la de 1D, sin cambios: dos operadores de transferencia que comparten un sitio conmutan cuando las flechas fluyen a través de él y anticonmutan cuando chocan. Sin embargo, un sitio tiene ahora hasta cuatro enlaces incidentes en lugar de dos, así que las flechas que pasan por él pueden chocar de más maneras, y un conjunto de flechas que nunca chocan ya no es un único camino. Leer las flechas horizontales de izquierda a derecha da uno de esos conjuntos, las verticales de arriba abajo dan otro, e invertir cada uno da dos más. Estos son los conjuntos de flujo este, oeste, norte y sur que organizan todo lo que sigue.

La codificación de Verstraete-Cirac

Verstraete-Cirac (VC) empareja cada modo fermiónico $j$ con su propio qubit auxiliar. Con $N$ sitios de red, el modo $j$ vive en el qubit $j$ y su auxiliar está en el qubit $N + j$, de modo que el registro contiene $2N$ qubits. Los operadores de vértice se mantienen tan simples como pueden ser,

$$V_j = Z_j \, ,$$

que es lo que hace que las ocupaciones se puedan leer directamente más adelante. Los operadores de transferencia dependen de la orientación de la arista. Para una arista horizontal $j \to k$ con $j$ el sitio occidental, y una arista vertical con $j$ el sitio septentrional,

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

un “triángulo” de peso 3 y un “cuadrado” de peso 4. Identificar cada flecha importa: con $j$ como sitio occidental, $j \to k$ apunta al este, pero con $j$ como sitio septentrional, apunta al sur. Las orientaciones inversas se deducen de la identidad $T_{kj} = -V_j V_k T_{jk}$, que convierte el triángulo XXZ en YYZ, y el cuadrado XYXY en YXXY. Esas cuatro etiquetas distintas son los cuatro conjuntos de flujo, y el resto de esta guía los sigue.

Nota

Estas no son exactamente la Ec. (A1) de la Ref. [1], que establece «hasta rotaciones de un solo qubit» y difiere en tres aspectos que conviene nombrar, ya que un lector que consulte el artículo lo notará:

• El prefactor: la Ec. (A1) no lleva ninguno. El $\tfrac12$ de aquí es lo que hace que $T_{jk}^2 = \tfrac14$, de acuerdo con la normalización $T_{jk} = -\tfrac12 (a^\dagger_j - a_j)(a^\dagger_k + a_k)$ que usa TransferVertexOperator. El artículo solo necesita sus operadores de transferencia dentro de $\exp(-i\,dt\,H)$, donde un factor de escala se absorbe en el ángulo de rotación; aquí tienen que ser operadores que satisfagan una identidad.

• La orientación: la Ec. (A1) selecciona las letras según la paridad de la fila de $j$, que en una arista horizontal es la misma para ambos extremos, por lo que no distingue $T_{jk}$ de $T_{kj}$. El artículo recupera la orientación por separado, conjugando con una capa de $R^Z(\pi/2)$. Esta guía, en cambio, la incorpora a la codificación, derivando la inversa a partir de $T_{kj} = -V_j V_k T_{jk}$, que es lo que hace que las dos orientaciones caigan en Paulis distintos en el paso 3.

• La forma vertical: la Ec. (A1) alterna XXXX y YYYY según la paridad de la fila, mientras que esta guía usa un único XYXY. Son dos maneras de lograr lo mismo: una asimetría de letras que hace que una arista vertical choque correctamente con las aristas horizontales que confluyen en ella. No interpretes la diferencia como una libre elección de base. Eliminar la alternancia de paridad y usar un único XXXX rompe el álgebra por completo, por lo que las letras distintas se comprueban en lugar de darse por supuestas.

Nota

El ancilla se asocia a un sitio, lo que da la razón $2$; la codificación de Derby-Klassen es la que coloca un qubit auxiliar por plaqueta (impar), para obtener la razón menor $3/2$. Aquí se usa VC porque sus operadores de transferencia son aquellos cuyos Cliffords diagonalizadores da la Ref. [1] en forma cerrada.

1. El hamiltoniano en una red cuadrada

La construcción es la misma que en 1D, pero con una lista de aristas bidimensional. Los sitios conservan la numeración por filas usada en la figura anterior, site(r, c) = r * cols + c, de modo que los enlaces horizontales conectan índices consecutivos y los enlaces verticales conectan índices que difieren en 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. Cuatro conjuntos de flujo en lugar de dos

Cada término de salto se clasifica según la dirección en que apunta su flecha en la red. Comparar las coordenadas de fila/columna divmod de los dos extremos separa las aristas horizontales de las verticales, y el signo de la diferencia determina la orientación. Todo lo diagonal (la identidad, los términos $V_j$ y $V_j V_k$) va a un quinto grupo:

Nota

Convención de orientación. Los nombres se toman simplemente de la figura anterior: una flecha es este u oeste según apunte a la derecha o a la izquierda a lo largo de una fila, y norte o sur según apunte hacia arriba o hacia abajo a lo largo de una columna. Como la numeración recorre las filas, este y sur son las direcciones en las que el índice aumenta. Nada de lo que sigue depende de los nombres (solo etiquetan grupos), pero es fácil confundirlos, así que conviene fijarlos con la imagen. La Ref. [1] nombra sus conjuntos de flujo de la red cuadrada de la misma manera.

>>> 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 un conjunto de flujo, cada flecha o bien se encuentra con otra cabeza con cola, o bien no se encuentra con ninguna, como en 1D, por lo que un conjunto de flujo sigue conmutando término a término y su evolución no acarrea error de Trotter. Lo nuevo es que ahora un conjunto de flujo es disconexo: el conjunto este de una red de rows x cols consta de rows caminos independientes de izquierda a derecha. Esa desconexión es una ventaja y no una complicación, como muestra el paso 6.

3. Codificar el operador

La codificación es una única función que asigna un operador de transferencia generalizado a una cadena de Pauli, que se entrega a map_transfer_vertex_generators(). En comparación con la versión 1D, el cuerpo tiene más casos: diagonal, horizontal, vertical y las orientaciones inversas mediante $T_{kj} = -V_j V_k T_{jk}$.

Hay un problema. map_transfer_vertex_generators() invoca la función de retorno solo con un generador y nada más, mientras que decidir la cadena de Pauli de un enlace también requiere cols para recuperar las coordenadas de fila y columna a partir de un índice de modo. En lugar de anidar las funciones para capturarlo, functools.partial() lo enlaza, lo que mantiene ambas funciones en el nivel superior y llamables de forma independiente.

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

El sentido de agrupar antes de asignar se ve de inmediato: cada conjunto de flujo se codifica a una única forma de Pauli, uniforme en todo el 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']

Las cuatro formas de salto son distintas, por lo que cada conjunto de flujo obtiene su propio Clifford en el paso 5. Las orientaciones inversas son Paulis diferentes, no los mismos con otro nombre; por eso la síntesis de conjuntos de flujo en 2D necesita cuatro circuitos, mientras que la síntesis en 1D necesita uno.

4. Transpilar con la codificación personalizada

Esta etapa es la de 1D con un registro de qubits más ancho: CustomF2QLayout asigna el registro de $N$ modos a $2N$ qubits, MapperFnEvolutionSynthesis lleva la codificación, y su argumento product_formula selecciona cómo se descompone la PauliEvolutionGate resultante. Compara con el paso 6 de la guía 1D, que construye el mismo flujo para un solo 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 compuertas, que llevan los 6, 6, 6, 6 y 22 términos contados en el paso 2. Cada una se asigna y sintetiza de forma independiente, de modo que la partición decidida en el nivel fermiónico sobrevive hasta el 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)

Esa es la línea base: la codificación y la agrupación ya están en su lugar, pero LieTrotter sintetiza cada conjunto de flujo término a término, y la profundidad crece con la red.

5. Un Clifford de profundidad 2 por conjunto de flujo

El mecanismo es el mismo que en 1D: conjugar todo un conjunto de flujo conmutativo con un Clifford que asigna cada uno de sus términos a un Pauli distinto de peso 1, de modo que la evolución del conjunto completo cuesta solo rotaciones de un qubit. En 1D, un entramado de compuertas CZ lo lograba.

En 2D, cada conjunto de flujo necesita su propio circuito, siguiendo la Sección IV C y la Fig. 7 de la Ref. [1]. Los cuatro actúan sobre una línea de la red a la vez (una fila para este/oeste, una columna para norte/sur) junto con los ancillas de esa línea, y los cuatro tienen la misma forma en tres partes:

• Una capa de rotaciones de un solo qubit, que coloca cada qubit en la base que necesita su conjunto de flujo.

• Una primera capa de CX, siempre de sitio a ancilla en el mismo sitio.

• Una segunda capa de CX, que acopla vecinos a lo largo de la línea.

Los pasos 2 y 3 son las únicas capas de dos qubits, por lo que la profundidad de dos qubits es $2$, sin importar cuán larga sea la línea. Se escriben a continuación, una función por conjunto de flujo. Las rotaciones son las que son porque se encontraron conjugando los términos de cada conjunto de flujo y exigiendo imágenes distintas de peso 1, así que tómalas como dadas y confía en la comprobación que sigue.

Los dos conjuntos horizontales son los más simples. Sus términos son de peso 3, XXZ y YYZ, y difieren solo por la $S$ adicional que convierte $X$ en $Y$ en los sitios:

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

Los dos conjuntos verticales son de peso 4, XYXY y YXXY, y son correspondientemente más complejos: las rotaciones alternan con la posición a lo largo de la línea, y la segunda capa de CX también alterna, acoplando sitios en posiciones pares y ancillas en las impares (o al revés):

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

Dibujado sobre una única línea de tres sitios (primero sus sitios, luego sus ancillas), todo el mecanismo cabe en una sola imagen. La forma triangular es la más clara de las dos: un cambio de base en los ancillas, luego dos capas de CX; la primera empareja cada sitio con su propio ancilla y la segunda llega hasta el siguiente sitio a lo largo de la línea.

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

La forma cuadrada necesita más trabajo de un solo qubit (las rotaciones alternan a lo largo de la línea) y su segunda capa de CX alterna entre acoplar sitios y acoplar ancillas:

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

Ambas tienen una profundidad de dos capas de CX, y añadir sitios a la línea las ensancha sin hacerlas más profundas; esa es toda la afirmación de profundidad constante, visible por inspección. Los Cliffords oeste y norte difieren solo en sus capas de un qubit, por lo que no se dibujan por separado.

Lo que la síntesis aún necesita es determinar dónde termina cada término. Una línea de $L$ sitios tiene $L - 1$ enlaces, cada uno con un término, así que la función auxiliar de abajo devuelve una entrada por enlace: el qubit en el que cayó el término y el Pauli de un solo qubit en que se convirtió. Es crucial que esto se lea de la conjugación en lugar de predecirse. A diferencia de 1D, donde un entramado de CZ asigna $X_{j+1} \mapsto Z_j X_{j+1} Z_{j+2}$ en forma cerrada, aquí no existe tal fórmula, y los ancillas de una línea de columna no son contiguos. Lo único que difiere entre los conjuntos de flujo es la cadena de Pauli de entrada, que el paso 3 ya imprimió:

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

Las imágenes deben salir distintas y de peso 1, de modo que cada término pueda rotarse de forma independiente en su propio qubit; la profundidad de dos qubits de $2$ ya es visible en las cuatro construcciones anteriores, cada una de las cuales emite dos capas de CX independientemente de la longitud de la línea. A continuación se comprueban ambas propiedades en todos los conjuntos de flujo, todas las líneas y varios tamaños de red. El paso 6 llama entonces a clifford_supports, de modo que es una parte operativa de la síntesis y no solo una demostración.

6. Síntesis consciente de los conjuntos de flujo

La síntesis parte con una ventaja. Como Evolution se divide por grupos, synthesize se invoca una vez por conjunto de flujo (cinco veces aquí, cuatro conjuntos de salto más el grupo diagonal), de modo que todo término de salto que ve comparte una forma de Pauli y un único Clifford los diagonaliza a todos. Solo tiene que leer esa forma del primer término y ordenar el resto según la línea en que se encuentran. Luego viene la observación que hace que dos dimensiones sean más baratas de lo que se podría temer:

Importante

Los componentes de un único conjunto de flujo son disjuntos por vértices: las filas del conjunto oeste no comparten qubits. Por lo tanto, todos los Cliffords de un conjunto de flujo pueden emitirse juntos y comparten las mismas dos capas de dos qubits, en lugar de ponerse en cola una fila tras otra. Un conjunto de flujo cuesta así un número fijo de capas (dos para el Clifford y dos para su inverso) sin importar cuántas filas o columnas tenga la red, que es la razón por la que la profundidad total resulta 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 en 1D, esto está escrito deliberadamente para esta codificación: reconoce las cuatro formas de Pauli y lanza un error ante cualquier otra en lugar de recurrir silenciosamente a una alternativa. También supone la división de un conjunto de flujo por compuerta descrita arriba, y lo dice; más de una forma de salto en una sola compuerta significa que no se aplicó la agrupación, lo cual es un error y no algo que haya que sortear.

Nota

Al igual que su equivalente en 1D, esta síntesis no recibe un número de pasos de Trotter. Un paso de Trotter debe alternar entre los conjuntos de flujo, pero transpile_trotter_step, más abajo, descompone el circuito fermiónico en una compuerta Evolution por conjunto, de modo que este método solo ve un conjunto a la vez y solo podría repetir dentro de él, lo cual no cambia nada: repetir un operador solo divide su propia evolución en partes iguales.

Por tanto, la trotterización de varios pasos necesita una fórmula de producto que respete los groups del operador dentro de una sola compuerta, antes de la descomposición. FermionicLieTrotter lo hace, y el paso 10 lo usa para medir la convergencia.

Ten en cuenta las dos clases de nombre similar que intervienen. LieTrotter de Qiskit, usada como línea base arriba, actúa sobre el hamiltoniano asignado, donde la agrupación ya se ha perdido. FermionicLieTrotter es la contraparte fermiónica, y es la que todavía puede ver los conjuntos de flujo.

7. El resultado: profundidad constante

Todo lo anterior se ejecutó sobre la única red de $3 \times 3$ construida en el paso 1. Comparar las dos rutas de síntesis en redes crecientes requiere repetir toda esa secuencia (operador, agrupación, codificación, circuito, gestor de pases) por cada tamaño, así que conviene encapsularla una 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)

Recopilar ambas métricas para un rango de tamaños de red y graficarlas frente al borde de la red deja clara la diferencia:

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

La profundidad de dos qubits es 24 y constante desde $3 \times 3$ hasta $8 \times 8$ (es decir, de 18 a 128 qubits), mientras que LieTrotter aplicado al hamiltoniano codificado de forma idéntica crece linealmente con el borde de la red, de 47 a 92. El número de compuertas crece en ambos casos, ya que hay $O(N)$ términos que aplicar de cualquier manera; lo que eliminan los conjuntos de flujo es el crecimiento en el número de capas que ocupan esas compuertas.

El presupuesto se divide en 16 para el salto y ocho para la interacción:

• Cada uno de los cuatro conjuntos de flujo cuesta $2 + 2$ capas de dos qubits: su Clifford de entrada y el inverso de salida, ambos de profundidad dos por construcción, con solo rotaciones de un qubit en medio. Eso es $4 \times 4 = 16$, la cifra que cita la Ref. [1], ya que cuenta solo los términos de salto. Establecer $U = 0$ en el código anterior la reproduce.

• El grupo diagonal aporta las ocho restantes. Sus compuertas Rzz se empaquetan en cuatro capas, y cada Rzz se descompone en dos compuertas CX alrededor de una rotación de un qubit, de modo que $4 \times 2 = 8$. Cuatro capas es lo óptimo: en una red cuadrada, la interacción solo alcanza a los vecinos más próximos, así que cada sitio interior lleva cuatro compuertas Rzz (una por enlace) y solo puede participar en una a la vez.

Nada en ninguna de las dos sumas hace referencia al tamaño de la red, que es precisamente el punto.

8. Verificar el circuito

La síntesis por conjuntos de flujo y LieTrotter describen el mismo paso de Trotter del mismo hamiltoniano codificado, por lo que sus unitarios deben coincidir salvo una fase global. Los vectores de estado aleatorios detectan cualquier discrepancia con mucho menos coste que construir el unitario completo:

Nota

Ambos circuitos se descomponen en compuertas de uno y dos qubits antes de aplicarlos a un estado. Si se dejan sin descomponer, cada instrucción Evolution de alto nivel obliga a Statevector a construir su operador denso completo de $2^{2N} \times 2^{2N}$; descomponer primero le permite aplicar las compuertas una a una, lo que aquí es varios cientos de veces más 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

Ambas mitades de esa comprobación importan. Un solapamiento de módulo unitario en un solo estado es necesario pero no suficiente: una síntesis que errara en una fase dependiente del estado seguiría mostrándolo. Exigir la misma fase en varios estados aleatorios es lo que cierra la brecha.

9. Verificar la física en el subespacio estabilizador

Este es el inconveniente de las dos dimensiones. En 1D, la codificación era una isometría y bastaba comparar espectros con una referencia de Jordan-Wigner. VC no es una isometría: usa $2N$ qubits para $N$ modos, de modo que un espacio de dimensión $2^{2N}$ contiene uno físico de dimensión $2^N$. Esto es genérico de las codificaciones locales en dos dimensiones (los qubits sobrantes están restringidos por simetrías locales conmutativas, una por plaqueta [2] [3]) y significa que el hamiltoniano codificado representa fielmente al fermiónico solo en el autoespacio conjunto $+1$ de esas restricciones.

Advertencia

La palabra estabilizador tiene doble uso en este tema, y los dos significados no guardan relación. La Ref. [1] la usa para los operadores de transferencia conmutativos de un conjunto de flujo, que forman un grupo estabilizador que los Cliffords del paso 5 diagonalizan (un recurso de compilación, y el sentido que se ha usado en todo lo anterior). Esta sección la usa en el sentido de la corrección de errores: operadores cuyo autoespacio $+1$ es el subespacio físico. La Ref. [1] señala explícitamente esta colisión y deja de lado el segundo significado; las restricciones de plaqueta de abajo son una propiedad de la codificación VC [2], no un resultado del marco de conjuntos de flujo.

Los generadores se construyen aquí en lugar de citarse: uno por cada cara cuadrada, y cada uno es simplemente el lazo de operadores de arista codificados alrededor de esa cara.

Los operadores de arista son el objeto natural para esto. Un lazo cerrado de ellos es una simetría de cualquier hamiltoniano que conserve el número de partículas, mientras que el lazo correspondiente de operadores de transferencia no lo es: ambos difieren en un operador de vértice en el sitio de origen de cada enlace,

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

En un lazo cerrado cada sitio es el origen de exactamente un enlace, por lo que el lazo recoge un $V_j$ por sitio, que es la paridad local de los cuatro sitios de la plaqueta. Esa paridad local es lo que marca la diferencia: el hamiltoniano hace saltar partículas a través del borde de la plaqueta, así que conserva solo la paridad global. Por tanto, el lazo de transferencia desnudo no conmuta con él, mientras que el lazo de aristas (que lleva el factor de paridad) sí conmuta.

Como la codificación está definida sobre operadores de transferencia, asignar un operador de arista necesita la relación anterior. Es una función de retorno de cinco líneas, y encaja en map_edge_vertex_generators() exactamente igual que vc_action encaja en 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 generador tiene peso $6$, repartido entre dos de los qubits físicos de su plaqueta y cuatro ancillas. Son involutivos y conmutan con el hamiltoniano codificado, de modo que la dinámica preserva el subespacio físico. Ambas propiedades son afirmaciones sobre operadores, por lo que se comprueban simbólicamente sobre el álgebra de SparseObservable, sin necesidad de matrices de $2^{18} \times 2^{18}$:

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

Con los estabilizadores en mano, la comparación no necesita correspondencia de bases, porque $V_j = Z_j$ significa que las ocupaciones viven en los mismos qubits físicos en ambas descripciones. Así: se escribe un estado producto directamente en los qubits físicos con los ancillas en $|0\rangle$, se proyecta sobre el subespacio estabilizador, se evoluciona con el hamiltoniano codificado y se leen las densidades de sitio $\langle n_j \rangle = (1 - \langle Z_j \rangle)/2$ en ambos lados.

Cada paso aquí se mantiene a nivel de vector. El proyector se aplica un factor $(1 + S)/2$ a la vez en lugar de multiplicarse por completo, y la evolución temporal usa expm_multiply(), que calcula $e^{-iHt} |\psi\rangle$ sin formar nunca $e^{-iHt}$. Como el hamiltoniano codificado es un operador disperso de dimensión $2^{2N}$, esa distinción es lo que mantiene barata esta sección.

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

Esas son las densidades codificadas; coinciden con la referencia de Jordan-Wigner hasta la precisión de máquina, en una red cuya codificación solo es fiel tras la proyección.

10. Medir la convergencia de Trotter

El paso 8 estableció que las dos síntesis coinciden entre sí, lo cual es una afirmación sobre equivalencia y no sobre exactitud: ambas describen el mismo paso de Trotter, así que ambas acarrean el mismo error de Trotter. Medir ese error requiere una referencia exacta, y el paso 9 acaba de construir una.

El número de pasos proviene de FermionicLieTrotter, la contraparte fermiónica de la fórmula a nivel de qubit usada antes como línea base. Aplicada a todo el operador agrupado antes de la descomposición, su repetición alterna entre los cuatro conjuntos de flujo y el grupo diagonal en lugar de hacerlo dentro de uno solo de ellos:

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

La fidelidad se acerca monótonamente a uno a medida que los pasos se acortan, que es el comportamiento que debe mostrar una fórmula de producto. Acortar los pasos es solo una de las dos palancas; FermionicSuzukiTrotter las equilibra entre sí a un orden superior, como mide Sintetizar una compuerta Evolution en el espacio fermiónico.

Un solo paso de este hamiltoniano en $t = 1$ es una mala aproximación (los cuatro conjuntos de salto y la interacción no conmutan), por lo que la primera fila está lejos de uno, y ese es el punto: lo que confirma la construcción es la convergencia, no una fila aislada.

Nota

Cada factor del calendario abarca todos los modos de esta red, por lo que reducir un factor a su soporte no tiene efecto aquí y los índices absolutos de qubit con los que coincide FlowSetSynthesis no cambian. Una síntesis escrita para un operador cuyos grupos tocan solo parte del registro tendría que tener en cuenta ese reetiquetado.

Qué conclusiones sacar

• Dos dimensiones cambian la contabilidad, no la receta. Todos los ingredientes de la guía 1D se trasladan sin modificación: groups para la agrupación, map_transfer_vertex_generators() para la codificación, CustomF2QLayout para el registro desajustado y MapperFnEvolutionSynthesis para conectarlo al transpilador. Lo que creció es el número de casos dentro de cada uno.

• La disyunción es lo que hace que 2D tenga profundidad constante. Un conjunto de flujo se descompone en líneas disjuntas por vértices, de modo que sus Cliffords se emiten en paralelo y comparten dos capas de dos qubits sin importar cuán grande sea la red. Cuatro conjuntos de flujo y un grupo diagonal dan entonces un presupuesto fijo: aquí $16$ para los cuatro conjuntos de flujo más $8$ para el grupo de interacción diagonal.

• Cada orientación es su propia forma de Pauli. $T_{kj} = -V_j V_k T_{jk}$ convierte XXZ en YYZ y XYXY en YXXY, de modo que los cuatro conjuntos de flujo necesitan cuatro Cliffords. Suponer que un Clifford cubre ambas orientaciones de un eje es la forma más fácil de equivocarse en una síntesis de conjuntos de flujo en 2D.

• Lee los soportes a partir de la conjugación. En 1D, un entramado de CZ da las imágenes en forma cerrada. En 2D no existe tal fórmula, y las líneas de columna tienen ancillas no contiguos, por lo que una suposición derivada a mano de dónde cae cada término es frágil. Conjugar e inspeccionar el resultado cuesta una construcción de Clifford y no puede fallar; además, significa que un Clifford igualmente válido que ubique los términos en otro lugar sigue funcionando.

• Deja que la agrupación haga la clasificación. Como Evolution se divide por grupo, la síntesis recibe un conjunto de flujo a la vez y nunca tiene que volver a derivar la partición: lee una forma de Pauli del primer término y elige el Clifford correspondiente. Suponerlo, y lanzar un error cuando no se cumple, es lo que mantiene breve synthesize.

• El número de pasos queda por encima de la división. Aquí la síntesis ve un conjunto de flujo por compuerta, por lo que no puede subdividir la evolución entre los conjuntos. FermionicLieTrotter lo hace en el espacio fermiónico, antes de la descomposición, que es lo que hace posible una prueba de convergencia.

• Verifica en el subespacio. Una codificación de $N$ modos con $2N$ qubits es fiel solo donde los estabilizadores lo dicen. Como VC mantiene $V_j = Z_j$, la comprobación fiable más barata es proyectar un estado producto y comparar densidades: sin entrelazador ni comparación espectral, que no es fiable cuando el hamiltoniano restringido es degenerado.

Nota

El hamiltoniano codificado conmuta con los estabilizadores, por lo que la proyección solo tiene que aplicarse al estado inicial. La evolución preserva el subespacio por sí misma. En hardware real, esto es también lo que hace útil la codificación como recurso de detección de errores: los estabilizadores se pueden medir.

Siguiente

Dos direcciones amplían lo anterior, y ambas son cambios de ingrediente y no de método. La primera es la codificación: aquí se eligió Verstraete-Cirac porque la Ref. [1] da sus Cliffords diagonalizadores en forma cerrada, pero Derby-Klassen gasta un ancilla por plaqueta impar en lugar de uno por sitio, para una razón qubit-modo de $3/2$ en vez de $2$. Sus operadores de transferencia tienen formas de Pauli diferentes, de modo que un portado consiste en reescribir la función de codificación del paso 3 y hallar el Clifford de cada nueva forma; la agrupación, la disposición y el andamiaje de síntesis se conservan intactos.

La segunda es la propia red. Todo lo de aquí supone fronteras abiertas y una sola especie fermiónica. Las fronteras periódicas cierran cada línea en un ciclo: las flechas siguen conmutando todas, ya que cada una se encuentra con la siguiente cabeza con cola, pero las dos formas de Clifford del paso 5 están construidas para una línea con extremos. Un ciclo de $L$ sitios tiene $L$ enlaces por orientación, donde una línea tiene $L - 1$, así que el enlace de cierre no tiene lugar en la construcción y requiere un tratamiento propio. Un modelo con espín es más fácil: duplica el registro y añade términos de interacción en el sitio al grupo diagonal, sin tocar los conjuntos de flujo. Ninguno de los dos perturba el argumento de profundidad constante, pero ambos cambian la contabilidad, que, como dice la primera conclusión, es la parte que dos dimensiones agranda.

Importante

Aquí no se aborda la conectividad del hardware. Todos los circuitos anteriores se construyen sobre un registro abstracto en el que cualquier par de qubits puede interactuar, y las profundidades indicadas se cuentan en ese contexto. Los dispositivos reales tienen un mapa de acoplamiento restringido, y asignar esta codificación a uno es un problema sustancial por sí mismo: el registro VC es un grafo de $2N$ qubits cuyas compuertas actúan entre un sitio y su propio ancilla, entre sitios vecinos y entre ancillas vecinos, por lo que una disposición debe mantener adyacentes los tres tipos de pares simultáneamente. Si no puede, el transpilador inserta compuertas SWAP, y las capas que estas añaden son lo que la construcción de profundidad constante había compensado. La Ref. [1] también cita su cifra de profundidad bajo un grafo de conectividad cuadrada supuesto con CX como único entrelazador nativo. Encontrar buenas disposiciones y estrategias de enrutamiento para codificaciones locales en una topología de dispositivo dada queda fuera del alcance de esta guía.

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