Rotacionar modos bosônicos

Uma ModeRotation mistura modos bosônicos com uma matriz unitária $U$, como faz uma rede de divisores de feixe e defasadores. Ela é chamada de passiva porque conserva o número total de bósons: move bósons entre modos e nunca cria nem destrói nenhum. (Squeezing e deslocamento são as operações ativas, e uma rotação de modos não consegue expressá-las.) Ela mapeia os operadores de criação como

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

o que, no espaço de Fock completo, é o unitário de muitos corpos

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

Esta é a porta por trás de um interferômetro de óptica linear: envie estados de Fock, rotacione os modos e leia como os bósons se redistribuem. Como toda porta de um BosonicCircuit (veja Trabalhar com circuitos bosônicos), ela não carrega truncamento nem codificação. Ambos entram no momento da transpilação, e este guia mostra por que o truncamento exige cuidado especial para uma rotação.

Construa uma rotação

O construtor recebe a matriz unitária $M \times M$. A porta atua em $M$ modos. A classe não verifica se a matriz é unitária, então isso é responsabilidade sua.

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

Rotações comuns têm construtores dedicados. ModeRotation.beamsplitter() recebe um ângulo de mistura $\theta$ e uma fase opcional $\phi$, e $\theta = \pi/4$ dá o divisor de feixe 50:50 balanceado. ModeRotation.phase_shift() recebe uma 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])

Você também pode partir do gerador anti-hermitiano $K$ e deixar a porta exponenciá-lo, o que dá a rotação $\exp(K)$. O gerador é recuperado com 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

A ponte para a álgebra de operadores é to_operator(), que retorna o gerador hermitiano $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ como um BosonOperator, de modo que a porta é $e^{-\mathrm{i} H}$. Para o divisor de feixe 50:50, isso é um par de termos de salto com coeficientes imaginários opostos. O operador é hermitiano, que é o que uma porta Evolution precisa, e conserva o 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 termo lê (is_creation, mode) da esquerda para a direita, então o segundo é $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$.

Atue em estados de Fock

Uma rotação move bósons entre modos, então seu efeito é mais fácil de ver em estados de Fock. O efeito Hong-Ou-Mandel é o exemplo padrão. Envie um bóson para cada modo de entrada de um divisor de feixe balanceado. Os dois caminhos que deixam os bósons em modos de saída diferentes interferem destrutivamente, então o estado $\lvert 1, 1 \rangle$ desaparece e os bósons sempre saem juntos, em $\lvert 2, 0 \rangle$ ou $\lvert 0, 2 \rangle$, cada um com probabilidade um meio.

Construa o circuito em modos bosônicos e reduza-o com o pass manager predefinido. O LevelSpec fixa o truncamento e a codificação neste ponto. Um corte de max_occupation=2 mantém os níveis 0, 1 e 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)

O pipeline retorna um circuito de qubits comum. Leia as probabilidades dos resultados a partir do estado final, aqui com o Statevector do Qiskit, que é exato neste tamanho. Qualquer simulador que retorne amplitudes ou contagens funciona. Um estado de Fock $\lvert n_0, n_1 \rangle$ é um único estado da base computacional, cujo índice é dado pelos padrões de bits de LevelEncoding.level_to_bits(); os qubits do segundo modo vêm depois dos do primeiro.

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

O resultado $\lvert 1, 1 \rangle$ tem probabilidade zero, e os dois resultados agrupados dividem a probabilidade igualmente.

Escolha o corte

O resultado anterior precisa de max_occupation=2. Execute o mesmo circuito com um corte de 1. Isso é suficiente para comportar o estado de entrada $\lvert 1, 1 \rangle$, mas não o estado que a rotação produz a partir dele, que tem peso em $\lvert 2, 0 \rangle$ e $\lvert 0, 2 \rangle$:

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

O circuito executa sem nenhum erro, e a resposta está errada: o estado permanece em $\lvert 1, 1 \rangle$ com probabilidade um. Uma rotação mistura setores de ocupação dentro de um modo. Dois bósons podem acabar no mesmo modo de saída, então o estado exato tem peso na ocupação 2 em um modo que começou com um bóson. Com um corte de 1, esse nível não existe, e o gerador truncado não consegue mover o estado para lá.

Um limite sobre a ocupação do estado inicial, portanto, não basta. O corte deve cobrir o número total de bósons que pode chegar a qualquer modo isolado da rotação. Para $N$ bósons no total, a escolha segura é um corte de $N$.

O mesmo raciocínio se aplica a cada modo separadamente, então uma mistura de cortes diferentes por modo (veja Descreva o truncamento com um LevelSpec) torna fácil mapear uma rotação de forma incorreta. Um modo que pode receber bósons precisa de um corte que cubra o que ele pode receber, e não apenas o que contém inicialmente. A verificação LevelSpec.truncation() compara uma spec com o limite inferior de ocupação de um operador. Ela não consegue sinalizar o caso acima. Cada termo do gerador do divisor de feixe contém um único operador de criação, então o limite é um e a verificação não relata nenhuma deficiência com um corte de 1, embora o resultado ali esteja errado:

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

