Dynamical Casimir Effect for Gaussian Boson Sampling

We show that the Dynamical Casimir Effect (DCE), realized on two multimode coplanar waveg-uide resonators, implements a gaussian boson sampler (GBS). The appropriate choice of the mirror acceleration that couples both resonators translates into the desired initial gaussian state and many-boson interference in a boson sampling network. In particular, we show that the proposed quantum simulator naturally performs a classically hard task, known as scattershot boson sampling. Our result unveils an unprecedented computational power of DCE, and paves the way for using DCE as a resource for quantum simulation.

The Dynamical Casimir Effect (DCE) consists in the generation of photons out of the vacuum of a quantum field by means of the abrupt modulation of boundary conditions -e.g. a mirror oscillating at speeds comparable to the speed of light. Predicted in 1970 1 , an experimental demonstration remained elusive until 2011, when it was implemented in a superconducting circuit architecture 2 . In addition to its fundamental interest, it has been shown that the radiation generated in the DCE displays entanglement and other forms of quantum correlations [3][4][5][6] , which poses the question of its utility as a resource for the heralded quantum technological revolution. As an example, small-scale continuous variable cluster states of four electromagnetic field modes have been shown to be in principle possible 7 . While this represents a preliminary step for a continuous variable one-way quantum computer, its scalability has not yet been demonstrated, hence the usefulness of DCE for quantum computing tasks remains unclear.
In this work, we establish a bridge between multimode parametric amplification induced by the modulation of boundary conditions -for which DCE is a paradigmatic case-and a non-universal quantum computing device, known as boson sampling (BS) 8 .
BS has recently gained a great deal of attention, as it solves a tailor-made problem-the problem of sampling from the output distribution of photons in a linear-optics network-that is widely believed to be intractable in any classical device. Thus it represents a promising avenue for proving the long-sought quantum supremacy 9 . We consider Gaussian BS (GBS) and in particular scattershot BS, a generalization of the original BS problem which is known to be equivalent in terms of computational complexity 10 . GBS have been proven to be of practical interest in reconstructing the Franck-Condon profile -a central problem in molecular spectroscopy,-both theoretically 11 and in a recent experimental trapped-ion implementation 12 . We show that GBS can be implemented in a superconducting circuit architecture by exploiting the possibility of multimode parametric amplification by means of the modulation of boundary conditions. We propose a setup consisting of two superconducting resonators coupled through a superconducting quantum interferometric device (SQUID) 13 (see Fig. 1). The resonators possess different lengths and thus different energy spectra and the SQUID plays the role of a shared tunable mirror-like boundary condition. The modulation of the external magnetic flux threading the SQUID implements an effective motion of the mirror whose corresponding Bogoliubov transformation results in multimode parametric amplification. We show that suitable choices of the SQUID pumping are able to implement the operations of a GBS -namely two-mode squeezers, beam-splitters and phase shifters. In this way, we show how the DCE can be exploited as a quantum simulator of GBS. Moreover, we will discuss how DCE is by itself a physical effect that is hard to simulate on a classical computer.

