Rotar modos bosónicos

Una ModeRotation mezcla modos bosónicos con una matriz unitaria $U$, como lo hace una red de divisores de haz y desfasadores. Se llama pasiva porque conserva el número total de bosones: mueve bosones entre modos y nunca crea ni destruye ninguno. (El estrujamiento y el desplazamiento son las operaciones activas, y una rotación de modos no puede expresarlas.) Mapea los operadores de creación como

$$b^\dagger_i \mapsto \sum_j U_{ji} \, b^\dagger_j \, ,$$

que en todo el espacio de Fock es el unitario de muchos cuerpos

$$\exp\left(\sum_{ij} \log(U)_{ij} \, b^\dagger_i b_j\right) \, .$$

Esta es la puerta que hay detrás de un interferómetro de óptica lineal: introduce estados de Fock, rota los modos y lee cómo se redistribuyen los bosones. Como toda puerta de un BosonicCircuit (consulta Trabajar con circuitos bosónicos), no lleva ni truncamiento ni codificación. Ambos entran en el momento de la transpilación, y esta guía muestra por qué el truncamiento requiere un cuidado especial en una rotación.

Construye una rotación

El constructor recibe la matriz unitaria de $M \times M$. La puerta actúa sobre $M$ modos. La clase no comprueba que la matriz sea unitaria, así que eso es responsabilidad tuya.

>>> import numpy as np
>>> from qiskit_bosons.circuit.library import ModeRotation
>>>
>>> unitary = np.array([[1, -1], [1, 1]]) / np.sqrt(2)
>>> rotation = ModeRotation(unitary)
>>> rotation.num_modes  # the gate acts on two modes
2

Las rotaciones habituales tienen constructores dedicados. ModeRotation.beamsplitter() recibe un ángulo de mezcla $\theta$ y una fase opcional $\phi$, y $\theta = \pi/4$ da el divisor de haz equilibrado 50:50. ModeRotation.phase_shift() recibe una fase por modo.

>>> beamsplitter = ModeRotation.beamsplitter(np.pi / 4)
>>> np.round(beamsplitter.rotation_unitary.real, 4)
array([[ 0.7071, -0.7071],
       [ 0.7071,  0.7071]])
>>> bool(np.allclose(beamsplitter.rotation_unitary, unitary))
True
>>> shift = ModeRotation.phase_shift([0.0, np.pi / 2])
>>> np.round(np.diag(shift.rotation_unitary), 4)
array([1.+0.j, 0.+1.j])

También puedes partir del generador antihermítico $K$ y dejar que la puerta lo exponencie, lo que da la rotación $\exp(K)$. El generador se recupera con generator().

>>> generator = np.array([[0, -0.3], [0.3, 0]], dtype=complex)
>>> small = ModeRotation.from_generator(generator)
>>> np.round(small.rotation_unitary.real, 4)
array([[ 0.9553, -0.2955],
       [ 0.2955,  0.9553]])
>>> bool(np.allclose(small.generator(), generator))
True

El puente hacia el álgebra de operadores es to_operator(), que devuelve el generador hermítico $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ como un BosonOperator, de modo que la puerta es $e^{-\mathrm{i} H}$. Para el divisor de haz 50:50, eso es un par de términos de salto con coeficientes imaginarios opuestos. El operador es hermítico, que es lo que necesita una puerta Evolution, y conserva el número de partículas.

>>> operator = beamsplitter.to_operator()
>>> # The stored order of the terms is not specified, so sort them before printing.
>>> sorted((term[0], round(term[1].imag, 4)) for term in operator.iter_terms())
[([(True, 0), (False, 1)], -0.7854), ([(True, 1), (False, 0)], 0.7854)]
>>> operator.is_hermitian()
True
>>> operator.conserves_particle_number()
True

Cada término se lee (is_creation, mode) de izquierda a derecha, así que el segundo es $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$.

Actúa sobre estados de Fock

