In December 2020, there was this claim of quantum advantage/supremacy by a team of UST China using Gaussian boson sampling. Here is the paper and here is an explanatory news article in Nature. To get an idea of how that works, I wanted to try to replicate a reduced version in Qiskit. Is this possible? Are there predefined functions I can exploit for doing photon-like circuits (if not how I can implement the boson-to-qubit mapping)?

The original paper on the theory of boson sampling by Scott Aaronson and Alex Arkhipov can be found here.

I looked in Qiskit documentation and I did not find any boson-qubit mapping (Holstein-Primakoff, Dyson-Maleev or Schwinger bosons). I just found out that there is a BosonicOperator in Qiskit Chemistry module but I do not know what is its purpose without mappings.

I would also like to have an estimation of the minimum number of qubits I will be needing. I am not trying to test quantum advantage. The paper reported a 76 photons coincidence which is a lot for even the best of supercomputers. The calculation involves permanents which are in practice exponentially difficult to calculate with ordinary computers. My goal is to simulate boson sampling to calculate the result of a very small permanent (if possible the smallest one) that I can verify later on my PC.

Note: according to Is it possible to "calculate" the absolute value of a permanent using Boson Sampling? you cannot calculate permanents directly using boson sampling.

  • 1
    $\begingroup$ +1. I asked a related question before: quantumcomputing.stackexchange.com/q/2244/2293, which I mention here since you said you wanted to calculate the permanent of a matrix. $\endgroup$ Nov 27, 2021 at 8:29
  • $\begingroup$ Were you able to find implementation of Boson Sampling on Qiskit ? Thanks. $\endgroup$ Aug 10, 2022 at 19:38
  • $\begingroup$ @ChetanWaghela not yet $\endgroup$
    – Mauricio
    Aug 11, 2022 at 9:13

1 Answer 1


I'm not sure if this is what you're looking for, but you can always describe the evolution of a many-boson system as the evolution of a set of qudits in the standard "kets and unitaries" formalism, at least provided the evolution is linear (as in, the total number of bosons is preserved). Probably also more generally, but I'm not sure.

Simple example

Consider a simple example. Two bosons in two modes, evolving with a unitary $U$. Say the initial state is $a_1^\dagger a_2^\dagger |0\rangle\equiv |1_1 1_2\rangle$, meaning one boson in each input mode. If the bosons evolve through the (single-particle) unitary $U$, then the overall output state is obtained evolving each creation operator via $U$, meaning applying the mapping $$a_i^\dagger \to \sum_j u_{ji} a_j^\dagger,$$ and thus reads $$|\operatorname{output}\rangle= (u_{11} a_1^\dagger + u_{21} a_2^\dagger)(u_{12} a_1^\dagger + u_{22} a_2^\dagger)|0\rangle = \Bigg[\sqrt2 u_{11} u_{12} \frac{(a_1^\dagger)^2}{\sqrt2} + (u_{12}u_{21}+u_{11}u_{22})a_1^\dagger a_2^\dagger + \sqrt2 u_{21} u_{22} \frac{(a_2^\dagger)^2}{\sqrt2} \Bigg] |0\rangle = \sqrt2 u_{11} u_{12} |2_1 0_2\rangle + (u_{12}u_{21}+u_{11}u_{22})|1_1 1_2\rangle + \sqrt2 u_{21} u_{22} |0_1 2_2\rangle.$$ But there's another way to describe similar situations. If we were not using the second quantisation formalism (that is, no creation/annihilation operators), then we would describe the initial state of two bosons in the two input modes with something like $|1,2\rangle$, where now the first (second) number tells you the state of the first (second) boson. But the bosons are indistinguishable aside for their mode state, thus we cannot say things like "the first boson is in the first mode and the second in the second mode", we have to say things like "one boson is in the first mode and the other one in the second mode". We can do this by symmetrising the input state, with something like $$\frac1{\sqrt2}(|1,2\rangle+|2,1\rangle).$$ Now each boson evolves independently with $U$, hence the output state reads $$(U\otimes U)\frac1{\sqrt2}(|1,2\rangle+|2,1\rangle).$$ You might now see that this allows you to simulate the situation with a two-qubit system. Upon redefinition of the labels, this is a Bell state evolving through a local unitary $U\otimes U$. There's however a catch, that you'll also want to perform non-local measurements on the output. This is because in boson sampling you want the output probabilities associated to the various occupation numbers, and the three (in this simple example) corresponding basis states read in this notation: $$|2_1 0_2\rangle \sim |1,1\rangle, \qquad |0_1 2_2\rangle \sim |2,2\rangle, \qquad |1_1 1_2\rangle \sim |1,2\rangle+|2,1\rangle.$$

You can verify explicitly that doing this produces exactly the same probability amplitudes the original description of the many-boson system did.

General statement

The above is fully general. Given an $n$-boson, $m$-mode state, you can describe it as a symmetrised $n$-qudit state with each qudit having dimension $m$. Linear evolution is then described by some $U\otimes\cdots\otimes U$, and the outcome probabilities for the many-boson system can be recovered projecting on the corresponding $n$-qudit states. This sort of thing you can do in a circuit, and I suppose if you can simulate $n$-qudit states with $nm$-qubit ones, you can build your circuit with qubits. Not that this will be easy, mind you. I'd expect the hardness of boson-sampling will translate into the hardness of creating the highly entangled states necessary to play this little game.

  • $\begingroup$ Thanks! This is very clear. Let me try this before validating the answer. $\endgroup$
    – Mauricio
    Aug 23, 2022 at 9:21
  • $\begingroup$ @Mauricio just as an aside though: you don't need to map the system into qubits to simulate it on your computer. You can just directly simulate the many boson system. The only nontrivial part is computing the various permanents, but if you use eg Mathematica you can find pretty efficient (relatively speaking) code to do it: mathematica.stackexchange.com/q/38177/27539 $\endgroup$
    – glS
    Aug 23, 2022 at 9:33
  • $\begingroup$ I wanted to try how good it run in a real backend. The solution works well, I have to think how this works on for higher number of qubits (not qudits) now. $\endgroup$
    – Mauricio
    Aug 23, 2022 at 10:14
  • $\begingroup$ What is the dimension of $U$ in the $n$ bosons $m$ modes case? $\endgroup$
    – Mauricio
    Aug 24, 2022 at 9:33
  • 1
    $\begingroup$ it's a unitary acting on a single qudit, so $m\times m$ $\endgroup$
    – glS
    Aug 24, 2022 at 9:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.