Faire pivoter les modes bosoniques

Une ModeRotation mélange des modes bosoniques avec une matrice unitaire $U$, comme le fait un réseau de séparateurs de faisceau et de déphaseurs. Elle est dite passive parce qu'elle conserve le nombre total de bosons : elle déplace des bosons entre modes sans jamais en créer ni en détruire. (Le compression d'état (squeezing) et le déplacement sont les opérations actives, et une rotation de modes ne peut pas les exprimer.) Elle mappe les opérateurs de création ainsi :

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

ce qui, sur tout l'espace de Fock, est l'unitaire à plusieurs corps

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

C'est la porte derrière un interféromètre à optique linéaire : envoie des états de Fock en entrée, fais pivoter les modes et lis comment les bosons se redistribuent. Comme toute porte d'un BosonicCircuit (voir Travailler avec des circuits bosoniques), elle ne porte ni troncature ni encodage. Les deux interviennent au moment de la transpilation, et ce guide montre pourquoi la troncature exige un soin particulier pour une rotation.

Construire une rotation

Le constructeur prend la matrice unitaire $M \times M$. La porte agit sur $M$ modes. La classe ne vérifie pas que la matrice est unitaire ; c'est donc ta responsabilité.

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

Les rotations courantes disposent de constructeurs dédiés. ModeRotation.beamsplitter() prend un angle de mélange $\theta$ et une phase optionnelle $\phi$, et $\theta = \pi/4$ donne le séparateur de faisceau équilibré 50:50. ModeRotation.phase_shift() prend une phase par mode.

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

Tu peux aussi partir du générateur anti-hermitien $K$ et laisser la porte l'exponentier, ce qui donne la rotation $\exp(K)$. Le générateur est retrouvé avec 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

Le pont vers l'algèbre des opérateurs est to_operator(), qui renvoie le générateur hermitien $H = \sum_{ij} \mathrm{i} \log(U)_{ij} \, b^\dagger_i b_j$ sous forme de BosonOperator, de sorte que la porte est $e^{-\mathrm{i} H}$. Pour le séparateur de faisceau 50:50, il s'agit d'une paire de termes de saut avec des coefficients imaginaires opposés. L'opérateur est hermitien, ce dont une porte Evolution a besoin, et il conserve le nombre de particules.

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

Chaque terme se lit (is_creation, mode) de gauche à droite ; le second est donc $\mathrm{i} \tfrac{\pi}{4} \, b^\dagger_1 b_0$.

Agir sur des états de Fock

Une rotation déplace des bosons entre modes ; son effet se voit donc plus facilement sur des états de Fock. L'effet Hong-Ou-Mandel en est l'exemple standard. Envoie un boson dans chaque mode d'entrée d'un séparateur de faisceau équilibré. Les deux chemins qui laissent les bosons dans des modes de sortie différents interfèrent destructivement, de sorte que l'état $\lvert 1, 1 \rangle$ disparaît et que les bosons sortent toujours ensemble, dans $\lvert 2, 0 \rangle$ ou $\lvert 0, 2 \rangle$, chacun avec une probabilité d'un demi.

Construis le circuit sur des modes bosoniques et abaisse-le avec le pass manager prédéfini. Le LevelSpec fixe ici la troncature et l'encodage. Une coupure de max_occupation=2 conserve les niveaux 0, 1 et 2 de chaque mode.

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

Le pipeline renvoie un circuit de qubits ordinaire. Lis les probabilités des résultats à partir de l'état final, ici avec le Statevector de Qiskit, qui est exact à cette taille. Tout simulateur qui renvoie des amplitudes ou des comptages convient. Un état de Fock $\lvert n_0, n_1 \rangle$ est un unique état de base computationnelle, dont l'indice est donné par les motifs de bits de LevelEncoding.level_to_bits() ; les qubits du second mode suivent ceux du premier.

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

Le résultat $\lvert 1, 1 \rangle$ a une probabilité nulle, et les deux résultats groupés se partagent la probabilité à parts égales.

Choisir la coupure

Le résultat précédent exige max_occupation=2. Exécute le même circuit avec une coupure de 1. Elle suffit pour contenir l'état d'entrée $\lvert 1, 1 \rangle$, mais pas l'état que la rotation en produit, qui a du poids sur $\lvert 2, 0 \rangle$ et $\lvert 0, 2 \rangle$ :

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

Le circuit s'exécute sans aucune erreur, et la réponse est fausse : l'état reste en $\lvert 1, 1 \rangle$ avec une probabilité de un. Une rotation mélange des secteurs d'occupation au sein d'un mode. Deux bosons peuvent finir dans le même mode de sortie ; l'état exact a donc du poids sur l'occupation 2 dans un mode qui avait démarré avec un boson. Avec une coupure de 1, ce niveau n'existe pas, et le générateur tronqué ne peut pas y déplacer l'état.

Une borne sur l'occupation de l'état initial ne suffit donc pas. La coupure doit couvrir le nombre total de bosons pouvant atteindre un mode quelconque de la rotation. Pour $N$ bosons au total, le choix sûr est une coupure de $N$.

Le même raisonnement s'applique à chaque mode séparément ; un mélange de coupures différentes par mode (voir Décrire la troncature avec un LevelSpec) rend donc une rotation facile à mapper incorrectement. Un mode qui peut recevoir des bosons a besoin d'une coupure qui couvre ce qu'il peut recevoir, et pas seulement ce qu'il contient au départ. La vérification LevelSpec.truncation() compare un spec à la borne inférieure d'occupation d'un opérateur. Elle ne peut pas signaler le cas ci-dessus. Chaque terme du générateur du séparateur de faisceau contient un seul opérateur de création ; la borne est donc de un et la vérification ne signale aucun déficit pour une coupure de 1, alors que le résultat y est faux :

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