Una rotación mueve bosones entre modos, así que su efecto se ve más fácilmente sobre estados de Fock. El efecto Hong-Ou-Mandel es el ejemplo estándar. Envía un bosón a cada modo de entrada de un divisor de haz equilibrado. Los dos caminos que dejan los bosones en modos de salida distintos interfieren destructivamente, de modo que el estado $\lvert 1, 1 \rangle$ desaparece y los bosones siempre salen juntos, en $\lvert 2, 0 \rangle$ o $\lvert 0, 2 \rangle$, cada uno con probabilidad un medio.

Construye el circuito sobre modos bosónicos y redúcelo con el gestor de pases preestablecido. El LevelSpec fija en este punto el truncamiento y la codificación. Un corte de max_occupation=2 conserva los niveles 0, 1 y 2 de cada modo.

>>> from qiskit_bosons.circuit import BosonicCircuit
>>> from qiskit_bosons.circuit.library import InitializeModes
>>> from qiskit_bosons.mappers import LevelEncoding, LevelSpec
>>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
>>>
>>> circuit = BosonicCircuit(2)
>>> circuit.append(InitializeModes([1, 1]), circuit.modes)
>>> circuit.append(ModeRotation.beamsplitter(np.pi / 4), circuit.modes)

El flujo devuelve un circuito de qubits ordinario. Lee las probabilidades de los resultados a partir del estado final, aquí con el Statevector de Qiskit, que es exacto a este tamaño. Sirve cualquier simulador que devuelva amplitudes o recuentos. Un estado de Fock $\lvert n_0, n_1 \rangle$ es un único estado base computacional, cuyo índice viene dado por los patrones de bits de LevelEncoding.level_to_bits(); los qubits del segundo modo siguen a los del primero.

>>> from qiskit.quantum_info import Statevector
>>> from qiskit.transpiler import PassManager
>>>
>>> def fock_probabilities(circuit, max_occupation):
...     """Maps each pair of occupations to the probability of that outcome."""
...     encoding = LevelEncoding.UNARY_UNPROJECTED
...     spec = LevelSpec(max_occupation, encoding)
...     pass_manager = generate_preset_occupation_pass_manager(spec)
...     pass_manager.qubit = PassManager()  # keep the evolutions as the pipeline emitted them
...     qubit_circuit = pass_manager.run(circuit)
...     state = Statevector(qubit_circuit)
...     width = encoding.num_qubits(max_occupation + 1)
...     probabilities = {}
...     for n0 in range(max_occupation + 1):
...         for n1 in range(max_occupation + 1):
...             index = encoding.level_to_bits(n0) | (encoding.level_to_bits(n1) << width)
...             probability = round(float(abs(state.data[index]) ** 2), 4)
...             if probability > 0:
...                 probabilities[(n0, n1)] = probability
...     return probabilities
>>>
>>> fock_probabilities(circuit, max_occupation=2)
{(0, 2): 0.5, (2, 0): 0.5}

El resultado $\lvert 1, 1 \rangle$ tiene probabilidad cero, y los dos resultados agrupados comparten la probabilidad a partes iguales.

Elige el corte

El resultado anterior necesita max_occupation=2. Ejecuta el mismo circuito con un corte de 1. Es suficiente para contener el estado de entrada $\lvert 1, 1 \rangle$, pero no el estado que la rotación produce a partir de él, que tiene peso en $\lvert 2, 0 \rangle$ y $\lvert 0, 2 \rangle$:

>>> fock_probabilities(circuit, max_occupation=1)
{(1, 1): 1.0}

El circuito se ejecuta sin ningún error, y la respuesta es errónea: el estado permanece en $\lvert 1, 1 \rangle$ con probabilidad uno. Una rotación mezcla sectores de ocupación dentro de un modo. Dos bosones pueden acabar en el mismo modo de salida, así que el estado exacto tiene peso en la ocupación 2 en un modo que empezó con un bosón. Con un corte de 1 ese nivel no existe, y el generador truncado no puede llevar el estado hasta allí.

