Extreme quantum nonlinearity in superfluid thin-film surface waves

We show that highly confined superfluid films are extremely nonlinear mechanical resonators, offering the prospect to realize a mechanical qubit. Specifically, we consider third-sound surface waves, with nonlinearities introduced by the van der Waals interaction with the substrate. Confining these waves to a disk, we derive analytic expressions for the cubic and quartic nonlinearities and determine the resonance frequency shifts they introduce. We predict single-phonon shifts that are three orders of magnitude larger than in current state-of-the-art nonlinear resonators. Combined with the exquisitely low intrinsic dissipation of superfluid helium and the strongly suppressed acoustic radiation loss in phononic crystal cavities, we predict that this could allow blockade interactions between phonons as well as two-level-system-like behavior. Our work provides a pathway towards extreme mechanical nonlinearities, and towards quantum devices that use mechanical resonators as qubits.


INTRODUCTION
Nonlinearities are widely used in quantum technologies. For instance, they allow the generation of nonclassical states 1-6 , twoqubit interactions [7][8][9] , and quantum nondemolition measurements [10][11][12][13][14] . Sufficiently strong nonlinearities can introduce resolvable anharmonicity in a resonator, so that when resonantly driven it can only absorb a single quantum of energy, mimicking the behavior of a two-level system. This provides the possibility of blockade-type interactions, where phonons (or photons, depending on the resonator) can only pass through the resonator one at a time [15][16][17] . It also allows artificial qubits to be engineered, such as the superconducting qubits widely used in quantum computing 18 .
Mechanical nonlinear resonators have quantum applications ranging from the preparation of nonclassical states [19][20][21][22] to quantum-enhanced force sensing [23][24][25][26][27] , quantum backactionevading measurement 28 , and mechanical quantum state tomography 29 . Achieving the single-phonon nonlinear regime in a mechanical resonator is of both fundamental and technological importance. It would allow artificial atoms to be built from massive objects consisting of billions of atoms, testing quantum physics in uncharted regimes of macroscopicity, and would provide a qubit for quantum computation among other quantum applications 17,30 .
Reaching the single-phonon nonlinear regime in a mechanical resonator requires an intrinsic nonlinearity far stronger than what has been achieved to date [31][32][33][34][35][36][37] , combined with exceptionally low dissipation so that the energy level shifts introduced by the nonlinearity are resolvable. Here, we propose to achieve this using a thin spatially confined superfluid helium film, similar to the ones used in recent experimental work on optomechanical cooling 38 , lasing 39 , and quantized vortex detection 40 . The superfluid resonator is a third-sound surface wave with restoring force provided by the Van der Waals interaction with the substrate.
We derive an analytical model of the cubic and quartic (Duffing) nonlinearities due to van der Waals forces for a film confined on a circular disk. We find that the nonlinearities depend strongly on the radius of confinement, and predict that the quartic nonlinearity in a 5-nm-thick film with 100 nm radius would manifest single-phonon frequency shifts three orders of magnitude larger than those seen in state-of-the-art nonlinear mechanical resonators including graphene sheets 36 , carbon nanotubes 36 , and molecule-coupled resonators 37 . Our analysis shows that the primary effect of the cubic nonlinearity is to modify the magnitude of the quartic term, consistent with previous work on classical mechanical resonators [41][42][43] . This result has broad relevance, since a wide variety of mechanical resonators exhibit both cubic and quartic nonlinearities 41 .
Even with extremely high nonlinearities, achieving sufficiently low dissipation to enter the single-phonon nonlinear regime is a significant challenge. Superfluid helium affords exceptionally low intrinsic dissipation. Indeed, sub-millihertz dissipation rates have been observed in third-sound resonators with millimeter dimensions 44 . However, radiative dissipation increases with increasing spatial confinement. We introduce the concept of superfluid thinfilm phononic crystal cavities to overcome this, showing that radiative dissipation can be greatly suppressed, even at hundrednanometer size scales. We further predict that dissipation due to thermalization and vortices should not be a barrier to reaching the single-phonon nonlinear regime. Together with the level of nonlinearity predicted by our model, this suggests that the single-phonon nonlinear regime can be reached, opening a path to probe quantum macroscopicity in a new domain and to build a new class of qubits for quantum computing and metrology.