Une autre précaution découle de la manière dont le transpileur construit une rotation (section suivante). Il applique la rotation comme une séquence de facteurs à deux modes, et la documentation du plugin avertit qu'un seul facteur peut peupler un niveau que la rotation complète n'atteint jamais. La séquence tronquée et l'évolution unique tronquée sont alors deux approximations différentes de la même unitaire. Une recherche sur 1 080 rotations aléatoires à trois modes, avec un boson dans chacun de deux modes d'entrée et sans double occupation en sortie, n'a pas révélé un tel cas ; considère donc cela comme une précaution et non comme une défaillance courante. Fais converger la coupure comme pour n'importe quel autre circuit, en l'augmentant jusqu'à ce que le résultat cesse de changer.

Comment le transpileur synthétise une rotation

L'étape optimization du pipeline prédéfini décompose d'abord une ModeRotation dans l'espace bosonique, avec la passe B2BSynthesis et la méthode BosonicClementsDecomposition (enregistrée sous ModeRotation.Clements). C'est un choix du preset et non une exigence du framework : une section ultérieure abaisse la même porte sans aucune décomposition. La décomposition n'exponentie pas le générateur complet en une seule étape. Elle factorise d'abord l'unitaire en au plus $M(M-1)/2$ rotations entre modes adjacents plus $M$ phases à un mode, en suivant le maillage rectangulaire de Clements et al. [1]. Chaque facteur devient une porte Evolution atomique, que l'étape boson-vers-qubit mappe comme n'importe quelle autre.

La factorisation elle-même est givens_decomposition(), que la décomposition appelle. L'application qui associe une unitaire à la rotation à N corps est un homomorphisme de groupes, donc le produit des rotations des facteurs est exactement égal à la rotation complète sur l'espace non tronqué. Le maillage n'ajoute donc aucune erreur propre. Son avantage est la localité. Chaque évolution émise n'implique que les termes de Pauli des deux modes qu'elle mélange (ou d'un seul mode pour une phase), au lieu de termes couvrant tous les modes.

Transpile une rotation aléatoire à trois modes et inspecte ce que le pipeline a émis. Une unitaire $3 \times 3$ générique nécessite au maximum trois rotations, plus trois phases, et chaque facteur est son propre PauliEvolutionGate. Tous sont déclarés sur le registre entier, mais les opérateurs qu'ils contiennent n'agissent que sur quelques 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]

Avec trois niveaux par mode, chaque mode occupe trois qubits. Un séparateur de faisceau entre deux modes touche les six qubits de ces modes, et une phase touche les qubits d'un seul mode.

Les facteurs ne sont exacts que si chaque évolution est bien synthétisée. Les termes du générateur d'un séparateur de faisceau ne commutent pas, donc la formule de produit par défaut en une seule étape laisse une erreur dans chaque facteur. Passe une formule de produit avec davantage de répétitions au plugin Evolution qui fait correspondre les facteurs, afin de la réduire :

>>> 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 boson unique entrant dans le mode 0 ressort dans le mode $j$ avec la probabilité $|U_{j0}|^2$. Compare le circuit synthétisé à ce résultat exact pour un nombre croissant de répétitions. Décompose d'abord le circuit afin que la formule de produit soit appliquée ; sans cela, un simulateur exact tel que Statevector applique la porte comme son exponentielle complète et masque l'erreur.

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

L'erreur diminue d'un facteur d'environ quatre pour chaque multiplication par quatre du nombre de répétitions, comme attendu pour une formule du premier ordre. Le coût augmente proportionnellement, car chaque facteur est répété.

Abaisser une rotation sans la décomposer

Une décomposition dans l'espace bosonique est une façon d'abaisser une rotation. (La bibliothèque n'en inclut pas d'autre, mais l'étape boson-vers-qubit n'en dépend pas.) B2QSynthesis recherche un plugin par le nom de la porte, et une ModeRotation laissée entière lui parvient comme n'importe quelle autre porte. Un plugin pour elle fait correspondre le générateur de ModeRotation.to_operator() en une seule fois et émet un unique PauliEvolutionGate à temps unitaire. C'est tout :

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

Enregistre-le pour le nom de la porte dans un gestionnaire de passes qui ne la décompose pas :

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

La rotation est maintenant une seule évolution sur les neuf qubits, là où le maillage en a six, dont aucune n'en couvre plus de six. Les deux ne sont pas la même approximation. La troncature est appliquée une seule fois au générateur de la rotation entière ici, et à chaque facteur séparément dans le maillage (d'où l'avertissement sur la coupure ci-dessus). Le chemin direct est égal à occupation_encoding() de l'opérateur pour toute coupure, c'est donc la référence à laquelle comparer un maillage. Son coût est la localité : ses termes couvrent tous les modes, et une formule de produit a davantage de termes à subdiviser.

Aucun n'est meilleur en général, ce qui justifie de garder le choix ouvert. Écris un plugin comme celui-ci pour comparer des abaissements, pour utiliser une rotation propre à un encodage, ou pour confier la porte à quelque chose qui la comprend nativement.

Étapes suivantes

• Transpiler des circuits bosoniques explique les étapes qu'exécute le gestionnaire de passes prédéfini, et comment remplacer le plugin qui traite une porte.

• Faire correspondre les représentations d'opérateurs couvre la troncature, les encodages, et la LevelSpec qui fixe la disposition de chaque mode.

• Travailler avec des circuits bosoniques montre comment assembler et combiner des circuits bosoniques.

• La documentation de l'API répertorie qiskit_bosons.circuit.library pour les portes, qiskit_bosons.linalg pour la décomposition, et qiskit_bosons.transpiler.passes pour les plugins de synthèse.

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