Methods
The DCE was observed in an open microwave coplanar waveguide interrupted by a single SQUID operated well below its plasma frequency 2 . Under the latter condition, the SQUID implements an effective mirror-like boundary condition. Ultrafast variation of the magnetic flux threading the SQUID amounts to motion of the mirror at Scientific RepoRts | (2018) 8:3751 | DOI:10.1038/s41598-018-22086-2 relativistic speeds, which generates a two-mode squeezing operation on the microwave field propagating along the transmission line. In particular, for an initial vacuum field state the modulation of the boundary condition results in generation of pairs of photons, a process which is resonantly enhanced if the mirror moves at a frequency matching the sum of the frequencies of the emitted photons.
Obviously, the DCE can be produced as well for different boundary conditions, such as the ones of a superconducting resonator interrupted by one 14 or two 15 SQUIDs. Moreover, we can think of the DCE as a particular instance of multimode parametric amplification induced by the modulation of boundary conditions, as we shall see in the following.
Let us consider a one-dimensional (1D) superconducting resonator in the presence of one or two movable boundary conditions. In the absence of any flux modulation, the resonator field φ is characterised by a set of creation and annihilation operators † where the Bogoliubov coefficients {α jl , β jl } are given by the inner product: jl j l jl j l Therefore, they depend on the particular initial boundary conditions -which determine the functions u l -and the particular type of boundary modulation -which determine the functions v l . In the case of small boundary oscillations characterised by a dimensionless amplitude δ, the Bogoliubov coefficients can be computed perturbatively in a variety of a cases including single 17 -and two-wall oscillations 18 . More general continuous motion of the two walls, such as the one required to mimic the motion of an accelerated cavity which is rigid in its proper frame, can also be addressed perturbatively 19 . In all these cases, the Bogoliubov coefficients depend on the features of the boundary modulation, for instance, the number of external pumps, together with their corresponding frequencies and durations.
Notice that the set {α jl } characterizes phase-shifting (j = l) and beam-splitting (j ≠ l) among the modes 20 , while {β jl } generates two-mode squeezing 2 . Therefore, we conclude that the modulation of the boundary conditions of a superconducting resonator is equivalent to a multimode parametric amplifier consisting of a set of tunable phase shifters, beam splitters and two-mode squeezers, which can be adjusted by suitably selecting the number, frequency and duration of external pumps.
A remarkable example is the DCE, where a modulation of frequency ω p generates Bogoliubov coefficients β lj growing linearly in time only for the modes in a resonance condition ω p = ω l + ω j -all the other Bogoliubov coefficients being negligible. Similarly, another resonance frequency ω p = |ω l − ω j | would make the corresponding α lj to increase linearly in time 19 . Clearly, a series of operations of this kind with suitable frequencies and duration times can generate a desired combination of beam splitters and two-mode squeezers.

Results
In the following, we will show how to use this scheme to implement a GBS protocol (see Fig. 2a) 10 .
An important problem is the well known lack of addressability in a harmonic oscillator. We need to generate beam splitters interactions only between a selected pair of modes, without further driving unwanted transition between modes. A way to overcome this problem is by considering a pair of resonators of different lengths -and hence different energy spectrum-confining two resonator fields φ r , φ l and sharing a tunable mirror (see Fig. 1a). In particular, to make sure this addressability condition holds for any pair of modes, it is convenient to choose resonators of incommensurate lengths.This condition can be relaxed, if the we set the appropriate cutoff in the energy spectrum of the resonators.
Therefore the collection of modes u consists of two collections of modes u l , u r given by the solutions to the Klein-Gordon 1D equation -plane waves-after imposing Dirichlet boundary conditions at points, say −L, 0 for u l and 0, L′ for u r : and v is the propagation speed of the field. The collection of modes v is given by the transformed solutions of the field confined in both resonators when the common mirror -initially placed at x = 0-undergoes an effective harmonic motion of frequency ω p and amplitude A, given by is a dimensionless small parameter. Then we can obtain the Bogoliubov coefficients as a perturbative expansion in δ. In particular, the first order of the expansion will be given by 19 : where ω l , ω′ j are the frequencies of the modes u l and v j respectively and α jl are the Bogoliubov coefficients associated to the transformation induced by a single change of mirror position x = δL′ at t = 0. Then: Now, if we select a mode u l in the left cavity and a mode v r in the right cavity: j r j i t 1 j we find that: By inspection of Eq. (4), we see that Bogoliubov coefficients oscillate in time, unless the pumping frequency ω p is either In the former case the corresponding α jl contain a term that grows monotonically in time, while in the latter the same happens for β jl . In both alternative scenarios, after a time ω ±  t 1 p , we can neglect all the oscillations, so we can assume that only the resonant terms are non-zero. In the first case (by simplicity, we assume L′ > L): while in the second one we obtain exactly the same expression for β jl , with a minus sign on the front. Notice that in our setup, the difference or sum of two mode frequencies when the modes correspond to different resonators is different for each pair. Thus Eq. (10) and the corresponding expression for the β jl illustrate the ability of implementing a beam splitter or a two-mode squeezer between a pair of selected modes by means of a suitable choice of the pumping frequency ω p . Then, the magnitude of the Bogoliubov coefficients can be adjusted with the choice of the amplitude A and duration t of the oscillation. In order to obtain random coefficients, we can randomize the value of A by means of a random number generator. Note that, while obviously α jl and β jl cannot be complex with the approach above-although additional random phases could be added by means of a rotation of the pump 2 -, it has been shown that arbitrary random real matrices are still classically hard to sample, as long as the entries are both positive and negative numbers 8,21 .