Uma ressalva adicional decorre da maneira como o transpilador constrói uma rotação (próxima seção). Ele aplica a rotação como uma sequência de fatores de dois modos, e a documentação do plugin avisa que um único fator pode popular um nível que a rotação completa nunca alcança. A sequência truncada e a evolução única truncada são, então, aproximações diferentes do mesmo unitário. Uma busca em 1.080 rotações aleatórias de três modos, com um bóson em cada um de dois modos de entrada e sem ocupação dupla na saída, não encontrou nenhum caso assim, então trate isso como uma precaução e não como uma falha comum. Faça o corte convergir como você faria para qualquer outro circuito, aumentando-o até que o resultado pare de mudar.

Como o transpilador sintetiza uma rotação

O estágio optimization do pipeline predefinido decompõe uma ModeRotation primeiro no espaço bosônico, com o pass B2BSynthesis e o método BosonicClementsDecomposition (registrado como ModeRotation.Clements). Essa é uma escolha do preset e não uma exigência do framework: uma seção posterior reduz a mesma porta sem nenhuma decomposição. A decomposição não exponencia o gerador completo em um único passo. Primeiro ela fatora o unitário em no máximo $M(M-1)/2$ rotações entre modos adjacentes mais $M$ fases de modo único, seguindo a malha retangular de Clements et al. [1]. Cada fator se torna uma porta Evolution atômica, que o estágio de bóson para qubit mapeia como qualquer outra.

A fatoração em si é feita por givens_decomposition(), que a decomposição chama. O mapeamento de uma unitária para a rotação de muitos corpos é um homomorfismo de grupos, então o produto das rotações dos fatores é exatamente igual à rotação completa no espaço não truncado. Portanto, a malha não adiciona nenhum erro próprio. Seu benefício é a localidade. Cada evolução emitida envolve apenas os termos de Pauli dos dois modos que ela mistura (ou de um modo, no caso de uma fase), em vez de termos que atravessam todos os modos.

Transpile uma rotação aleatória de três modos e inspecione o que o pipeline emitiu. Uma unitária $3 \times 3$ genérica precisa do máximo de três rotações, mais três fases, e cada fator é seu próprio PauliEvolutionGate. Todos são declarados sobre o registrador inteiro, mas os operadores internos atuam apenas em alguns 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]

Com três níveis por modo, cada modo ocupa três qubits. Um beamsplitter entre dois modos toca os seis qubits desses modos, e uma fase toca os qubits de um modo.

Os fatores só são exatos quando cada evolução é bem sintetizada. Os termos de um gerador de beamsplitter não comutam, então a fórmula de produto padrão de um único passo deixa um erro em cada fator. Passe uma fórmula de produto com mais repetições ao plugin Evolution que mapeia os fatores para reduzi-lo:

>>> 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(),
...     )

Um único bóson que entra no modo 0 sai no modo $j$ com probabilidade $|U_{j0}|^2$. Compare o circuito sintetizado com esse resultado exato para contagens de repetição crescentes. Decomponha o circuito primeiro para que a fórmula de produto seja aplicada; sem isso, um simulador exato como o Statevector aplica a porta como sua exponencial completa e esconde o erro.

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

O erro cai por um fator de cerca de quatro a cada quadruplicação das repetições, como esperado para uma fórmula de primeira ordem. O custo cresce em proporção, porque cada fator é repetido.

Reduza uma rotação sem decompô-la

Uma decomposição no espaço bosônico é uma forma de reduzir uma rotação. (A biblioteca não inclui outra forma, mas a etapa de bóson para qubit não depende dela.) B2QSynthesis procura um plugin pelo nome da porta, e uma ModeRotation que foi deixada inteira chega até ele como qualquer outra porta. Um plugin para ela mapeia o gerador de ModeRotation.to_operator() de uma só vez e emite um único PauliEvolutionGate em tempo unitário. É só isso:

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

Registre-o para o nome da porta em um pass manager que não a decomponha:

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

A rotação agora é uma única evolução sobre todos os nove qubits, enquanto a malha tem seis, das quais nenhuma abrange mais de seis. As duas não são a mesma aproximação. A truncagem é aplicada uma vez ao gerador da rotação inteira aqui, e a cada fator separadamente na malha (daí o aviso sobre o corte acima). O caminho direto é igual a occupation_encoding() do operador em qualquer corte, então ele é a referência para comparar uma malha. O que ele custa é a localidade: seus termos abrangem todos os modos, e uma fórmula de produto tem mais deles para subdividir.

Nenhum dos dois é melhor em geral, e essa é a razão para manter a escolha em aberto. Escreva um plugin como este para comparar reduções, para usar uma rotação específica de codificação ou para entregar a porta a algo que a entenda nativamente.

Próximos passos

• Transpile circuitos bosônicos explica as etapas que o pass manager predefinido executa e como substituir o plugin que trata uma porta.

• Mapeie representações de operadores aborda a truncagem, as codificações e o LevelSpec que fixa o layout de cada modo.

• Trabalhe com circuitos bosônicos mostra como montar e combinar circuitos bosônicos.

• A documentação da API lista qiskit_bosons.circuit.library para as portas, qiskit_bosons.linalg para a decomposição e qiskit_bosons.transpiler.passes para os plugins de síntese.

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