Por tanto, una cota de la ocupación del estado inicial no basta. El corte debe cubrir el número total de bosones que pueden llegar a cualquier modo de la rotación. Para $N$ bosones en total, la opción segura es un corte de $N$.

El mismo razonamiento se aplica a cada modo por separado, de modo que combinar cortes distintos por modo (consulta Describe el truncamiento con un LevelSpec) hace fácil mapear incorrectamente una rotación. Un modo que puede recibir bosones necesita un corte que cubra lo que puede recibir, no solo lo que contiene inicialmente. La comprobación LevelSpec.truncation() compara una especificación con la cota inferior de ocupación de un operador. No puede señalar el caso anterior. Cada término del generador del divisor de haz contiene un único operador de creación, así que la cota es uno y la comprobación no informa de ningún déficit con un corte de 1, aunque el resultado allí es erróneo:

>>> LevelSpec(1, "unary_unprojected").truncation(beamsplitter.to_operator())
[0, 0]

Otra precaución se deriva de la manera en que el transpilador construye una rotación (siguiente sección). La aplica como una secuencia de factores de dos modos, y la documentación del complemento advierte de que un único factor puede poblar un nivel que la rotación completa nunca alcanza. La secuencia truncada y la evolución única truncada son entonces aproximaciones distintas del mismo unitario. Una búsqueda sobre 1080 rotaciones aleatorias de tres modos, con un bosón en cada uno de dos modos de entrada y sin doble ocupación en la salida, no encontró ningún caso así, de modo que trátalo como una precaución y no como un fallo frecuente. Haz converger el corte como en cualquier otro circuito, aumentándolo hasta que el resultado deje de cambiar.

Cómo sintetiza el transpilador una rotación

La etapa optimization del flujo preestablecido descompone primero una ModeRotation en el espacio bosónico, con el pase B2BSynthesis y el método BosonicClementsDecomposition (registrado como ModeRotation.Clements). Es una elección del preajuste y no un requisito del marco: una sección posterior reduce la misma puerta sin ninguna descomposición. La descomposición no exponencia el generador completo en un solo paso. Primero factoriza el unitario en como máximo $M(M-1)/2$ rotaciones entre modos adyacentes más $M$ fases de un solo modo, siguiendo la malla rectangular de Clements et al. [1]. Cada factor se convierte en una puerta Evolution atómica, que la etapa de bosón a qubit mapea como cualquier otra.

La factorización en sí es givens_decomposition(), a la que llama la descomposición. La aplicación que lleva una unitaria a la rotación de muchos cuerpos es un homomorfismo de grupos, de modo que el producto de las rotaciones de los factores es exactamente igual a la rotación completa en el espacio sin truncar. Por tanto, la malla no añade ningún error propio. Su ventaja es la localidad. Cada evolución emitida involucra solo los términos de Pauli de los dos modos que mezcla (o de un modo en el caso de una fase), en lugar de términos de todos los modos.

Transpila una rotación aleatoria de tres modos e inspecciona lo que emitió el pipeline. Una unitaria genérica de $3 \times 3$ necesita el máximo de tres rotaciones, más tres fases, y cada factor es su propio PauliEvolutionGate. Todos están declarados sobre el registro completo, pero los operadores internos actúan solo sobre unos pocos qubits:

>>> from qiskit.quantum_info import random_unitary
>>>
>>> random = np.asarray(random_unitary(3, seed=7).data)
>>> spec = LevelSpec(2, "unary_unprojected")
>>> rotation_circuit = BosonicCircuit(3)
>>> rotation_circuit.append(InitializeModes([1, 0, 0]), rotation_circuit.modes)
>>> rotation_circuit.append(ModeRotation(random), rotation_circuit.modes)
>>> pass_manager = generate_preset_occupation_pass_manager(spec)
>>> pass_manager.qubit = PassManager()
>>> qubit_circuit = pass_manager.run(rotation_circuit)
>>> dict(qubit_circuit.count_ops())
{'PauliEvolution': 6, 'x': 3}
>>> for instruction in qubit_circuit.data:
...     if instruction.operation.name == "PauliEvolution":
...         terms = instruction.operation.operator.to_sparse_list()
...         print(len(terms), sorted({qubit for _, indices, _ in terms for qubit in indices}))
64 [0, 1, 2, 3, 4, 5]
64 [3, 4, 5, 6, 7, 8]
64 [0, 1, 2, 3, 4, 5]
2 [1, 2]
2 [4, 5]
2 [7, 8]