DCE for Gaussian Boson Sampling. Circuit QED implementation. Gaussian Boson Sampling
and scattershot boson sampling. GBS 10 is a modification of the original BS problem, where the initial state is Gaussian, as opposed to the initial state of the original BS, which is a Fock state. An example of Gaussian initial state would be a product of several two-mode squeezed states. In particular, we can think of a setup in which half of the output of n 2 two-mode squeezers are input into a linear network of n 2 optical modes, while the other half is sent directly to single photon detectors. Then, n single photons are detected in the latter half. As shown in 10 this device is able to solve a randomized version of BS known as scattershot BS, which possesses similar computational complexity as the original problem 8 , and therefore is widely believed to be out of reach classically in the same limit, namely approximately 20 single photons and 400 modes. Therefore, the necessary optical operators are: (I) two-mode squeezers for the preparation of the initial state, (II), beam splitters and phase shifters for the implementation of the linear network 22 and (III), single photon detectors.
Circuit QED implementation. The results of the previous section suggest that superconducting circuits are particularly suitable for implementing the above operators. In particular, we propose multimode coupled resonators as laid out in Fig. 1b, to prepare two-mode squeezed states between arbitrary modes i, j, and realize beam-splitting operations among the modes. This is done by time-evolving the coupled resonators under circuit Hamiltonian The charges q k are canonical conjugate variables of the flux modes φ k , which have frequencies ω k = kπv/L, where L is the size of the resonator and = v lc 1/ the speed of light in the coplanar waveguide-l and c are the inductance and capacitance per unit length, respectively. Naturally, a similar expression holds for modes φ′ k . E J± is the Josephson energy of the junctions in the SQUID loop, and ϕ ± represents the superconducting phase drop across each junction. The second line of Eq. (11) accounts for the nonlinear inductive energy of the SQUID, and is responsible of the coupling between resonators. This is obvious after imposing the fluxoid quantization relation ext 0 , being Φ 0 the flux quantum. Introducing the boson creation and annihilation operators, a k , † a k , and assuming small phase slips ϕ − across the SQUID, the circuit Hamiltonian can be written in the interaction picture as † † lj lj 1 0 -where J 1 is the Bessel function of the first kind-results from Jacobi-Anger expansion, which is indeed the Bogoliubov coefficient β lj given by equation (10). The sequential evolution under (13) for distinct pair of modes i, j and the different frequencies ω + p results in the preparation of the desired initial state |ψ〉 = S † |ψ 0 〉, where |ψ 0 〉 = |0〉 ⊗2M , and the squeezing operator S = exp(iHt) results from the time evolution under the Hamiltonian 23 .
Similarly, pumping the SQUID with a field that implements beam splitting and phase shifting operations between modes l, j. Note that arbitrary relative phase shifts on each mode are implemented either by a period of free evolution -no pumps-, since all the resonator modes possess a different frequency, or by phase-shifting the pump itself with an external phase φ lj . The coupling coefficient lj lj 1 0 is indeed the Bogoliubov coefficient α lj in equation (7). As shown by Reck et al. 22 , any unitary operator R U can be decomposed in these passive linear operations, 2 , and U k connects nearest neighbor modes. Interestingly, the circuit depth can be reduced to  = K M M ( log ( )) by implementing non-local beam-splitters 8 . This can be naturally implemented in our proposal, as we are dealing with beam splitters in frequency space and then we can connect any pair of modes i, j by choosing the right frequency for the pump |ω i − ω j |. This has the advantage of reducing the number of operations considerably.
The final step of the protocol, after evolving under red and blue sideband Hamiltonians, is reading out the number of photons on each mode of the resonators. In the diluted limit of boson sampling, where the number of modes M = N 2 , more than one excitation per mode is unlikely to occur, for which parity measurements would be enough. In this case, one can perform Ramsey-type measurements, where ancilla qubits are coupled to a different mode of the resonator. For each qubit j of frequency Ω j , one applies two π/2 pulses separated by a time π = Δ t g / n j n n j , 2 , , where Δ n,j = ω n − Ω j is the qubit-nth cavity mode detuning. With the qubits initially in the ground state, the Ramsey-type measurement maps the even (odd) parity onto the excited (ground) state of the qubit. The qubit state can be finally measured by a projective measurement, revealing the n-th mode state partity, and thus whether there is 0 or 1 photon in such mode 24 . Alternatively, one could also perform number-resolving measurments on each mode of the resonator by measuring a photon number-dependent energy splitting on the ancillary qubits, as described elsewhere 25 .
We would like to highlight that a similar system of coupled resonators with tunable coupling has already been implemented in the laboratory 13 . Using a SQUID as a coupler would not increase the experimental requirements 6,13 . Multimode parametric amplification by means of SQUID boundary condition modulation has already been reported as well 14 . Note that the experimental state -of-the-art in boson sampling with optical setups is still within the regime of low number of photons (three photons up to six modes) [26][27][28][29] , while integrated photonic circuits have achieved a maximum number of three photons in thirteen modes probabilistically generated from six SPDC sources 30 . Using this number as benchmark, our setup would require addressing 4 modes per resonator, which is completely within experimental reach in multimode circuit QED setups 31 . While there has no been experiments implementing boson sampling with superconducting circuits so far, at least another realistic proposal exists 25 . Our setup might be less resource-consuming, since it only involves two resonators instead of a large array.
In order to remain within the employed perturbative approximations, we should have |β jl | < 0.1 for all the relevant j,l. Using Eq. (10) with the realistic parameters  A 1 mm and  ω 10 GHz p , this implies a feasible time duration of the pulses of around 100 ns, for each pulse involved in both the state preparation and the unitary evolution. Since the average number of photons in a two-mode squeezed state is given by the square of the beta coefficient, this means that we would need around 100 repetitions in order to achieve successful single-photon detections. Putting everything together and considering a measurement time of around 1 μs 24 we can predict an event rate of approximately kHz in the low photon number regime  n 3, which could be improved with faster readout times, as in 32 .
It is important to remark that errors stemming from noisy state preparation, imperfect implementation of the unitary U and measurement, will occur. This is one of the main limitations in any practical implementation of boson sampling, as the error rapidly scales with the size of the simulator. The implementation proposed in this work is not error-free either and, while the noise analysis may yield promising results as other circuit QED implementations of boson sampling 25 , a careful study to figure out error thresholds for scalability is still needed.