RESULTS
The anharmonic superfluid oscillator potential Liquid helium exhibits superfluidity below a critical temperature T λ . It can then be described as an effective mixture of a normal fluid with density ρ n and a superfluid with density ρ, with total density ρ He = ρ + ρ n 45,46 . At temperatures well below the critical temperature the ratio ρ/ρ He approaches 1. For instance, for a critical temperature T λ = 2.2 K, ρ/ρ He > 0.98 for temperatures beneath 1 K 47 . In this low-temperature limit, superfluid helium has a combination of traits often sought after in mechanical resonators: low mechanical dissipation arising from near-zero viscosity and ultralow optical absorption. Indeed, it has been used as the mechanical resonator in several recent optomechanical platforms 38,39,[48][49][50][51][52][53][54][55] and in experiments that study the physics of quantum fluids 40,56 . In these references, the superfluid fills a cavity [51][52][53] or channel 54 , is levitated as a droplet 50 , or spread out on a surface as a few nanometers thin film 38,39,44,57 . The latter case is investigated here. Due to the thinness of the film, the normal component can be considered viscously clamped to the surface 45 , while the superfluid component exhibits thickness fluctuations that resemble shallow water waves, as illustrated in Fig. 1a. These waves are named "third sound" and are unique to twodimensional superfluid helium films 44,57,58 .
In this work, we consider the superfluid film to be confined to a circular surface of radius R. This geometry is quite general, and can be realized for instance by condensing the film on the surface of a microscopic silica disk (see Fig. 1b) 38,39,49,59 . That design is attractive because it allows laser light to circulate in the disk. These "whispering-gallery" light waves interact strongly with the third-sound waves in the superfluid, and can serve as a tool to observe and control the superfluid motion 38 . In this study, however, we focus on the film's dynamics: while constrained here to a circular disk, we expect our predictions to be qualitatively mirrored in other superfluid thin-film geometries.
A helium atom at height z is attracted to the substrate atoms via the van der Waals force 45,46 . This leads to a height-dependent potential energy per unit mass stored in a film. The analysis here is restricted to films with mean thicknesses d between 1 and 30 nm, and with diameters much larger than their thickness. In this case the potential is well approximated by 60,61 with a vdw the substrate-dependent van der Waals coefficient characterizing the attraction strength. The scaling with height is modified for thicker and lower-aspect-ratio films 46,62 , while for thinner films (on the order of a few atomic layers) corrections due to the approximately one inactive atomic layer must be taken into account [63][64][65] . The potential provides a restoring force for fluctuations of the film surface. Turning to Fig. 1a, the circularly confined film somewhat resembles a drumhead-and in fact, the helium surface undulates like the skin of a resonating drum. While it is clear that Eq. (1) is nonlinear and therefore does not describe a Hookean potential, in the small amplitude limit (where nonlinearities can be neglected) the resonances of the surface can be described by Bessel modes 59 . These eigenmodes are, strictly speaking, valid only for a linear oscillator. However, they are a good approximation for the high-quality (see "Results") mechanical resonances considered here where the nonlinearity only shifts the mechanical resonance frequency by a small fraction. The time-dependent Bessel mode amplitude h r; θ; t ½ that quantifies the deviation of the film height from the mean thickness d with polar coordinates r and θ is given by 59 Here, A is the amplitude, J μ is the Bessel function of the first kind of integer order μ, and Ω m = (ζ μ,ν c 3 )/R is the mechanical resonance frequency in the absence of nonlinearity. Here the resonance frequency depends on the superfluid speed of sound c 3 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 3a vdw ðρ=ρ He Þ d À3 q and a parameter ζ μ,ν , which depends on the boundary conditions. In the absence of flow across the resonator boundary, the film is described by volume-conserving Bessel modes, i.e., Bessel functions with free boundary conditions 59 . These mode amplitudes have an extremum at r = R: a condition met by choosing ζ μ,ν as the νth zero of J 0 μ . We restrict our analysis to these boundary conditions since they match previous observations with disk resonators in ref. 40 . However, fixed boundary conditions could alternatively be used by choosing the zeroes of J μ for the coefficient ζ μ,ν . η[r, θ] is the time-independent amplitude of a drumhead mode typically specified by its mode numbers (μ, ν) with the order μ the number of nodal diameters, also called azimuthal mode number (i.e., 2μ is the number of zeros in the azimuthal direction for θ = 0-2π) and ν the number of nodal circles, also called radial mode number (i.e., ν is the number of zeros in the radial direction for r = 0 to R). The Bessel mode function J μ ζ μ;ν r R Â Ã is graphed in Fig. 2, and the values ζ μ,ν for the first three mode numbers are tabulated in Table 1.
The van der Waals potential energy stored in the surface deformation U is obtained by integrating the potential V over the deviation from equilibrium 59 : with Eq. (1) and the Taylor series expansion, we obtain The majority of previous work on third sound considers only the first two terms in this expansion, which describe the potential of an harmonic oscillator. The higher-order terms introduce nonlinearities. Experiments have investigated a variety of nonlinearities 58,64,[66][67][68][69][70] in millimeter and centimeter-scale devices, far away from the regime we consider here. The first term (proportional to η) averages out to zero for volume-conserving Bessel modes with R and fourth-order terms, respectively, represent cubic and quartic nonlinearities. We neglect terms of fifth order in η and higher, expecting them to be small compared to these first two nonlinear terms. Thus, Eq. (4) becomes where we have identified the quadratic, cubic, and quartic potential energies associated with linear, quadratic, and cubic spring terms. (In this work we label the nonlinearities by the orders in which they appear as energies and potentials, because we work exclusively in the Hamiltonian formalism. In some other work, especially early work on spring forces, they may be labeled by their lower order in force and acceleration equations.) By introducing a reference point x: = η[R, 0]-that is, the displacement at the periphery of the disk at an angular location θ = 0-we can rewrite the potential energy as This is the potential energy of an anharmonic oscillator (see Fig.  1c and β and α are the nonlinear spring constants. The cubic nonlinearity is given by (10) and the quartic (also known as Duffing) nonlinearity is given by By evaluating the integrals (see Supplementary Methods), the strong dependence of the (non)linear spring constants on the film Ã of (azimuthal) orders μ = 0, 1, and 2 for radial mode number ν = 1-3. thickness d and confinement radius R is exposed: and Here we have introduced the Kronecker delta function δ, and integrals of the Bessel function ϕ ðpÞ μ;ν :¼ Table 2. It is worth noting that the cubic nonlinearity β vanishes for all but the rotationally invariant (μ = 0) modes. The film thickness d can be independently determined and tuned in situ, as we demonstrated in previous work 38 , while the confinement radius R can be changed through choice of device geometry. By controlling these, Eqs. (14) and (15) reveal that it is possible to access a wide range of cubic and quartic nonlinearities.
Although the relationships found here suggest that the nonlinear coefficients are larger for thinner films and larger radii, the desired parameter regime depends on the relative magnitudes of the nonlinear and linear coefficients, and the mass. In the following section, we will derive what exactly that regime is and what platform parameters one should aim for in order to reach it.

Single-phonon transition resonances
We consider the Hamiltonian for the oscillator with natural frequency Ω m , effective mass m eff ¼ k=Ω 2 m and zero-point with a (a † ) the phonon annihilation (creation) operators satisfying the commutation relation [a, a † ] = 1, the canonical position and momentum operators x and p satisfying x = x zpf (a + a † ) and [x, p] = i_, and n = a † a the phonon number. The nonlinearities modify the eigenenergies of the oscillator (see Fig. 3). To second order in perturbation theory we find that the transition between states n and n + 1 occurs at frequency (see Supplementary Methods) with the single-phonon nonlinear shift and α eff [α, β] the effective quartic nonlinearity absorbing the cubic nonlinearity: The frequency shift in Eq. (18) is the mechanical analog of the Kerr nonlinearity in quantum electrodynamics 71 , but with the quartic nonlinear coefficient α reduced to α eff due to the presence of the cubic nonlinearity. The modification, Eq. (19), agrees with the classical result found in refs. [41][42][43] .
When the mechanical decay rate (or linewidth) Γ of the oscillator is smaller than the spectral splitting, i.e., the resonator is sufficiently anharmonic that a single absorbed phonon shifts the frequency off resonance for phonons that arrive later, causing phonon blockade: it behaves more like a two-level system than an harmonic oscillator 1,2 .
The criterion in Eq. (20) can be re-expressed in terms of the , which defines the amplitude at which a classical Duffing resonator become bistable. In this form the criterion is x zpf > x crit , that is, the mechanical zero-point motion must exceed the critical amplitude to reach the singlephonon nonlinear regime 37 .
Open quantum system spectrum To determine the expected resonator spectrum in the presence of decoherence and validate our perturbation theory calculation, we numerically solve the Lindblad master equation for a resonator with cubic and quartic nonlinearities (see Supplementary Methods). We compare the resulting correlation spectra S xx for three cases: a superfluid resonator with finite quartic and cubic nonlinearity S xx [α, β], a pure Duffing resonator It is evident at a glance that as the single-phonon nonlinear strength δΩ/Γ exceeds unity, the nonlinearity lifts the oscillator's spectral degeneracy. The spectra for pure Duffing resonators S xx [α, 0] manifest transition resonances at the values Ω[n] = Ω m + (n + 1)δΩ[α] obtained analytically from first-order perturbation theory. (The small discrepancy stems from the numerical error on the bare energy eigenstates E n , which increases with n. This absolute error is inversely proportional to the size of the basis spanning the Hilbert space in the numerical algorithm; the smaller α and δΩ, the larger the basis must be chosen to mitigate the relative error on the transition frequencies.) The pure Duffing spectrum does, however, differ significantly from the full-Hamiltonian third-sound resonator spectrum S xx [α, β], and from that of the effective Duffing third-sound resonator S xx [α eff , 0]. The latter two are nearly identical in all cases, both in terms of amplitude and transition resonance frequencies. This indicates that the third-sound resonator with a quadratic nonlinearity β is accurately described by the effective Duffing resonator according to Eq. (21), and its single-phonon transition frequencies are correctly analytically predicted as Ω[n] = Ω m + (n + 1)δΩ[α eff ] from Eq. (18).
Single-phonon nonlinear shift for a superfluid resonator Using Eqs. (12) to (15) and (19) we can express δΩ in terms of the adjustable parameters R and d. With one obtains This equation is valid in the limit of a third-sound wavelength large compared to the film thickness d, and a motional amplitude η ≪ d, which is the case for all examples considered here. It identifies every parameter available to the researcher who seeks to maximize the single-phonon nonlinear strength of the superfluid resonator: the shift grows strongly with decreasing confinement radius R, scales with the inverse of the film thickness d, but is independent of the van der Waals coefficient a vdw between the helium film and the substrate. The shift is always positive, so the oscillator is effectively "spring-hardened". Its dependence on the mode numbers (μ; ν) is rather intricate due to the lack of a closed form of the coefficients ζ μ,ν and ϕ ð4Þ μ;ν (see Supplementary Methods). Therefore, we have graphed the frequency shift as a function of the mode numbers (μ; ν) in Fig. 5. In this figure and henceforth we take ρ/ρ He → 1, considering only the regime where the temperature is well beneath T λ . It can be seen that the nonlinear strength δΩ can grow by four orders of magnitude when the Bessel mode order μ and radial mode number ν vary from 0 (for the mode order) and 1 (for the radial mode number) to 20. The predicted nonlinearity of the superfluid film is compared to other systems in Table 3 for a range of confinement radii R and radial nodes ν. From the table it is clear that the simultaneous attainment of large α and x zpf is nontrivial for many platforms, while for the superfluid third-sound resonator both are intrinsically large and tunable. When a 5-nm-thick superfluid film is confined to a radius of 10 μm, its predicted single-phonon frequency shift surpasses those in membranes 31 , cantilevers 32 , and Si 3 N 4 beams [33][34][35] . For R = 1 μm, it exceeds levitated nanoparticles that inherit their strong nonlinearity from an optical trapping potential 72 and resonators with an engineered chemical bond 37 . Finally, in the submicron regime, single-phonon nonlinear shifts might surpass by orders of magnitude those of carbon nanotubes and graphene sheets 36 , exceeding the single-phonon nonlinear threshold.
Reaching the single-phonon regime in thin superfluid helium Having obtained the parameter space required to bring a thin-film superfluid resonator into the single-phonon nonlinear regime, the question becomes: Are there platforms available that may facilitate these requirements? Can one reasonably engineer an on-chip superfluid resonator whose damping Γ approaches the single-phonon frequency shift δΩ?
In the earliest third-sound resonators 70,73 , the superfluid film was adsorbed on the inside of two parallel metalized silica disks. This approach allowed capacitive detection of the film's dynamics: film thickness variations change the capacitance between the plates. In previous works, we have used a few-nanometer-thick film adsorbed on a single on-chip silica microdisk allowing a reduction in confinement radius of three orders of magnitude together with optical detection of the film's dynamics. [38][39][40]49 . Such disks, however, become unsuitable for confinement below 10 μm. Using high refractive index materials such as silicon or gallium arsenide instead, the disks can be miniaturized down tõ 1 μm radii 74 . These set-ups span the hatched region in Fig. 6a; it can be seen that reaching the single-phonon nonlinear regime then requires damping rates in the millihertz regime. This is a challenging condition, yet such low damping rates-with corresponding mechanical quality factors Q = Ω m /Γ in excess of 10 5 -have been experimentally observed in ref. 44 , albeit for larger millimeter dimensions. For bulk modes in superfluid 4 He, quality factors in excess of 10 8 have been measured 75 .
Further size reduction could be achieved through the use of phononic crystal cavities [76][77][78] . When condensed on a patterned silicon substrate, the superfluid third-sound wave experiences a periodic modulation of its speed of sound and a periodic potential 79 . This provides the possibility to engineer phononic crystal band structures that contain frequency ranges, or band gaps, for which wave propagation in the crystal is not allowed. In analogy to photonic crystals and solid-state phononic crystals, we find that third-sound modes can be trapped or confined in a defect in the crystal as long as its frequency lies in the band gap, as shown in Fig. 6b. Thereby, phononic crystal cavities could enable confinement down to and below hundreds of nanometers (Fig. 6a, gray).
Here, as an example, we calculate the band structure for a 11nm-thick superfluid film condensed on a suspended silicon slab perforated with 55 nm diameter holes and a 100-nm lattice constant, as shown in Fig. 6b. The hydrodynamic third-sound equations are solved through finite-element method (FEM) simulation (see Supplementary Methods and ref. 56 ), showing a complete band gap at around 30 MHz. When a hole is removed to form a central defect, it confines a 30 MHz acoustic mode. This trapped mode closely resembles the (μ = 0, ν = 1) R = 56 nm circular Bessel mode (see Supplementary Fig. 2). From Eq. (24), the corresponding single-phonon frequency shift is δΩ[α eff ]/2π = 35 Hz (indicated with an X-mark in Fig. 6a), so that a quality factor of 9 × 10 5 is required to reach the single-phonon nonlinear regime.
Unambiguous identification of all damping mechanisms in third sound remains an open problem in the field 70,80-83 . The known main dissipation channels are thermal dissipation that arises due to the temperature gradient between the peaks and troughs of the sound wave, dissipation due to interactions with pinnedvortices 70,84 , and radiation [38][39][40] (or clamping 85 ) losses.
Aside from allowing submicron confinement, a virtue of the phononic crystal architecture is the ability to control the radiation loss associated with imperfect reflection of the wave at the confining boundary. Our model confirms (see Supplementary Fig.  3) that the radiation damping decreases strongly as the number of cells in the phononic crystal lattice increases. Indeed, for the 30 MHz mode in Fig. 6 the threshold for single-phonon resolution Q = 9 × 10 5 is reached for just 10-cell lattices.
Thermal dissipation occurs as the result both of evaporation and recondensation within the superfluid wave itself 45 , and of irreversible heat flow through the film and substrate 45,82,86 (similar to thermoelastic losses in micromechanical resonators 87 ). These mechanisms were modeled for a plane wave by Atkins in ref. 45 . Using this model (see Supplementary Fig. 4) we find that these mechanisms are greatly suppressed at low temperatures (where   the temperature gradients in the wave are reduced) and also for high mechanical frequencies. The suppression with increasing mechanical frequency is perhaps surprising. However, it can be understood since if the frequency is significantly higher than the rate of thermal equilibration, then heat has insufficient time to flow between troughs and peaks. A similar phenomenon has been observed for thermoelastic losses 88 . Combined, we predict that thermal losses will not be a significant source of dissipation for strongly confined films at low temperatures. For the specific example of the phononic crystal structure modeled in Fig. 6b, we find that the thermal dissipation rate is suppressed beneath the frequency shift due to a single-phonon even at temperatures as high as 0.5 K. Indeed, sub-hundred millikelvin thermalization is not uncommon for third sound 44,58,63,64,66,68,70,89,90 , in which case Atkins' model predicts that the thermal dissipation dominated superfluid quality factor Q thermal would exceed 10 10 . Dissipation resulting from interactions with quantized vortices pinned to the resonator surface is thought to arise due to both vortex-normal fluid interactions and vortex dimple drag 70,84 . These damping mechanisms require large pinned vortex densities in order to account for the observed dissipation, with densities on the order of 10 13 cm −2 estimated in ref. 84 . Recent experimental work, however, shows that coherent vortex-vortex interactions dominate pinning when the superfluid film is confined at hundred Fig. 6 Radial dependence of single-phonon shift and band structure phononic crystal. a Resonance frequency shift induced by a single phonon δΩ ¼ 3α eff x 4 zpf =_ as a function of mode confinement radius R for various superfluid film thicknesses d for the (μ = 0; ν = 1) mode. The quantized nature of the anharmonic oscillator energy spectrum is resolvable if the mechanical linewidth Γ < δΩ. Hatched: confinement radii for microdisks [38][39][40]49 . Gray shaded: confinement radii achievable through phononic crystal (PnC) trapping. X-mark: mode confined in PnC cavity with lattice constant 100 nm. b PnC band structure and complete band gap (gray shaded) for a hexagonal honeycomb lattice with holes 55 nm in diameter, lattice constant 100 nm, and the van der Waals coefficient for silicon a vdw = 3.5 × 10 −24 m 5 s −2 (ref. 60 ). A central defect (absence of hole) introduces a Bessel-like mode at 30 MHz-within the band gap-as shown by the finite-element-method simulation. micron-scales on smooth microfabricated resonators 40 . The vortex distribution then evolves towards its lowest energy, vortex-free, configuration over minute periods. As the resonator radius is scaled down, it is expected that the influence of pinning sites will be further reduced since the sound-vortex coupling scales inversely with the resonator area 56 , causing vortices to be dislodged from pinning sites and their subsequent annihilation.
Liquid helium only remains superfluid as long as the fluid particle velocity is less than the superfluid critical velocity. It is interesting then to ask whether there are any fundamental limits to how strongly the superfluid film may be confined, both in thickness and radius, before superfluidity breaks down, and whether these constrain the possibility of entering the single phonon nonlinear regime. In the Supplementary Discussion we consider two cases, first a third-sound mode cooled to its motional ground state, and second a third-sound mode in thermal equilibrium with its environment. For both cases, we find that a 10-nm-thick film would need to be confined to a radius beneath its thickness for the fluid particle velocity to exceed the critical velocity (a parameter regime outside the validity of our model). We conclude, therefore, that the breakdown of superfluidity places no constraints on reaching the single phonon strong coupling regime with the hundred-nanometer-and-above diameter superfluid resonators considered here.

DISCUSSION
We have shown that third-sound resonances (surface oscillations of two-dimensional superfluid helium) are intrinsically strongly nonlinear, and have identified the specific parameters that allow maximization of the nonlinearities: film thickness, confinement radius, and radial and rotational mode symmetry. We showed that the cubic nonlinearity can be treated analytically as an effective reduction of the quartic nonlinearity, which diminishes rapidly for surface waves with an increasing number of radial nodes.
We calculated the expected output spectrum in the presence of decoherence and predict that single-phonon nonlinear frequency shifts exceeding even those of graphene sheets and carbon nanotubes by orders of magnitude may be possible. Combined with the intrinsically low dissipation of motional states in superfluid 44 , and phononic crystal cavities providing sub-100-nm confinement and strongly suppressed radiation damping, this may well open the door to the single-phonon nonlinear regime where a single phonon can shift the resonance frequency by more than the mechanical linewidth.
Our results dovetail recent theoretical proposals that lay out how strong Duffing nonlinearities can be used for quantum control and metrology 19,24,25,27,91,92 . Ultimately, they provide a pathway towards tests of quantum macroscopicity and new tools for quantum technologies, where mechanical resonators can function not only as oscillators, memories, and interfaces, but also as qubits.

DATA AVAILABILITY
The computational code and data sets generated and analyzed during the current study are available from the corresponding author upon reasonable request. Numerical simulations were done using the Quantum Optics Toolbox for MATLAB 93 , available at https://qo.phy.auckland.ac.nz/toolbox/.