Con tres niveles por modo, cada modo ocupa tres qubits. Un divisor de haz entre dos modos afecta a los seis qubits de esos modos, y una fase afecta a los qubits de un modo.

Los factores solo son exactos cuando cada evolución se sintetiza bien. Los términos del generador de un divisor de haz no conmutan, así que la fórmula de producto predeterminada de un solo paso deja un error en cada factor. Pasa una fórmula de producto con más repeticiones al plugin de Evolution que asigna los factores para reducirlo:

>>> from qiskit.passmanager import MultiStagePassManager
>>> from qiskit.synthesis import LieTrotter
>>> from qiskit_bosons.circuit.library.synthesis import BosonicClementsDecomposition
>>> from qiskit_bosons.mappers.library import occupation_encoding
>>> from qiskit_bosons.transpiler import (
...     BosonicCircuitToDAG,
...     BosonicPassManager,
...     QuantumDAGToCircuit,
... )
>>> from qiskit_bosons.transpiler.passes import (
...     B2BSynthesis,
...     B2QSynthesis,
...     MapperFnEvolutionSynthesis,
...     OccupationB2QLayout,
...     OccupationInitializeModesSynthesis,
... )
>>>
>>> def with_product_formula(product_formula):
...     synthesis = B2QSynthesis()
...     synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(
...         lambda operator, num_qubits: occupation_encoding(operator, spec, num_qubits),
...         product_formula,
...     )
...     synthesis.methods["InitializeModes"] = OccupationInitializeModesSynthesis(spec)
...     return MultiStagePassManager(
...         input=BosonicCircuitToDAG(),
...         optimization=BosonicPassManager(
...             B2BSynthesis({"ModeRotation": BosonicClementsDecomposition()})
...         ),
...         layout=OccupationB2QLayout(spec),
...         synthesis=synthesis,
...         output=QuantumDAGToCircuit(),
...     )

Un solo bosón que entra en el modo 0 sale en el modo $j$ con probabilidad $|U_{j0}|^2$. Compara el circuito sintetizado con ese resultado exacto para un número creciente de repeticiones. Descompón primero el circuito para que se aplique la fórmula de producto; sin ello, un simulador exacto como Statevector aplica la puerta como su exponencial completa y oculta el error.

>>> encoding = LevelEncoding.UNARY_UNPROJECTED
>>> width = encoding.num_qubits(3)
>>> exact = np.abs(random[:, 0]) ** 2
>>>
>>> def one_boson_probabilities(qubit_circuit):
...     """Probability that the single boson sits in mode j, for each j."""
...     state = Statevector(qubit_circuit)
...     probabilities = []
...     for j in range(3):
...         occupations = [1 if mode == j else 0 for mode in range(3)]
...         index = sum(
...             encoding.level_to_bits(n) << (width * mode) for mode, n in enumerate(occupations)
...         )
...         probabilities.append(abs(state.data[index]) ** 2)
...     return np.array(probabilities)
>>>
>>> for repetitions in (1, 4, 16):
...     lowered = with_product_formula(LieTrotter(reps=repetitions)).run(rotation_circuit)
...     error = np.max(np.abs(one_boson_probabilities(lowered.decompose(reps=2)) - exact))
...     print(f"{repetitions:2d} repetitions: error {error:.1e}")
 1 repetitions: error 1.2e-03
 4 repetitions: error 7.6e-05
16 repetitions: error 4.7e-06

El error disminuye en un factor de aproximadamente cuatro por cada cuadruplicación de las repeticiones, como se espera de una fórmula de primer orden. El coste crece en proporción, porque cada factor se repite.