Complexity of Dynamical Casimir Effect.
Let us discuss the computational complexity inherent to a randomized DCE-like evolution. So far we have seen how DCE implemented in a two-coupled superconducting resonator system acts a quantum simulator of GBS by virtue of simple red-and blue-sideband dynamics. However, one could think of a more general scenario, in which the SQUID is fed with a multimode magnetic field l j M l lj j l lj j , 1 which resembles a generalized boson sampling Hamiltonian 23 . It is not difficult to realize that the generalized Anger-Jacobi expansion yields a one-to-one relation between the external field amplitudes Φ ± lj and the coefficients g lj , ξ lj in terms of multivariate normal moments 11,[33][34][35] . Using this mapping, and provided that the magnetic field amplitudes Φ ± lj -more precisely the dimensionless ratio Φ Φ ± / lj 0 -are drawn from a random Haar measure, we conclude that a randomized DCE-like evolution lies outside the complexity class P, thus implementing a task that it is widely believed to be classically hard.
While cutting edge signal generators are capable of creating train pulses with hundreds of frequencies, and resonators allocating hundreds of modes have been developed 36 , a sufficiently large randomized DCE experiment that exhibit quantum supremacy 9,25 seems challenging due to the current limitation in resonator lifetimes and frequency-resolving measurements. A promising idea could be replace our standard transmission line resonators by left-handed transmission line metamaterials, where a very dense mode spectrum has already been reported [37][38][39][40] . However, we believe that this remarkable implication about the DCE computational complexity will trigger the forthcoming development of DCE-like experiments.

Discussion
In summary, we have shown how to use the DCE in order to implement a GBS device. We propose a setup consisting of two superconducting transmission line resonators with different energy spectra that are coupled by a SQUID. The ultrafast modulation of the magnetic field fed into the SQUID by an external pump plays the role of relativistic motion of a mirror shared by the two resonators. The corresponding Bogoliubov transformation results into multimode parametric amplification. We show how a suitable choice of parameters allows to implement GBS and in particular scattershot BS, thus demonstrating that DCE can be used to implement a task that it is widely believed to be classically hard. Moreover, we show that randomized DCE-like dynamics should be itself classically hard.