Reducir una rotación sin descomponerla

Una descomposición en el espacio bosónico es una forma de reducir una rotación. (La biblioteca no incluye otra forma, pero la etapa de bosón a qubit no depende de ella.) B2QSynthesis busca un plugin por el nombre de la puerta, y un ModeRotation que se dejó entero le llega como cualquier otra puerta. Un plugin para él asigna de una vez el generador de ModeRotation.to_operator() y emite un único PauliEvolutionGate a tiempo unitario. Eso es todo:

>>> from qiskit.circuit.library import PauliEvolutionGate
>>> from qiskit_bosons.transpiler.passes.utils import map_node_single_register
>>>
>>> class DirectModeRotationSynthesis:
...     """Maps the generator of the whole rotation, with no bosonic decomposition."""
...
...     def __init__(self, mapper_fn):
...         self.mapper_fn = mapper_fn
...
...     def run(self, in_node, out_dag, *, b2q_layout):
...         modes, qreg = map_node_single_register(in_node, b2q_layout)
...         generator = in_node.op.to_operator().relabel_modes(list(modes))
...         observable = self.mapper_fn(generator, len(qreg))
...         out_dag.apply_operation_back(PauliEvolutionGate(observable, time=1.0), tuple(qreg))

Regístralo con el nombre de la puerta en un pass manager que no la descomponga:

>>> def mapper_fn(operator, num_qubits):
...     return occupation_encoding(operator, spec, num_qubits)
>>>
>>> synthesis = B2QSynthesis()
>>> synthesis.methods["Evolution"] = MapperFnEvolutionSynthesis(mapper_fn)
>>> synthesis.methods["InitializeModes"] = OccupationInitializeModesSynthesis(spec)
>>> synthesis.methods["ModeRotation"] = DirectModeRotationSynthesis(mapper_fn)
>>> direct = MultiStagePassManager(
...     input=BosonicCircuitToDAG(),
...     layout=OccupationB2QLayout(spec),
...     synthesis=synthesis,
...     output=QuantumDAGToCircuit(),
... )
>>> direct_circuit = direct.run(rotation_circuit)
>>> dict(direct_circuit.count_ops())
{'x': 3, 'PauliEvolution': 1}
>>> terms = direct_circuit.data[-1].operation.operator.to_sparse_list()
>>> len(terms), sorted({qubit for _, indices, _ in terms for qubit in indices})
(198, [0, 1, 2, 3, 4, 5, 6, 7, 8])

Ahora la rotación es una sola evolución sobre los nueve qubits, mientras que la malla tiene seis, ninguna de las cuales abarca más de seis. Las dos no son la misma aproximación. Aquí el truncamiento se aplica una sola vez al generador de la rotación completa, y en la malla a cada factor por separado (de ahí la advertencia sobre el corte más arriba). El camino directo es igual a occupation_encoding() del operador con cualquier corte, así que es la referencia con la que comparar una malla. Su coste es la localidad: sus términos abarcan todos los modos, y una fórmula de producto tiene más de ellos que subdividir.

Ninguna es mejor en general, y esa es la razón para mantener abierta la elección. Escribe un plugin como este para comparar reducciones, para usar una rotación específica de la codificación o para entregar la puerta a algo que la entienda de forma nativa.

Siguientes pasos

• Transpilar circuitos bosónicos explica las etapas que ejecuta el pass manager predefinido y cómo reemplazar el plugin que gestiona una puerta.

• Asignar representaciones de operadores cubre el truncamiento, las codificaciones y el LevelSpec que fija la disposición de cada modo.

• Trabajar con circuitos bosónicos muestra cómo ensamblar y combinar circuitos bosónicos.

• La documentación de la API enumera qiskit_bosons.circuit.library para las puertas, qiskit_bosons.linalg para la descomposición y qiskit_bosons.transpiler.passes para los plugins de síntesis.

[1]

W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016), arXiv:1603.08788.