Method and unit for characterizing a quantum system

Abstract

The invention provides a method and a quantum characterisation unit that (i) prepares the quantum system in at least two (preferably three or more) distinct dynamical configurations, each describable by its own stochastic master equation (SME) but sharing a subset of physical parameters; (ii) records continuous detector signals, bins or filters them, and computes low order correlation functions (one, two, three point, optionally higher order) from repeated measurements; and (iii) jointly fits simulated correlation functions — generated directly from the SME without trajectory sampling — to the experimentally obtained ones, using a least squares (or similar) optimisation to assign values to the shared parameters. The method exploits the fact that different configurations break parameter degeneracies that would be hidden in a single configuration.

Description

[0001]METHOD AND UNIT FOR CHARACTERIZING A QUANTUM SYSTEM FIELD OF THE INVENTION The invention concerns a method of characterizing a quantum system measured continuously, and in particular characterizing a quantum system for stabilizing one or more bosonic qubits, as well as a quantum characterization unit for doing same. BACKGROUND Characterising the dynamics of a quantum system is a fundamental task in experimental quantum physics. Typically, a series of hand-crafted measurements are performed in a precise order to estimate the value of each system parameter in turn. However, as quantum technologies mature, this approach is no longer efficient to characterize increasingly large and complex systems, such as recent demonstrations of quantum error correction embedding multiple superconducting qubits in a surface code or stabilisation of bosonic qubits using continuous variable encodings which leverage relatively large Hilbert space dimensions. In particular, continuous variable encodings, involving superpositions of specific harmonic oscillator states with infinite-dimensional Hilbert spaces (e.g., bosonic qubits such as GKP qubits, Fock states, cat qubits), show promise for protecting and processing quantum information efficiently with lower overhead resources compared to discrete variable systems with finite-dimensional spaces. Specifically stabilized cat qubits (which can be hosted in superconducting resonators) benefit from a noise bias, where bit-flip errors are exponentially suppressed with the average number of bosons (normally photons) in the so-called Schrödinger cat states hosted in a resonator. This suppression is effective for a wide range of physical noise processes, including photon loss, thermal excitations, photon dephasing, and nonlinearities from Josephson junctions. Recent experiments in quantum superconducting circuits have demonstrated this exponential suppression of bit-flip errors. However, such experiments have also demonstrated difficulties in accurately characterizing the physical parameters of such highly dynamically non-linear quantum systems. For instance, Berdou, Camille, et al. "One hundred second bit-flip time in a two-photon dissipative oscillator." PRX Quantum 4.2 (2023): 020350 discloses a method of characterizing the photon number ^^2in a two-photon dissipative oscillator, the associated two-photon dissipation rate ^^2and a term ^^2arising from and proportional to the amplitude of an external drive. However, a specifically designed sequence of separate experiments, including measuring the variance of the measured signal Ε[^^^2^] for various values of ^^2(where ^^^^is the discrete-time digitised signals resulting from detector signals of the continuous- time analogue signals from the oscillator) was required to characterise these three parameters, with a resulting uncertainty on each estimate being larger than 10%. Bayesian inference is a method for estimating parameters from continuous measurement data, such as disclosed Mabuchi, Hideo. "Dynamical identification of open quantum systems." Quantum and Semiclassical Optics: Journal of the European Optical Society Part B 8.6 (1996): 1103. These make optimal use of the available information by selecting parameters that maximise the likelihood of the observed data, e.g. by computing the trace of the solution to the linear stochastic master equation (SME), which is solved using the measured signal. However, the numerical cost is prohibitive for estimating many parameters, and they do not take into account signal filtering and digitisation. Tilloy, Antoine. "Exact signal correlators in continuous quantum measurements." Physical Review A 98.1 (2018): 010104 discloses an exact formula for calculating correlation functions of detectors continuously measuring an arbitrary quantum system, in the presence of detection imperfections, for diffusive measurement. Guilmin, Pierre, Pierre Rouchon, and Antoine Tilloy. "Correlation functions for realistic continuous quantum measurement." IFAC-PapersOnLine 56.2 (2023): 5164-5170 discloses similar exact formula for calculating correlation functions, but generalized to quantum jump measurements. However, both disclosures only present theoretical derivations of formulae for computing correlation functions directly from the theoretical dynamical model (the stochastic master equation) and the theoretical initial state of an arbitrary quantum system. Thus, a major challenge to streamline the development of any of the above technologies is to find characterization methods that (i) facilitate the automated estimation of many parameters, and (ii) remain practical for everyday use in the laboratory. Moreover, successful operation of a quantum processor requires regular re- calibration routines to bring the delicate quantum systems back into operability regimes. Such routines require rapid and accurate characterization of the quantum systems in order to apply effective correcting calibration operations. This invention aims to address one or more of the aforementioned problems to provide an improved method of characterizing a quantum system. SUMMARY In a first aspect, there is provided a method of characterizing a quantum system, wherein the quantum system comprises one or more qubit-hosting structure(s), each of the one or more qubit-hosting structure(s) comprising: (I) at least one portion configured to have at least one physical oscillatory mode, (II) one or more signal generator(s) coupled to the at least one portion and configured to input control signals to physically prepare the quantum system in a dynamical configuration, and (III) a detector coupled to the at least one resonant portion and configured to measure continuous signals output from the at least one oscillatory mode to provide measurement signals from the quantum system. The method comprises: (i) physically preparing, with the one or more signal generator(s), the quantum system in a first dynamical configuration, wherein the first dynamical configuration is describable via a first stochastic master equation having a first plurality of physical parameters; (ii) physically measuring, with the detector, a continuous signal from the quantum system in the first dynamical configuration to provide a first measured continuous signal; (iii) repeating steps (i) and (ii) a plurality of times to provide a plurality of first measured continuous signals, and calculating one or more first correlation function(s) from the plurality of first measured continuous signals; (iv) physically preparing, with the one or more signal generator(s), the quantum system in a second dynamical configuration which is different to the first physical configuration, wherein the second dynamical configuration is describable via a second stochastic master equation having a second plurality of physical parameters, and wherein the quantum system is physically prepared in the second dynamical configuration such that a subset of the second plurality of physical parameters is a shared subset with a subset of the first plurality of physical parameters; (v) physically measuring, with the detector, a continuous signal from the quantum system in the second dynamical configuration to provide a second measured continuous signal; (vi) repeating steps (iv) and (v) a plurality of times to provide a plurality of second measured continuous signals, and calculating one or more second correlation function(s) from the plurality of second measured continuous signals; (vii) assigning values to the shared subset of physical parameters which characterizes the quantum system by jointly numerically fitting: (a) a first simulation of the one or more first correlation function(s), numerically simulated from the first stochastic master equation, to the calculated one or more first correlation function(s), and (b) a second simulation of the one or more second correlation function(s), numerically simulated from the second stochastic master equation, to the one or more second correlation function(s). Thus, the method according to the present invention measures continuous signals of a quantum system in different dynamical configurations, and subsequently fits the statistics of the measured signals (in the form of correlation functions, which leverages temporal information from the correlation functions at non-coincident time) to efficient numerical simulations so as to extract the values of the physical parameters which characterize the quantum system. In contrast to prior art methods which aim to make optimal use if the available information from observed data (e.g. via computationally expensive maximum likelihood methods), the present method instead directly fits statistics of the observed experimental data with numerically simulated statistics in the form of correlation functions from a continuous signal, and without necessarily trying to reconstruct intermediate quantities such as the quantum state. Although the present method does not make optimal use of the available information in a mathematical sense, the present Inventors have discovered that the experimental processing (namely the physical measurements) is particularly straightforward and the fitting procedure is both easy to implement numerically, and extremely efficient in that it allows for the simultaneous estimation of multiple parameters and is suitable for large Hilbert space dimensions. Moreover, the present method physically prepares and measures the quantum system in different dynamical configurations so as to extract physical parameters which characterize the quantum system and is shared between the different dynamical configurations. The present Inventors have discovered that, when some physical parameters are not observable through the continuous signal of a given dynamical configuration (e.g. have no effect on the measured signals average), then they can be hidden by some symmetry of the quantum system. These occur when the physical parameter is fundamentally unobservable, or when it has too small an effect on the continuous signal to be estimated with good precision. However, the present Inventors have recognized that for typical quantum systems, any generic physical parameter will influence in different ways the continuous signal between different dynamical configurations of the quantum system. By physically preparing and measuring the quantum system in different dynamical configurations, wherein the difference between second dynamical configuration and the first dynamical configuration is known (e.g. the second dynamical configuration may be another known “initial” state, or a configuration obtained from the first dynamical configuration by modifying a signal with a known proportionality constant), the present method crucially ensures that physical parameters are extracted which otherwise would be hidden. Thus, the synergistic combination of both (i) measuring continuous signals of a quantum system and numerically fitting correlation functions thereto in an efficient manner, and (ii) doing so in different dynamical configurations of the quantum system to ensure all physical parameters are observable, presents a new tool for characterizing quantum systems. This method may be particularly suitable for automating the characterization phase of a quantum system due to the increased speed of the characterization method relative to prior art processes, and because it allows the operator of the quantum system to characterize many parameters in a single physical experiment. Furthermore, the present method enables an inherently interpretable characterization process, in that an operator can build intuition on how each physical parameter affects the average trajectories of the continuous signal of the quantum system. For the avoidance of doubt, the first simulation is a simulation of the one or more first correlation function(s) for the first dynamical configuration and the second simulation is a simulation of the one or more second correlation function(s) for the second dynamical configuration. Assigning values to the shared subset of physical parameters may thus comprise: (viii) numerically simulating, from the first stochastic master equation, a plurality of the one or more first correlation function(s), wherein at least one physical parameter of the first plurality of physical parameters has a different value when numerically simulating each of the simulated one or more first correlation function(s); (iv) numerically simulating, from the second stochastic master equation, a plurality of simulations of the one or more second correlation function(s), wherein at least one physical parameter of the second plurality of physical parameters has a different value when numerically simulating each of the simulations of the one or more second correlation function(s); and (x) assigning values to the shared subset of physical parameters which characterizes the quantum system by fitting the plurality of simulations of the one or more first correlation functions of step (viii) to the calculated one or more first correlation functions of step (iii) and fitting the plurality of simulations of the one or more second correlation functions of step (ix) to the calculated one or more first correlation functions of step (vi). Said otherwise, numerically fitting a simulated correlation function with the corresponding measured correlation function may be achieved by (a) simulating the correlation function for an initial choice of physical parameters, (b) comparing this simulated correlation function with the corresponding measured correlation function and quantifying how close it fits, (c) simulating the correlation function for a different choice of physical parameters, (d) comparing this simulated correlation function with the corresponding measured correlation function and quantifying how close it fits, and (e) iterating until the quantification of how close one of the simulated correlation functions is to the corresponding measured correlation function is below a certain fit threshold. The choice of physical parameters for the final fitted simulated correlation functions is then taken to be the physical parameters of the measured correlation function. Thus, values for the shared subset of physical parameters is extracted and assigned by (A) fitting the plurality of simulations of the one or more first correlation functions to the calculated one or more first correlation functions from the plurality of first measured continuous signals and (B) fitting the plurality of simulations of the one or more second correlation functions to the calculated one or more second correlation functions from the plurality of second measured continuous signals. Numerically fitting the first and second simulations may comprise performing a least-squares fitting procedure. Preferably, numerically simulating, from the first (respectively second) stochastic master equation, a plurality of the one or more first (respectively second) correlation function(s) comprises numerically simulating the plurality of the one or more first (respectively second) correlation function(s) directly from the first (respectively second) stochastic master equation without sampling a plurality of trajectories. Here, sampling a plurality of trajectories means explicitly calculating the correlation function(s) by averaging over all realizations. Preferably, each of the plurality of first (respectively second) measured continuous signals is a filtered or binned signal, and numerically simulating, from the first (respectively second) stochastic master equation, a plurality of the one or more first (respectively second) correlation function(s) comprises numerically simulating the plurality of the one or more first (respectively second) correlation function(s) of filtered or binned signals. The one or more signal generator(s) may be directly coupled to the at least one portion, or may be indirectly coupled to the at least one portion (for instance by being coupled to another element or structure of the qubit-hosting structure such as a non- linear element, wherein said other element or structure is itself directly coupled to the at least one portion). The detector may be directly coupled to the at least one resonant portion, or may be indirectly coupled to the at least one portion (for instance by being coupled to another element or structure of the qubit-hosting structure such as a non-linear element, wherein said other element or structure is itself directly coupled to the at least one portion). By continuous signal herein it is meant time-continuous signal. The method may comprise: physically preparing, with the one or more signal generator(s), the quantum system in a third or further different dynamical configurations, wherein each of the third or further different dynamical configurations is describable via a third or further stochastic master equation having a third or further plurality of physical parameters, and wherein the quantum system is physically prepared in the third or further different dynamical configuration such that a subset of the third or further second plurality of physical parameters is part of the shared subset of physical parameters; physically measuring, with the detector, a continuous signal from the quantum system in the third or further different dynamical configurations to provide a third or further measured continuous signal; repeating said steps of physically preparing and physically measuring the third or further different dynamical configurations a plurality of times to provide a plurality of the third or further measured continuous signals, and calculating one or more third or further correlation function(s) from the plurality of the third or further measured continuous signals; and wherein assigning values to the shared subset of physical parameters which characterizes the quantum system further comprises numerically fitting a third or further simulation, numerically simulated from the third or further stochastic master equation, to the one or more third or further correlation function(s). The present inventors have discovered that characterizing the quantum system in at least 3 different dynamical configurations improves the numerical quality of the fit, improves numerical resolution, and may lift further statistical degeneracies amongst the physical parameters. Moreover, the greater the number of different configurations, the less likely the scenario of overfitting. One or more correlation function(s) may be comprised of: one-point correlation functions, two-point correlation functions, and / or three-point correlation functions. Using such low-order correlation functions enables the experimental sampling to be performed in an affordable measurement time and / or numerical sampling to be numerically affordable. The one or more correlation function(s) may comprise fourth- or higher-order correlation functions. A fourth-order correlation function is a four-point correlation function. In some embodiments, using four-point correlation functions may result in higher precision / resolution in fitting parameters, and in particular may further help to distinguish so-called low-order degeneracy which are not a fundamental non-identifiability (wherein the chosen set of correlation functions is not sufficient to estimate all the parameters e.g. indistinguishable due to the error tolerances of physical measurement). At least one of the physically prepared dynamical configurations may be a steady state of the quantum system, and wherein repeating said steps of physically preparing and physically measuring the quantum system in the steady state a plurality of times may comprise: physically maintaining, with the one or more signal generator(s), the quantum system in the steady state for a single trajectory; physically measuring, with the detector, a continuous signal from the quantum system in the steady state to provide a steady state measured continuous signal; and dividing the steady state measured continuous signal into a plurality of sub-portions, and calculating the corresponding one or more correlation function(s) from the plurality of sub-portions. As will be appreciated, this may advantageously reduce the physical measurement time. By at least one of the physically prepared dynamical configurations being a steady state of the quantum system means that the quantum system is physically prepared in a dynamical configuration which is a steady state. For instance, the first dynamical configuration may be a steady state, and / or the second dynamical configuration may be a steady state, and / or one or more of the third or further dynamical configurations may be a steady state. At least one of the physically prepared dynamical configurations may be a steady state of the quantum system such that the stochastic master equation describing said at least one of the dynamical configurations has a time-constant Liouvillian the method may comprise: defining a plurality of ^^ time-bins of width Δ^^ from the measured continuous signal from said at least one of the dynamical configurations, wherein the width Δ^^ is a digitisation time of the acquisition chain of the detector; wherein numerically fitting the simulations of the one or more correlation function(s) of the steady state, the one or more correlation function(s) being ^^-point correlation function(s), comprises calculating the trace of an operator obtained by numerically evolving a numerically simulated initial state of the quantum system ^^0with a Lindblad evolution interspersed with the application of a correlation superoperator at the center ^^^′^of each ^^-th time-bin from = wherein ^^ is the jump operator of the stochastic master equation describing said at least one of the dynamical configurations, ^^′^^ = ^^Δ^^ + Δ^^ / 2,and ^^^^(^^) is the correlation superoperator for an operator ^^ and defined either by: (i)^^^^(^^) = ^^^^ + ^^^^^^^^† if the detector has binary output with dark count rate ^^ ≥ 0 anddetector efficiency and 0 < ^^ ≤ 1, such that the stochastic master equation describingthe steady state is of the jump type; or (ii) ^^^^(^^) = √^^(^^^^ the detector hascontinuous output with detector efficiency and 0 < ^^ ≤ 1, such that the stochasticmaster equation describing the steady state is of the diffuse type. This approximation is particularly numerically efficient, with the numerical cost being the same a numerically solving a single Lindblad master equation, up to the final “correlation time” ^^^′^being the center of the final time-bin. Alternatively, at least one of the physically prepared dynamical configurations need not be a steady state of the quantum system, but rather a given prepared non- steady state dynamical configuration. In this case, that the stochastic master equation describing said given prepared non-steady state dynamical configuration has a time- dependent Liouvillian having a particular form at each ^^-th time-bin from ^^1′to ^^^′^(ℒ^^1′,…,ℒ^^^′^ ), such that the simulated correlation function is given by ^^ …^^^^ ^^1′ ℒ1^^ ], with ^^ ^ ^^ is the correlation superoperator being ^^^^′1(0)^^^^′ ^( )also calculated at each ^^-th time-bin in view of any time dependence of the jump operators, dark count rates, or detector efficiencies. At least one of the physically prepared dynamical configurations may be a steady state of the quantum system such that the stochastic master equation describing said at least one of the dynamical configurations has a time-constant Liouvillian the method may comprise: defining a plurality of ^^ time-points from the measured continuous signal from said at least one of the dynamical configurations; wherein numerically fitting the simulations of the one or more correlation function(s) of the steady state, the one or more correlation function(s) being ^^-point correlationfunction(s), comprises computing the ^^-dimensional integral of ^^[^^ … against transfer function ^^^^(^^^^) of the acquisition chain of the detector for the ^^-th time-point, wherein ^^[^^^^1 … ^^^^^^] is obtainedby calculating the trace of an operator obtained by numerically evolving a numerically simulated initial state of the quantum system ^^0with a Lindblad evolution interspersed with the application of a correlation superoperator at each time-point from to ^^^^: … ^^^^^^^^1ℒ(^^0)], wherein ^^ is the jump operator of the stochastic master equation describing said at least one of the dynamical configurations, and ^^^^(^^)is the correlation superoperator for an operator ^^ and defined either by: (i)^^^^(^^) = ^^^^ + ^^^^^^^^† if the detector has binary output with dark count rate ^^ ≥ 0 anddetector efficiency and 0 < ^^ ≤ 1, such that the stochastic master equation describingthe steady state is of the jump type; the detector hascontinuous output with detector efficiency and 0 < ^^ ≤ 1, such that the stochasticmaster equation describing the steady state is of the diffuse type. In particular, this provides numerically exact results, particularly suitable for numerically simulating one-point and two-point correlation functions. The sharp signal ^^-dimensional correlation function ^^[^^^^1 … may for instancebe integrated against a rectangular approximation of the . That is wherein transfer function ^^^^of the acquisition chain of the detector at time-point ^^ is approximated by a rectangular function of width Δ^^, wherein the width Δ^^ is a digitisation time of the acquisition chain of the detector. Integrating against the rectangular approximation may for instance comprise integrating using Gaussian quadratures. Alternatively, at least one of the physically prepared dynamical configurations need not be a steady state of the quantum system, but rather a given prepared non- steady state dynamical configuration. In this case, that the stochastic master equation describing said given prepared non-steady state dynamical configuration has a time- dependent Liouvillian having a particular form at each ^^-th time-point from ^^1to ^^^^correlation function is given = , with ^^^^^^^^(^^) is the correlation superoperator being also calculated at each ^^-th time-point in view of any time dependence of the jump operators, dark count rates, or detector efficiencies. In embodiments, numerically fitting at least one of the simulations of the one or more correlation function(s) of a given one of the dynamical configurations, the one or more correlation function(s) being ^^-point correlation function(s), comprises: defining a plurality of ^^ time-points from the corresponding measured continuous signal from the given dynamical configuration; assigning a transfer function ^^^^of the acquisition chainof the detector at time-point ^^; defining ^^ = realnumbers; deriving a system of 2^^coupled ordinary differential equations by forward differentiating an ordinary differential equation ^^^^ ^^ / ^^^^ = ℒ^^ ^^ ^^ such that the system ^^^^1,…,^^ of 2^^coupled ordinary differential equations is describable by:^^^^^^ = wherein ^^ ^^ ^ is a generating density matrix at time-point ^^, and ℒ ^^ ^^^is a superoperator defined either by: (i) ℒ ^^ ^^ = ℒ^^ + (^^ ^^^^ − 1)^^^^ if the stochastic masterequation describing the given dynamical configuration is of the jump type, ^^^^(^^)is acorrelation superoperator for an operator ^^ defined ^^ (^^) = ^^^^ + ^^^^^^^^†if the detector^^has binary output with dark count rate ^^ ≥ 0 and detector efficiency and 0 < ^^ ≤ 1, or^^2= ℒ + ^^ ^^ + ^^(ii) ℒ / 2 ℐ if the stochastic master equation describing the given^^ ^^ ^^ ^^^^dynamical configuration is of the diffusive type, with ℐ the identify operator, and ^^ ^^( )^^( ) is a correlation superoperator for an operator ^^ defined ^^ ^^ = ^^(^^^^the√^^ detector has continuous output with detector efficiency and 0 < ^^ ≤ 1; solving the^^ system of 2 coupled ordinary differential equations from time ^^ = 0 to time ^^, where ^^is chosen to be greater than any time in the support of ^^ , … , to thereby 1,…,^^and calculating the trace of ^^ to provide the simulations of the one or^^correlation function This is a particularly numerically efficient approach, and for low-order correlation functions has a numerical cost which is the same as solving a few Lindblad master equations up to time ^^. The transfer function ^^ of the acquisition chain of the detector at time-point ^^^^may be approximated by a rectangular function of width Δ^^, wherein the width Δ^^ is a digitisation time of the acquisition chain of the detector. In embodiments wherein transfer function ^^ of the acquisition chain of the^^detector at time-point ^^ is approximated by a rectangular function of width Δ^^, the height of the rectangular function may be defined by gain ^^ of the detector acquisition chain. This parameter gain ^^ may be estimated independently from the characterization method described herein (i.e. without explicitly preparing more than one dynamical configuration). For instance, estimating gain ^^ may comprise: preparing a single dynamical configuration by letting the quantum system relax to vacuum state (e.g., by not applying any control signals at all with the one or more signal generator(s)); detecting a continuous signal with the detector; and fitting autocorrelation functions of the measured signal (wherein an autocorrelation function is correlation function of the signal with itself at different times). Thus, the method may comprise independently estimating gain ^^, wherein the gain ^^ is therefore not a physical parameter of the shared subset of physical parameters. The quantum system, in one of the physically prepared dynamical configurations, may be described by a Liouvillian ℒ. Moreover, the quantum system may be physically prepared in said dynamical configuration such that has a symmetry superoperator ^^which commutes with the system’s Liouvillian [ℒ, ^^] = 0 and which anti-commutes withthe signal correlation superooperator {^^^^ , ^^} = 0. In this embodiment, the methodcomprises: calculating only one or more non-null correlation function(s) from said dynamical configuration which are even-order correlation functions; and numerically simulating from the stochastic master equation of said dynamical configuration only one or more even-order correlation function(s). For instance, the symmetry superoperator may be the parity superoperator ^^ ≡^^(^^) in the case of the cat-qubit systems where ^^ denotes the annihilation operator of the physical oscillatory mode a. In embodiments, for each of the one or more qubit-hosting structure(s), the at least one portion is a resonant portion configured to have at least one resonant physical oscillatory mode, each of the one or more qubit-hosting structure(s) further comprises a non-linear element coupled to the at least one resonant portion so as to non-linearly couple to the at least one resonant physical oscillatory mode, and wherein the one or more signal generator(s) are coupled to the at least one resonant portion and / or to the non-linear element, and the detector is coupled to the at least one resonant portion or to the non-linear element and configured to measure continuous signals output from the at least one resonant physical oscillatory mode to provide measurement signals from the quantum system. The quantum system may be for hosting cat qubit; wherein at least one of the one or more qubit-hosting structure(s) is a cat-qubit-hosting structure; wherein the at least one resonant portion is configured to have a first resonant physical oscillatory mode and a second resonant physical oscillatory mode; wherein the non-linear element coupled to the at least one resonant portion is configured to non-linearly couple the first resonant physical oscillatory mode and the second resonant physical oscillatory mode such that, under input control signals to the non-linear element, a two-to-one boson exchange is engineered between the first resonant physical oscillatory mode and the second resonant physical oscillatory mode; wherein a signal generator of the one or more signal generator(s) is coupled to the at least one resonant portion, said signal generator being configured to drive the second resonant physical oscillatory mode; and wherein the detector is coupled to the at least one resonant portion and configured to measure continuous signals output from the second resonant physical oscillatory mode. As will be understood, by engineering a two-to-one boson exchange between the first resonant physical oscillatory mode and the second resonant physical oscillatory mode, under further input control signals to the at least one resonator portion, the quantum system may be physically prepared in a dynamical configuration wherein a cat- qubit is hosted in the first resonant physical oscillatory mode of the least one resonant portion. That is, the first mode may be the “memory” mode discussed herein, and the second mode may be the “buffer” mode discussed herein. Thus, said signal generator and said detector share the same input / output line, thereby reducing the number of control or transmission lines required to operate the quantum system, and further extending the lifetime of the first resonant physical oscillatory mode (that which may host a cat qubit) by reducing the amount of coupling of the first resonant physical oscillatory mode to the external environment. Physically preparing, with the one or more signal generator(s), the quantum system in the first dynamical configuration may comprise physically preparing a cat qubit hosted in the first resonant physical oscillatory mode; wherein step (iii) comprises calculating only one or more non-null correlation function(s) from the first dynamical configuration which are even-order correlation functions; and wherein, step (vii), numerically fitting the first simulation of the one or more first correlation function(s) comprises numerically simulating only one or more even-order correlation function(s). That is, the present Inventors have recognised that a symmetry exists in a cat- qubit system such that only even-order correlation functions are non-null, such that the present characterization method can be made more efficient by only numerically simulating even-order correlation functions. Physically preparing different ones of the dynamical configurations may comprise inputting, with the one or more signal generator(s), respectively different control signals to the quantum system, wherein the different control signals are configured to drive the first resonant physical oscillatory mode and / or second resonant physical oscillatory mode. The present Inventors have discovered that changing dynamical configurations by inputting different drives on the first “memory” physical oscillatory mode, the second “buffer” physical oscillatory mode, or both, is particularly convenient, as signal generator configured to drive these modes are already required for stabilizing and manipulating a cat-qubit (e.g. via a so-called Zeno drive to perform a Z-gate requires inputting a drive on the first “memory” physical oscillatory mode). Thus, the present invention enables full characterization of the cat-qubit system in a single characterization process, using only input / output lines already required for eventually manipulating the cat qubit. Of course, any number of different dynamical configurations may be prepared, for instance different dynamical configurations may comprise stabilizing a cat-qubit having different “sizes”, namely average boson (e.g. photon) number hosted in the first “memory” physical oscillatory mode. The non-linear element may be an asymmetrically threaded SQUID (ATS) comprised of a first superconducting loop and a second superconducting loop, wherein a signal generator of the one or more signal generator(s) is coupled to the ATS and configured to thread the first and second superconducting loops together with a common magnetic flux ^^^^and a differential magnetic flux ^^Δ; and wherein physically preparing different ones of the dynamical configurations comprises inputting, with said signal generator, respectively different pairs of common ^^^^and differential ^^Δmagnetic fluxes(^^^^, ^^Δ), and wherein at least one of the pairs is not at a saddle point (^^^^ , ^^Δ) ≠ . The quantum system may comprise a plurality of qubit-hosting structures, wherein each of the qubit-hosting structures is respectively coupled to at least one other of the qubit-hosting structures by a multi-qubit linear coupling element or by a multi-qubit non-linear coupling element; wherein the method may comprise: in at least one of physically prepared dynamical configurations, physically performing, with the one or more signal generator(s), a multi-qubit gate operation between a qubit hosted in one of the qubit-hosting structures and at least one other qubit hosted in at least one other of the qubit-hosting structures; wherein at least one physical parameter of the shared subset of physical parameters is a physical parameter which characterizes the multi- qubit gate operation. For instance, in embodiments wherein a given one of the qubit-hosting structures is a cat qubit-hosting structure configured to host a cat qubit, the at least one resonant portion of the cat qubit-hosting structure may be coupled to the at least one resonant portion of another qubit-hosting structure by a multi-qubit linear coupling element such as a capacitive coupler or an inductive coupler. The multi-qubit gate operation may be a CNOT gate, and a physical parameter which characterizes the multi-qubit gate operation may be selected from the group comprising: an amplitude or phase of a CNOT pump drive ^^^^^^at the control qubit frequency; an amplitude or phase of a coupling strength ^^^^^^; an amplitude or phase of a linearly compensating counter drive; a CNOT gate operation time ^^^^^^; and a coupling strength of the multi-qubit linear or non-linear coupling element. According to some embodiments, at least one, and optionally each, of the one or more qubit hosting structure(s) comprises a plurality of detectors, each of the detectors being respectively coupled to one of the at least one resonant portion or to the non-linear element so as to measure signals output from the first and / or second physical oscillatory modes to provide measurement signals from the quantum system. Calculating a given set of one or more correlation function(s) may comprise calculating from a plurality of measured continuous signals physically measured with different ones of the plurality of detectors. In a second aspect, there is also provided a quantum characterization unit comprising a quantum system and a command circuit coupled to the quantum system. The quantum system comprises one or more qubit-hosting structure(s), each of the one or more qubit-hosting structure(s) comprising: (I) at least one portion configured to have at least one physical oscillatory mode, (II) one or more signal generator(s) coupled to the at least one portion and configured to input control signals to physically prepare the quantum system in a dynamical configuration, and (III) a detector coupled to the at least one resonant portion and configured to measure continuous signals output from the at least one oscillatory mode to provide measurement signals from the quantum system. The command circuit is configured to: (i) physically prepare, with the one or more signal generator(s), the quantum system in a first dynamical configuration, wherein the first dynamical configuration is describable via a first stochastic master equation having a first plurality of physical parameters; (ii) physically measure, with the detector, a continuous signal from the quantum system in the first dynamical configuration to provide a first measured continuous signal; (iii) repeat steps (i) and (ii) a plurality of times to provide a plurality of first measured continuous signals, and calculate one or more first correlation function(s) from the plurality of first measured continuous signals; (iv) physically prepare, with the one or more signal generator(s), the quantum system in a second dynamical configuration which is different to the first physical configuration, wherein the second dynamical configuration is describable via a second stochastic master equation having a second plurality of physical parameters, and wherein the quantum system is physically prepared in the second dynamical configuration such that a subset of the second plurality of physical parameters is a shared subset with a subset of the first plurality of physical parameters; (v) physically measure, with the detector, a continuous signal from the quantum system in the second dynamical configuration to provide a second measured continuous signal; (vi) repeat steps (iv) and (v) a plurality of times to provide a plurality of second measured continuous signals, and calculating one or more second correlation function(s) from the plurality of second measured continuous signals; (vii) assign values to the shared subset of physical parameters which characterizes the quantum system by jointly numerically fitting: (a) a first simulation of the one or more first correlation function(s), numerically simulated from the first stochastic master equation, to the calculated one or more first correlation function(s), and (b) a second simulation of the one or more second correlation function(s), numerically simulated from the second stochastic master equation, to the one or more second correlation function(s). The second aspect may comprise any of the features and embodiments described above in relation to the first aspect. BRIEF DESCRIPTION OF THE DRAWINGS Various embodiments of the invention will now be described, by way of example only, and with reference to the accompanying drawings in which: Figure 1 shows a diagram of a generic quantum system to be characterized according to an embodiment of the invention; Figure 2 shows components of a quantum system for stabilizing a cat qubit by implementing a parametric dissipative stabilization with an ATS, according to an embodiment of the invention; Figure 3A shows components of a quantum system for stabilizing a cat qubit arranged to perform a resonant dissipative stabilization of a cat qubit, according to an embodiment of the invention; Figure 3B shows an exemplary three-wave mixing nonlinear element which may be used in the embodiment of Figure 3A; Figure 4 shows components of a quantum system on which a CNOT gate may be implemented, according to an embodiment of the invention; Figure 5 shows two longitudinal coupling establishment schemes for a CNOT gate, according to an embodiment of the invention; Figure 6 shows a flow diagram of method steps for characterizing a quantum system, according to an embodiment of the invention; Figure 7 shows how the filtering and digitisation of a sharp diffusive signal is modelled, according to an embodiment of the invention; Figure 8 shows how a three-point correlation function is calculated from physical measurements of filtered or binned signal, according to an embodiment of the invention; Figure 9 shows a method for extracting physical parameters by fitting correlation functions, according to an embodiment of the invention; Figures 10A-10E shows the systems of coupled ODEs to solve, in matrix form, to compute correlation functions up to order three for a single signal for jump or diffusive SME, according to an embodiment of the invention; Figure 11 shows two-point correlation functions for three different combinations of parameters for an exemplary quantum system in a dynamical configuration, according to an embodiment; Figure 12 shows two-point correlation functions for an exemplary quantum system in two different dynamical configurations, according to an embodiment; Figure 13A shows one-point correlation functions of X- and P-quadrature binned signals for an exemplary quantum system in three different dynamical configurations, according to an embodiment; Figure 13B shows a table which summarises the fitted parameters resulting from Figure 13A which characterizes the exemplary quantum system, according to an embodiment; and Figure 14 shows how the one-point correlation of the modelled sharp signal in the example of Figure 13A varies with each parameter. DETAILED DESCRIPTION Quantum system to be characterized Figure 1 shows an exemplary quantum system 1 which is arranged to stabilize or host one or more qubits. The quantum system 1 comprises one or more qubit-hosting structures 6, each of which is configured to host or stabilize a qubit therein. Figure 1 shows two qubit-hosting structures 6 as an example, which are coupled to each other via a linear or non-linear coupler 4. As will be appreciated, there could instead be a single qubit-hosting structure 6 or there could be more than two. Each qubit-hosting structure 6 has at least one resonant portion 9 which is configured to host at least one physical oscillatory mode. In embodiments this could be a first physical oscillatory mode ^^ and a second physical oscillatory mode ^^. As will be understood, the oscillatory mode(s) could be any physical dynamical oscillation suitable for hosting a qubit, such as electromagnetic, magnonic (spin), phononic (acoustic), plasmonic etc. The at least one resonant portion is a physical structure configured to enable the physical oscillatory mode(s) to live and oscillate therein, such as a superconducting resonator in which electromagnetic wave(s) oscillate, a variety of electrodes and lasers configured to trap ions in potential wells defined thereby. Coupled to the at least one resonant portion 9 there may be a non-linear element 7, which for instance couples to the physical oscillatory mode. The at least one resonant portion 9 and the non-linear element 7 coupled thereto together thus provide a non-linear quantum circuit 3. Each qubit-hosting structure 6 also has one or more signal generator(s) 11,13 coupled to the at least one resonant portion 9. In Figure 1, two signal generators are shown both coupled to the at least one resonant portion 9, however there could be only a single signal generator, or more than two. Additionally or alternatively, one or more signal generator(s) could be coupled to the non-linear element 7. Regardless of which part of the non-linear circuit 3 they are coupled to, the one or more signal generator(s) 11, 13 are configured to physically prepare the quantum system 1 in a dynamical configuration by inputting control signals to the at least one resonant portion 9 and / or the non-linear element 7. Preparing the quantum system 1 in a dynamical configuration could be by preparing a qubit of a particular state in at least one of the one or more qubit-hosting structures 6. Similarly, this could be preparing some coherent state hosted in one of the at least one physical oscillatory mode. What is important is that the one or more signal generator(s) 11, 13 enable inputs to the quantum system 1 such that quantum system 1 is a dynamical system. That is, in general, a new dynamical configuration is achieved by either changing the initial state of the quantum system 1 and / or changing the control signals input thereto. Moreover, as will be appreciated, the specific form of the control signals will depend on the specific architecture of quantum system 1. For instance, as discussed more below, if the quantum system 1 is based on superconducting circuits, the control signals may be electromagnetic pulses and drives (in particular microwave radiation), or current biases, or DC voltage biases etc. Figure 1 also shows a detector 8 coupled to the at least one resonator portion 9 in each of the qubit-hosting structures. The detector 9 could instead be coupled to the non-linear element 7. What is important, is that the detector 9 is configured to measure continuous signals output from the at least one oscillatory mode hosted by the at least one resonant portion 9 to provide measurement signals from the quantum system 1. The detector 9 may be of any suitable type, for instance the detector could have a discrete- valued output (e.g. a binary output zero or 1), or the detector 9 could have a continuous -valued output. Figure 1 further shows that each qubit-hosting structure 6 comprises a command circuit 5, which includes a classical processor (i.e. conventional computer) configured to control the signal generator(s) 11,13 and receive the output from the detector 8. Of course, a single command circuit 5 may globally control all of the elements of the various structures of quantum system 1, such as all of the qubit-hosting structures 6 and the non- linear couplers 4 where appropriate, as indicated via the dashed lines in Figure 1. The command circuit 5 is therefore not “quantum”, but rather controls the elements such as the signal generator(s) and the detector which themselves may be coupled to or have some influence on the qubit states hosted in the one or more qubit- hosting structures 6. In this sense then, the command circuit 5 and the quantum system 1 may together be considered quantum characterization unit. The quantum characterization unit may itself be part of a quantum processor, a quantum communication system, or a quantum detection system, as the case may be according to the specific technical use the qubits are put to. The method of the present invention may be applied to a host of different quantum systems, provided a time-continuous signal output from the system is observable by a detector. For instance, the quantum system could be a driven anharmonic oscillator under heterodyne detection. Such a quantum system could for instance be realised by… Another example could be a driven two-level system (i.e. a qubit) whose loss channel is monitored by a photodetector. Such a quantum system could be realised by… The present Inventors have discovered that the present invention is particularly suited to characterising so-called cat-qubits, which typically face difficulties in swiftly and accurately extracting physical parameters which characterize the quantum system. A two-legged cat qubit is defined as a two-dimensional manifold spanned by the so-called cat states |^^^±^ which are superpositions of two coherent states |^^^ and |−^^^: |^^±^ ^^^ = Stabilized cat qubits are known to benefit from a high noise bias, which means that the bit-flip probability is exponentially smaller than the phase-flip probability. More precisely, an effective error channel (e.g., bit errors or “bit-flips”) is suppressed in anexponential way with the “size” – i.e. the average number of photons ^̅^ = |^^|2 – of theSchrödinger cat states of the cat qubits. As previously mentioned, this exponential suppression of bit-flips is only at the cost of linear increase of phase-flips. The cat qubits can be stabilized or confined by the following exemplary schemes: (A) A parametric dissipative stabilization, with jump operator ^^2 = √^^2(^^2 − ^^2),where ^^2is the two-photon dissipation rate, ^^ is the photon annihilation operator of the memory mode a and ^^ is a complex number defining the cat qubit. This jump operator can be realized by coupling a lossy buffer mode b with dissipation rate ^^^^, and a four- wave mixing non-linear element – typically a Josephson junction or an asymmetrically threaded SQUID (ATS) – to the memory mode a and by engineering the Hamiltonian c., where ^^ is the photon annihilation operator of the buffermode b and ^^2is the two-photon coupling rate, by applying to the four-wave mixing non-linear element a pump at frequency |2^^^^ − ^^^^| and a drive of the buffer mode b atfrequency ^^^^ provided ^^2 < ^^^^.(B) A Kerr Hamiltonian where ^^ is the amplitude of the Kerr Hamiltonian, ^^ is the photon annihilation operator, and |^^|2is the mean photon number. (C) A detuned Kerr Hamiltonian where ^^is the amplitude of the Kerr Hamiltonian, ^^ is the photon annihilation operator, ^^ is a complex number defining the cat qubit, and ^^ is the detuning factor. (D) A two-photon exchange (TPE) Hamiltonian ^^⁄ ℏ = ^^2(^^2 − ^^2)^^+ + h. c.,where ^^2is the complex two-photon coupling rate, ^^ is the photon annihilation operator, ^^ is a complex number defining the cat qubit, and ^^±are the lowering and raising operators of the two-level system. This Hamiltonian can be engineered in the same way as the parametric dissipative stabilization (A). (E) A dissipative squeezing stabilization, with jump operator ^^ =^^^^+ sinh(^^)^^^^^^ †2 ^^ )− ^^2), where ^^^^ is the squeezed two-photon^^ dissipation rate, ^^ is the photon annihilation operator of the memory mode a, ^^ is a complex number defining the cat qubit, ^^ and ^^ are the modulus and argument of thecomplex squeezing parameter ^^ = ^^^^^^^^. This jump operator can be realized by couplinga lossy buffer mode b with dissipation rate ^^^^, and a four-wave mixing non-linear element – typically a Josephson junction or an ATS –, to the memory mode a and by engineeringthe Hamiltonian = ^^^^^^ ((cosh(^^)^^ + sinh(^^) ^^† + h. c., where ^^ is thephoton annihilation operator of the buffer mode b and ^^^^^^is the squeezed two-photoncoupling rate, with several pumps at frequencies |2^^^^ − ^^^^|, ^^^^ and 2^^^^ + ^^^^, and a driveof the buffer mode b at frequency ^^^^ provided ^^^^^^ < ^^^^.(F) A variant of the previous stabilization scheme e), for which the Applicant filed the European patent application EP 23175147.0, in which a bosonic qubit – called “moon cat qubit” since the two blobs of the Wigner function have a crescent moon shape – isstabilized by engineering the Hamiltonian ^^⁄ ℏ = − ^^2)^^† + h. c., where g2is the amplitude of the two-photon coupling rate as described above achieved via a pumpat frequency |2^^^^ − ^^^^|, ^^ is the annihilation operator of the memory mode a, ^^ is acomplex number which phase and amplitude result from the amplitude of longitudinal coupling produced by a pump at frequency ^^^^, ^^ is a complex number resulting from a drive of the buffer mode b at frequency ^^^^and ^^ is the annihilation operator of the buffer mode b; a comparison between the moon cat qubit and the squeezed cat qubit could beestablished by expressing ^^ as a function of the complex squeezing parameter ^^ = ^^^^^^^^as follows: ^^ = 2tanh (^^).(G) A DC dissipative stabilization, for which the Applicant filed the European patent application EP 23306839.4, in which a cat qubit is stabilized in the spirit of the previous stabilization scheme a), except that the two-photon pump that engineers the non-linear conversion between two photons of the memory mode a and one photon of the buffer mode b is replaced with a DC voltage source which biases a non-linear element formed exclusively of one or more Josephson junctions, such that the DC-biased non-linear element acts as a voltage-to-frequency converter which, at the appropriatevoltage bias, provides the required parametric interaction at the frequency |2^^^^ − ^^^^|necessary to achieve dissipative stabilization. In particular, the two-photon coupling rate ^^ g^^2is therefore not limited by the amplitude of the two-photon pump: ^^2= 4 ^^^2^^^^^, where ^^^^is the Josephson energy of the one or more Josephson junctions, ^^^^is the zero-point fluctuation of the phase of the memory mode a, and ^^^^is the zero-point fluctuation of the phase of the buffer mode b. (H) A resonant dissipative stabilization, for which the Applicant filed the Europeanpatent application EP 21306965.1, with jump operator ^^2 = √^^2(^^2 − ^^2), where ^^2 isthe two-photon dissipation rate, ^^ is the photon annihilation operator of the memory mode a and ^^ is a complex number defining the cat qubit. This jump operator can be realized by coupling a lossy buffer mode b with dissipation rate ^^^^, and a three-wave mixing non-linear element to the memory mode a which engineers the Hamiltonian =^^2 + h. c., where ^^ is the photon annihilation operator of the buffer mode bprovided the mode frequencies verify substantially 2^^^^ = ^^^^ and ^^2 < ^^^^ to which a driveof the buffer mode at frequency ^^^^is added. In the context of Figure 1, in embodiments wherein the quantum system 1 is for hosting a cat qubit via dissipative stabilization, the quantum circuit 3 may be a circuit fabricated out of superconducting material comprising of various resonators, transmission lines, capacitors, inductors, and Josephson junctions, and may be arranged to make possible three- or four-wave mixing between a first physical oscillatory mode a and a second physical oscillatory mode b. In the following, the first physical oscillatory mode a is used as a memory hosting a cat qubit, while the second physical oscillatory mode b is used as a buffer in between the cat qubit and the external environment. Hereafter, the identifier “physical oscillatory” will be dropped when referring to modes. The first mode a and the second mode b each correspond to natural resonant frequencies of the non-linear superconducting quantum circuit 3. Consequently, the first mode a and the second mode b each have a respective resonant frequency. The first mode a has a resonant frequency ^^^^= and the second mode b has a resonant frequency ^^^^= where ^^^^and ^^^^are the respective angular frequencies of the first mode a and the second mode b. By “having” a first mode and a second mode, it should be understood here that the non-linear superconducting quantum circuit 3 comprises components operating in a superconducting regime which host the modes independently of each other or concurrently. In other words, the first mode a and the second mode b may be hosted in different subsets of components of the superconducting circuit or on the same subset of components. The memory mode a has a high-quality factor ^^^^while the buffer mode b has a low-quality factor ^^^^. As will be appreciated, the quality factor (Q-factor) can be determined in various ways, for instance by: (a) spectroscopic linewidth measurement,wherein the quality factor is given by ^^ = ^^ / Δ^^, wherein ^^ is resonant frequency and Δ^^is the spectroscopic linewidth; or (b) via a time-domain measurement wherein a tone issent in, and a return signal is measured after a pre-determined time, wherein ^^ = ^^ ∗ ^^with ^^ the characteristic decay time. Of course, the skilled person would be aware of various other methods for determining the quality factor of a particular mode. The non-linear superconducting quantum circuit 3 is intended to be subject to electromagnetic radiation delivered by the command circuit 5 in order to engineer various non-linear interactions between the first mode a and the second mode b. The frequency of each electromagnetic radiation is tuned to select specific terms within the rotating wave approximation. Typically, for superconducting materials such as aluminium, the electromagnetic radiation used herein may be in the microwave regime. The non-linear superconducting quantum circuit 3 comprises a non-linear element 7 and at least one resonant portion 9. The non-linear element 7 may be a four- wave mixing non-linear element, such as an ATS or Josephson junction suitable for parametric dissipative stabilization schemes such as scheme (A) above, or one or more Josephson junctions suitable for a DC dissipative stabilization scheme corresponding to scheme (G) above. The non-linear element 7 can also be a three-wave mixing non-linear element, for example a superconducting non-linear asymmetric inductive element (SNAIL) or other three-wave mixing non-linear elements in the case of a resonant dissipative stabilization which corresponds to stabilization scheme (H) above. The resonant portion 9 is arranged to be connected to the non-linear element 7 to provide the non-linear superconducting quantum circuit 3 with the first mode a and the second mode b having respective resonant frequencies ^^^^and ^^^^. More particularly, the first mode a and the second mode b “participate” in the non-linear element 7, which means that a portion or the entirety of the mode magnetic energy is stored in the non- linear element 7. Such a participation can be quantified by the zero-point fluctuation of the superconducting phase across the non-linear element 7, noted ^^^^for the first mode a and ^^^^for the second mode b. In the schematic diagram of the quantum system 1 illustrated in Figure 1, the non-linear superconducting quantum system 3 comprises only one resonant portion, i.e. the resonant portion 9. This may be achieved via design of the resonator (e.g. various capacitances and inductances as defined for instance by the resonator geometry) such that it hosts modes of different frequencies. The design may also result in the modes having different engineered quality factors. However, in general it should be understood here that the non-linear superconducting quantum system 3 comprises at least one resonant portion, and typically two or more resonant portions to form the two electromagnetic modes, which act as a cat qubit “memory” mode and a buffer mode. The first resonator potion may be provided by a 3D resonator, such as a cavity machined out of high purity aluminium (e.g. above 99.99% purity). Similarly, the second resonator portion may provided as an antenna-like cavity. Alternatively, the first and second resonator potions may be provided by a 2D resonator, such as a co-planar waveguide design. Some resonator portions may be 3D, some may be 2D. In either case, such superconducting resonators are considered electromagnetic resonators and the modes hosted therein are electromagnetic modes. Although electromagnetic resonators are convenient for design and control purposes, the resonator portions need not necessarily be electromagnetic resonators. In particular, for the first resonator portion which hosts the cat qubit, it could be a nanomechanical resonator, such as phononic-crystal-defect resonators (PCDRs), similar to those described in P. Arrangoiz-Arriola, E. A. Wollack, Z. Wang, M. Pechal, W. Jiang, T. P. McKenna, J. D. Witmer, R. Van Laer, and A. H. Safavi-Naeini's work, "Resolving the energy levels of a nanomechanical oscillator," published in Nature 571, 537 (2019). These PCDRs are periodically patterned suspended nanostructures that support localized acoustic resonances in the gigahertz range. They are made from a piezoelectric material like LiNbO3, enabling the coupling of these resonances to superconducting circuits with nearly the same efficiency as standard electromagnetic cavities. Specifically, quasi-one-dimensional PCDRs made from lithium niobate, a piezoelectric crystalline material, could be used, with modes localized within a volume less than 1 μm3 of a suspended nanostructure. In such acoustic resonators, the cat qubit is encoded in the phonons of the resonator. Other types of acoustic resonators, such as acoustic membranes, also exist. The primary requirement is an acoustic-to-electromagnetic coupler, to couple the acoustic first resonator portion to the non-linear element 7. Furthermore, other types of resonators may be used, such as a magnetic resonator wherein the quanta of oscillations are magnons (collective spin excitations). An example could be a designed yttrium-ion-garnet (YIG) particle, as described in Kounalakis, Marios, Gerrit EW Bauer, and Yaroslav M. Blanter. "Analog quantum control of magnonic cat states on a chip by a superconducting qubit." Physical review letters 129.3 (2022): 037205. Again, there is a requirement to couple magnetic first resonator portion to the non-linear element 7. Yet another type of resonator suitable for hosting cat qubits could be semiconductor double quantum dots embedded in a cavity such as a split-ring resonator cavity as disclosed in Kozin, Valerii K., et al. "Quantum phase transitions and cat states in cavity-coupled quantum dots." Physical review research 6.3 (2024): 033188. Thus, herein a cat qubit may be encoded in the bosons of the resonator. Herein, for the sake of simplification, the 2-to-1 boson exchange is a 2-to-1 photon exchange, however it will be appreciated that all reference to photons and photon exchange regarding the first mode which hosts the cat qubit can be replaced with bosons or boson exchange (e.g. acoustic phonons, or magnons, or indeed plasmons). The command circuit 5 is arranged to deliver control signals, such as electromagnetic radiation of DC voltage biases. To this end, as illustrated in figure 1, the command circuit 5 comprises signal generators in the form of at least a first mode drive 11 and a second mode drive 13. For instance, the command circuit 5 may be arranged at least to drive the second mode b by delivering radiation at a frequency substantially equal to the second resonant frequency ^^^^to the resonant portion 9. By “a frequency substantially equal to”, it should be understood that, ideally, the frequency is exactly equal to the desired value. However, in practice, the frequency value deviates from the desired value, typically by 1 or even 5%, due to the inherent precision of the hardware used. Moreover, as will be understood, the frequency of a mode may be shifted in use. For instance, under certain dynamical interactions such as application of an electromagnetic pump to the resonant portion 9, the frequencies of the first mode a and the second mode b may be effectively shifted (e.g. Stark shifted) due to the presence of a static potential or dynamical interaction. However, the skilled person would be well aware of characterising such effective frequencies, and for simplicity herein, by the first mode a and the second mode b “having” respective resonant frequencies ^^^^and ^^^^, it will be understood as static meaning “un-dressed” frequencies or shifted “dressed” frequencies when in the presence of a dynamical interaction, depending on the situation. Figure 2 illustrates an embodiment in which a single qubit-hosting structure 6 is arranged to perform the parametric dissipative stabilization to stabilize a cat qubit (i.e. a cat-qubit-hosting structure 6), wherein the non-linear element 7 is an ATS. An ATS is an inductive dipole element formed by a pair of Josephson junctions (indicated by the crosses and ^^^^label in Figure 2) being shunted by an inductance (indicated by the ^^^^ label in Figure 2). The inductance may be formed for instance from a linear array of Josephon junctions. The Josephson junctions are each arranged on a different superconducting path such that they each form a distinct loop with the inductance, such that the ATS can be viewed as two superconducting loops connected by the shared edge which comprises the inductance. The ATS is thus sensitive to the magnetic flux threaded through each of the two loops. The at least one resonant portion 9 is coupled via a linear coupler 15 to the ATS 7 which acts as an inductive element, such that the first mode a and the second mode b participate in the ATS 7. Example of embodiments for the linear networks which form the at least one resonant portion are known. For instance, among these embodiments, the non-linear superconducting quantum circuit 3 can be a two-mode hybridized system (also known as "galvanic cat", for which the Applicant has filed patent application EP22306815.6 and EP22306816.4) for which two lumped modes couple strongly via the ATS. It is known to the skilled person that an ATS can be used to engineer in general a 2-to-1 boson conversion, i.e. the first term ^^2of the jump operator in order to perform a parametric dissipative stabilization More particularly, a 2-to-1 photon conversion between the first mode a and the second mode b is obtained byparametrically pumping the ATS at the frequency ^^^^ = |2^^^^ − ^^^^|, as successfully shownby Lescanne, Raphaël, et al. "Exponential suppression of bit-flips in a qubit encoded in an oscillator." Nature Physics 16.5 (2020): 509-513. Advantageously, the pumpfrequency satisfies ^^^^ ≫ ^^2 to make this parametric pumping work as well as possible.It is known to the skilled person that, in the case of the non-linear element 7 beingan ATS 7, when biased at its flux working point also known as the “saddle point” (0 − ^^,namely 0 flux threaded through one of the loops and ^^ flux in unit of the flux quantum threaded through the other of the loops), the Hamiltonian of the ATS 7 has the following “sin-sin” form: where ^^ = ^^^^(^^ + ^^†) + ^^^^(^^ + ^^†) is the total superconducting phase difference across the ATS 7, ^^^^is the zero-point fluctuation of the phase of the first mode a across the ATS 7 and ^^^^is the zero-point fluctuation of the phase of the second mode b across the ATS 7, ^^^^is the Josephson energy of the side junctions and ^^^^is the inductive energy of a central inductance, ^^^^(^^) corresponds to a common flux modulation of the two loops of the ATS 7 and ^^Δ(^^)corresponds to a differential flux modulation of the two loops of the ATS 7. The former can be implemented by delivering electromagnetic radiations through flux lines of the ATS 7 out of phase, the latter can be implemented by delivering electromagnetic radiations through the two flux lines of the ATS 7 in phase. Parametric pumping of the ATS 7 is typically done by pumping the common flux as pumping the differential flux merely displaces the modes coupled to the ATS 7. As anexample, by pumping the common flux at the frequency ^^^^ = |2^^^^ − ^^^^|, ^^^^(^^) =^^2^^ℎcos (2^^^^^^^^), the non-linear resonant part of the Hamiltonian writes in the rotatingframe: ^^^^^^^^ = + h. c. ), which is typically the two-to-one photonexchange Hamiltonian needed to engineer the two-photon stabilization. When an external DC magnetic field is set such that a 0 mod 2^^ magnetic flux threads one of the loop and a ^^ mod 2^^ flux threads the other loop, it is ensured that the ATS Hamiltonian has its “sin-sin” form. For clarity, the set-up applying the external DC magnetic field is not drawn on Figure 1 but can be applied via the two bottom mutual inductances of the ATS 7. A typical implementation consists in interleaving a bias-tee connected to a DC current source to input DC current into the system while letting the electromagnetic radiations go through. As will be described below, the components which generate the required magnetic fields are current-carrying flux bias lines close to the ATS 7. The cat-qubit-hosting structure 6 comprises a further signal generator in the form of an electromagnetic source 17 which is set-up to modulate the common flux in the ATS 7. For this purpose, a electromagnetic network 19 is used to split the radiation emitted by the electromagnetic source 17 and apply it to each node of the ATS 7 with the correct phase. Alternatively, two different electromagnetic sources could be used, each being simply coupled to a single node of the ATS 7 and their relative phase and amplitude being set so as to achieve the desired flux modulation. When the electromagnetic source 17 is set at the frequency |2^^^^ − ^^^^|, the non-linear superconducting quantum circuit 3 performs the 2-to-1 photon conversion between the first mode a and the second mode b. To convert this 2-to-1 photon conversion into two-photon dissipation, the second mode b is selectively coupled to a load 21 via a linear coupler 23 and a electromagnetic filter 25 configured as a band pass filter with a frequency ^^^^. Alternatively, the electromagnetic filter 25 may be configured as a band stop filter at a frequency ^^^^and may be placed in between, on the one hand, the external environment and, on the other hand, the first mode a and the second mode b to isolate the first mode a and thus prevent the first mode a from suffering additional losses coming from unwanted coupling to the load 21. Alternatively, it may be configured as a low-pass (respectively high-pass) filter if^^^^ > ^^^^ (resp ^^^^ > ^^^^). In other embodiments, the electromagnetic filter 25 can be omittedwhen coupling between the load 21 and substantially only the second mode b can be established. The skilled person thus understands that the first mode a has a high-quality factor while the second mode b has a low-quality factor. As previously explained, the second mode b is driven at its resonant frequency ^^^^. This two-photon drive is performed by a electromagnetic source 13 set at frequency ^^^^. In the above, the load 21 can be seen as part of the command circuit 5 of figure 2, while the linear coupler 23 and the electromagnetic filter 25 can be seen as part of the non-linear superconducting quantum circuit 3. The command circuit 5 of cat-qubit-hosting structure 6 may further comprise signal-generators in the form of a electromagnetic source 11 and a electromagnetic source 27. In addition to stabilizing the cat qubit, the command circuit 5 also enables measurement an observable of the cat qubit or application of a quantum gate to the cat qubit. The electromagnetic source 11 may be arranged to drive the first mode a by delivering electromagnetic radiation at the frequency ^^^^to the linear electromagnetic network b / a. Such a drive of the first mode a causes the non-linear superconducting quantum circuit 3 to engineer a Hamiltonian HZ expressed as ^^^^⁄ ℏ = ^^^^^^ + h. c., wherethe complex rate ^^^^results from the amplitude and phase of the drive of the first mode a. An additional signal generator in the form of another electromagnetic source may be set-up, in addition to the electromagnetic source 17, to contribute to the modulation of the common flux in the ATS 7. The additional electromagnetic source may be arranged to deliver electromagnetic radiation for causing the non-linear superconducting quantum circuit 3 to engineer a Hamiltonian. For instance, such a Hamiltonian may yield a so- called longitudinal coupling term between the first mode a and the second mode b, which may be used to apply a Z gate. The electromagnetic source 27 is set-up to modulate the differential flux in the ATS 7 through the electromagnetic network 19. The electromagnetic source 27 can be used to reduce spurious Hamiltonian terms induced by the additional electromagnetic source. It is to be noted that the electromagnetic network 19 is used for convenience but can be omitted and the electromagnetic source 27 and any additional electromagnetic source could be applied directly to both nodes of the ATS 7, their relative phase and amplitude being set so as to achieve the desired flux modulation. For the sake of completeness, it may also be noted that the electromagnetic source 27 can be used, instead of the electromagnetic source 13, to deliver electromagnetic radiation at a frequency ^^^^to the linear electromagnetic network b / a to drive the second mode b. Typically, the circuit has flux lines through which radiation can be delivered to provide (and optionally modulate) the common flux (herein also referred to as ^^^^or “sigma flux”) and provide (and optionally modulate) the differential flux (herein also referred to as ^^Δor “delta flux”). These are generally referred to as the two bias modes. As such, the flux lines may deliver both a DC-bias (setting the working point described above by applying DC current to thread the ATS loops with the relevant magnetic fluxes) and an RF-bias (superposed modulation to permit the time-dependent magnetic flux terms to “pump” the ATS as discussed above). Finally, detector 8 is shown as being connected to the electromagnetic source 13 which delivers electromagnetic radiation at a frequency ^^^^to the linear electromagnetic network b / a. Thus, as will be appreciated, a single line (e.g. transmission line) may be used to both input electromagnetic radiation and output a continuous signal. In this case, the continuous signal is output from the second mode b. In particular, the detector 8 may be configured for heterodryne detection of the fluorescence from the second mode b. As will be appreciated, other implementations of the detector is possible, for instance it could primarily be coupled for detection of the first mode a. Figure 3A illustrates an embodiment in which the quantum system 1 is arranged to perform resonant dissipative stabilization to stabilize a cat qubit, i.e. wherein the non- linear element 7 is a three-wave mixing non-linear element 7. Thus, in comparison with Figure 1, the non-linear element 7 and the at least one resonant portion 9 are brought together to form part of the non-linear superconducting quantum circuit 3, which explains why the reference signs “7” and “9” are absent from Figure 3A. Here, the non-linear superconducting quantum circuit 3 is arranged to perform intrinsically the 2-to-1 photon exchange – symbolized here by a back-and-forth single arrow 57 and double arrows 59 – between the first mode a and the second mode b (denoted by reference numerals 29 and 31 to indicate hosting of the two modes). Like reference numerals indicate features which are the same as described above in Figure 2. For instance, the electromagnetic radiation sources 11 and 13 are configured in the same manner as in Figure 2 to respectively drive the first mode a and the second mode b. The command circuit 5 further comprises a current source 61 or applying a current to the non-linear element 7 and the at least one resonant portion 9. The current source 61 may be connected by wires (i.e. galvanically) to the non-linear element 7 and the at least one resonant portion 9, or may be a flux line for flux biasing a loop (e.g. of the non-linear element 7 as in Figure 3B) to thereby induce a current flowing around the loop due to the magnetic flux therethrough. The current source 61 is arranged to deliver current which flows through one or more components of the non-linear superconducting quantum circuit 3. The current source 61 is configured to both allow three-wave mixinginteraction and to tune the frequency matching condition 2^^^^ = ^^^^, in particular byinducing a particular current flow across the three-wave mixing non-linear element 7 which enables tuning of the resonant frequencies of first mode a and the second mode b as they participate in the non-linear element 7. Detector 8 in Figure 3A is shown as being coupled to the part 29 of the at least one resonant portion which predominantly hosts the first mode a. Specifically, a tomography transmon 10 may be coupled (e.g. capacitively coupled) to the part 29 of the at least one resonant portion, to which the detector 8 is coupled to read the state of the tomography transmon (e.g. reading a continuous signal through an associated intermediate readout resonator and Purcell filter). Of course, other implementations of the detector 8 again is possible, for instance it could primarily be coupled via the same line which inputs from the electromagnetic source 13 in a similar manner as in Fig.2. Figure 3B shows a possible circuit to realize the three-wave mixing non-linear element 7, which in this embodiment is formed by at least one loop 63 including a first Josephson junction 65, a central inductive element 67 and a second Josephson junction 69. Alternative embodiments of the three-wave mixing non-linear element 7 are described in EP21306965.1 filed by the Applicant. Returning to Figures 3A and 3B, when a predetermined current of a constant intensity is applied by the current source 61, the resonant frequency ^^^^is substantially equal to twice the resonant frequency ^^^^. Said otherwise, what is important is that the superconducting circuit 3 has these structural features so linked, and that (regardless of the specific shape of the components) the resonance matching condition is satisfied when the constant current is applied. This current can be determined by the skilled person any number of known ways. In particular, what is important is that the non-linear element comprises at least one loop 63 including at least one Josephson junction therein, and that the resonant modes a and b at least partially participate in the loop. The circuit shown in Figure 3B is particularly configured to discriminate symmetrically the memory mode a and the buffer mode b. The high symmetry of such a circuit achieves an improved quality of the 2-to-1 photon exchange, as well as to make the function of filter 25 integral to the flux line used to flux bias the loop 63 (which is also the second mode electromagnetic source 13). The relative position with respect to the loop 63, and the geometry, of the flux line which acts as the second mode electromagnetic source 13 enables coupling to the second mode whilst inherently filtering the first mode to prevent or substantially minimize coupling of the first mode a to the load (environment) 21. The central inductive element 67 can be an inductance, a single Josephson junction or an array of Josephson junctions. The central inductive element 67 may thus be arranged between the first Josephson junction 65 and the second Josephson junction 69 as a loop in series. The arrangement in series may comprise a first inner node connecting a pole of the first Josephson junction 65 with a pole of the central inductive element 67. The arrangement in series may also comprise a second inner node connecting a pole of the second Josephson junction 69 with another pole of the central inductive element 67. The arrangement in series may also comprise a closed-loop node connecting another pole of the first Josephson junction 65 with another pole of the second Josephson junction 69. The at least one loop 63 may be connected to a common ground via the closed- loop node. The circuit may also comprise a first capacitive element 71 and a second capacitive element 73. The first capacitive element 71 may be connected in parallel with the first Josephson junction 65 between the common ground and the first inner node of the loop. The second capacitive element 73 may be connected in parallel with the second Josephson junction 69 between the common ground and the second inner node of the loop. The first Josephson junction 65 and the second Josephson junction 69 are substantially identical and the capacitive elements 71 and 73 are also substantially identical. Hence, the symmetry of the circuit implies that the memory mode a is the symmetric superposition of the two resonators (as shown by the full arrows) and the buffer mode b is the anti-symmetric superposition of the two resonators (as shown by the dashed arrows). It can be noticed that only the buffer mode b has a contribution across the central inductive element 67 which is advantageously used to preferentially couple the external environment to this buffer mode b while isolating the memory mode a from the external environment via the design of the electromagnetic source 13 as described above. Figure 4 shows a physical system on which a CNOT gate may be implemented. Here, the quantum system 1 comprises a target cat-qubit-hosting structure 300 (which is similar to that of Figure 2) and a control qubit-hosting structure 302 connected by a linear electromagnetic coupler 304. In the example described herein, the target cat-qubit hosting structure 300 comprises a non-linear superconducting circuit 306 to which are connected several signal generators in the form of electromagnetic (e.g. microwave) sources 310, 311, 316 and a load 314. In a similar manner as described above in relation to Figure 2, the non-linear superconducting circuit 306 comprises at least one resonant portion (denoted by reference numerals 320 and 322 to indicate hosting of the two modes) such that when coupled via a linear coupler 309 to an ATS 308 which acts as an inductive element the non-linear superconducting circuit comprises at least 2 normal modes (or eigenmodes) a and b at frequency ^^^^and ^^^^which participate in the ATS. As above, this participation means that a portion or the entirety of the mode magnetic energy is stored in the ATS. This participation can be quantified by the zero-point fluctuation of the superconducting phase across the ATS, noted ^^^^for first mode a and ^^^^for second mode b. In the context of the present embodiment, the physical realization of a CNOT gates between a control qubit with annihilation operator ^^ and a stabilized cat qubit with annihilation operator ^^ usually relies on the use of the following two ingredients. (1) The confinement of the control qubit such that a microwave drive at its resonant frequency results in a Rabi oscillation. This is typically native in two-level system qubits (e.g. such as transmons) and engineered via parametric interactions for cat qubits. (2) The addition in the circuit of a ‘CNOT’ Hamiltonian or ‘longitudinal’ Hamiltonianwith the formula ^^^^^^⁄ ℏ = ^^^^^^(^^ where ^^^^^^ is the amplitude of theHamiltonian (it is chosen real without loss of generality). This longitudinal coupling can be seen as a drive on the control qubit (first factor) which amplitude depends on the photon number of the target cat qubit (second factor). For this Hamiltonian to be effective on the target cat qubit, one also needs to turn-off the confinement on the target cat qubit, hence the target cat qubit confinement strength should be controllable. To engineer the CNOT Hamiltonian between the control qubit and the target cat qubit one needs to: (i) couple the control qubit to the ATS such that the phase differenceacross the ATS writes ^^ = ^^ + ^^†^^ ) + ^^^^(^^ + ^^†) ; and (ii) pump thecommon flux at the control qubit frequency ^^, = ^^ ^^^^^^ (2^^^^^^) with ^^ ^^^^ ^^ ^^^^amplitude of the pump drive. In the rotating frame, the parametric part of the Hamiltonian writes ^^ =^^^^^^†−^^^^ ^^ ^^ + ^^(. This Hamiltonian can be written to highlight the desired dynamics)^^ ^^^^ ^^Since second mode b (the “buffer”) is a lossy mode coupled to a cold environment, one can †further assume that ^^ ^^ = 0, such that the engineered Hamiltonian is3 ^^^^† 2 2 †2 † 2( ) ^^ = ^^ − ^^^^ ^^ (^^ + ^^ ) 1 − ^^ ^^ + ^^^^ (^^ ^^ + ^^ ^^ ).^^^^^^ ^^^^ ^^ ^^^^ ^^ ^^ ^^ ^^^^2 This Hamiltonian is close to the CNOT Hamiltonian except for 2 additional terms. 22 ⁄ ( ) The first corresponds to a linear drive with strength ^^^^ ^^ 1 − ^^ ^^ ℏ . This^^ ^^^^ ^^ ^^linear drive can readily be compensated by sending a counter drive directly on the control qubit with the correct relative phase and amplitude that can both be tuned experimentally. The Applicant also found that the accuracy and scope of the compensation can be greatly increased by pumping the differential flux of the ATS with the correct phase and 2^^^^2 2( ) ( )( amplitude which writes, at first order in ^^ , ^^ ^^ = ^^ ^^ 1 − ^^)^^ . Although the^^^^ ^^ ^^ ^^^^^^amplitude and phase of the drive can be computed analytically, it is fine-tuned experimentally by ensuring the control qubit remains undergo no drive in the right circumstances. Experimentally, this compensation is much more appropriate than the counter drive on the control qubit because, it directly compensates the spurious linear drive where it originates from (i.e., at the ATS) and does not only try to compensate its main consequences (i.e., the displacement of the control qubit). The second term cannot typically be fully compensated in a simple manner in view of the design choices made for the present embodiment. Instead, the system may be defined such that the amplitude of this term is much smaller than the amplitude of the ⁄2 ≪CNOT Hamiltonian. In other words, which simplifies into^^ 2⁄. Typically, the Applicant has found that ^^ ≤ ^^ 2 is sufficient minimize the^^^^ ^^ ^^detrimental impacts of this second term. In summary, provided these two conditions (compensation and smaller amplitude for the second term) are met, the ATS Hamiltonian ^^ can be engineered such that it^^^^^^is close enough to the perfect CNOT Hamiltonian. For clarity, the set-up applying the external DC magnetic field is not drawn on Figure 4, but can be applied via the 2 bottom mutual inductances of the ATS. A typical implementation consists in interleaving a bias-tee connected to a DC current source to input DC current into the system while letting the electromagnetic radiations go through. The electromagnetic sources 310 and 311 are set-up to modulate respectively the common and differential flux in the ATS which is required to activate parametric interactions. To clearly distinguish the roles of the two sources in the Hamiltonian, an electromagnetic network 324 which applies the correct phase offset is provided in the schematic. Alternatively, each electromagnetic source can be simply coupled to a single node of the ATS and their relative phase and amplitude can be set so as to get the desired flux modulation. In that case, to modulate the common flux the two sources need to address the circuit out of phase and to modulate the differential flux, the two sources need to address the circuit in phase. The 2-to-1 photon conversion necessary to dissipatively stabilize a cat qubit is performed by setting the frequency of theelectromagnetic source 310 to ^^^^ = |2^^^^ − ^^^^|, and stabilization is further realised bydriving the electromagnetic source 316 at the second mode b frequency ^^^^. When the CNOT gate is not performed, i.e., in a so-called "idle mode", the command circuit controlling the entire system is configured to perform the data (target) cat qubit stabilization exclusively. The control qubit and the configuration required to perform the CNOT gate according to an embodiment will now be described. In the example of Figure 4, the control qubit-hosting structure 302 comprises a resonant portion 305 which has a mode ^^ with resonant frequency ^^^^which hosts the control qubit. Although not shown, qubit-hosting structure 302 will inevitably also include some non-linear element. In various embodiments, the control qubit-hosting structure 302 can be any superconducting qubit hosted in a qubit-hosting structure such as a transmon qubit, a flux-qubit or a fluxonium qubit or any bosonic qubit encoded in a resonator such as a Kerr cat qubit (detuned or not), another dissipative stabilized cat qubit device (squeezed or not), or a cat qubit confined via a two-photon exchange Hamiltonian. As described above, the control qubit 305 is coupled to the target cat-qubit- hosting structure 300 via a small linear coupler 304. The linear coupler 304 is arranged such that the control qubit 305 slightly hybridizes with the cat qubit device 300 which leads to a small participation of the control qubit in the cat qubit device ATS 308. This participation is denoted ^^^^. As described above with respect to the second spurious member of the engineered Hamiltonian which is not compensated, the fact that this participation remains small compared to the participation ^^^^of the target cat qubit mode is critical to accurately implement the CNOT Hamiltonian. In various embodiments, coupler 304 can be capacitive, inductive, galvanic or mediated via a resonating bus coupler or an additional linear electromagnetic network. In the embodiment shown in Figure 4, the control qubit-hosting structure 302 is coupled to the first mode ^^ 320. In other embodiments, control qubit-hosting structure 302 can be coupled to the second mode ^^ 322. This coupling location is not critical as mode ^^ and ^^ are typically delocalized in the electromagnetic network comprising the at least one resonant portions^^⁄ ^^ 320-322 and coupling to a specific location does not necessarily mean coupling toa specific mode unless the linear microwave network and the ATS are specifically designed to do so. The control qubit-hosting structure 302 also comprises a signal generator in the form of control signal source 303 which is coupled to the control qubit 305 in order to have the ability to drive it. For instance, source 303 may be an electromagnetic radiation source, which may be used to compensate the linear drive arising from the CNOT Hamiltonian engineering. The control qubit-hosting structure 302 also comprises a detector 8 similar to those described above, here shown coupled to resonant portion 305. As explained above, a further limitation on the control qubit-hosting structure 302 is that the non-linear superconducting circuit 306, the first mode (hereinafter cat qubit mode) ^^ participates strongly in the ATS 308 as compared to the control qubit mode 305. When the control qubit mode ^^ participates into the ATS via a small linear coupling with the cat qubit mode ^^ in order to perform the CNOT gate – typically, capacitive coupling with capacitance value small compared to the either mode capacitance, inductive coupling with inductance value small compared to the either mode inductance or coupling mediated by a detuned bus resonator –, this coupling is moderate in general. More precisely, the two coupled modes (here ^^ and ^^, alternatively ^^ and ^^) detuning(i.e., the frequency difference ^^ = |^^^^ − ^^^^|) also plays a role. Indeed, if the two modeshave the same resonant frequency, any slight linear coupling would lead to fullhybridization and to ^^^^ ≈ ^^^^. However, with typical detunings ^^⁄ 2 ^^ greater than a fewtens of MHz, the participation asymmetry is achieved with standard linear couplings. As explained earlier, the weaker coupling of the control qubit ^^ to the ATS 308 (compared to the target cat qubit mode ^^ coupling to the ATS 308) is key. If the linear coupling is characterized by a strength ^^ as known in the field, then one also needs to ensure the detuning ^^ between the control qubit mode ^^ and the mode it is coupled to(^^,320 or ^^,322) is such that ^^ < ^^. Although this provides design rules, a fullelectromagnetic (e.g. microwave) simulation or circuit diagonalization of the circuit layoutis required to precisely compute the values of ^^^^ , ^^^^ , ^^^^. This feature helps to ensurethat the CNOT gate can be performed accurately and that the circuit uses a single ATS in its elements for the CNOT gate, and that parametric pumping of the ATS of the target cat qubit has a single purpose at all times. Figure 5 shows exemplary signal timings for operations 800, 810 and 820. During idle mode, the data cat qubit is stabilized with two-photon dissipation ^^2. When the gate starts, this dissipation is turned off such that the CNOT Hamiltonian can be effective.Finally, after a time = ^^⁄ (4^^^^^^^^) the CNOT Hamiltonian is turned off and the two-photon dissipation is turned back on again. At no point in time does the ATS serve two purposes at once, which guarantees experimental robustness. The limitations for the CNOT pulse strength and shape (dotted or dashed lines) will now be explained. The CNOT Hamiltonian effectively acts as a linear drive on the control qubit ^^ which strength depends on the number of photons in the target cat qubit ^^. The nature of the control qubit ^^ sets upper bounds on the maximal effective drive strength ^^^^^^^^. There exist two main cases. In the first case, the control qubit ^^ is defined by areal or Hamiltonian gap – for a two-level system this gap is the anharmonicity |^^12 − ^^^^|where ^^^^is the qubit frequency or the frequency of ground to first excited state transition and ^^12is the frequency of first excited to second excited state transition (this may for instance be a transmon). In that case, the adiabatic theorem applies and states that the qubit stays exponentially confined despite the action of the effective drive provided^^^^^^^^ < |^^12 − ^^^^|. For instance, if the ancilla qubit is a cat qubit with amplitude ^^^^confined by a Kerr Hamiltonian (scheme (B) above) the condition writes ^^^^^^^^ < ^^^^^2^. If it is confined by a detuned Kerr Hamiltonian (scheme (C) above) the condition writes^^^^^^^^ < ^^(^^^2^ + √^^). If it is confined by a TPE Hamiltonian (scheme (D) above), thecondition writes ^^^^^^^^ < ^^^^^^2.In the second case, the control qubit ^^ is defined by an imaginary gap or equivalently is stabilized by dissipation. In that case, the adiabatic theorem does not apply and the gate strength has to be much smaller than the imaginary gap in order to avoid undesired decoherence of the ancilla qubit during the gate. In practice, for an ancilla cat qubit stabilized by two-photon dissipation with rate ^^2(a) the condition writes^^^^^^^^ ≪ ^^2^^^2^, which applies also to schemes (G) and (H) above. For a squeezed ancillacat qubit (scheme (E) above), the condition writes ^^^^^^^^ ≪ ^^2^^^2^^^2^^where ^^ is the squeezing parameter, and a similar condition may be derived for the “moon cat” stabilization (scheme (F) above). On top of the gate strength requirement, the nature of the control qubit confinement sets constraints on how the gate should be applied. In the case of a dissipative ancilla qubit, the CNOT Hamiltonian can be turned on instantaneously as shown by the dashed lines of the above figure as soon as the data cat qubit confinement is turned off (solid line ^^2(^^)). In the case of a Hamiltonian ancilla qubit, the CNOT Hamiltonian may be turned on smoothly such that the spectral content of the pulse ^^^^^^(^^) (dotted lines) does not contain frequency components above the gap. Typically, a gaussian pulse can be used. On the above figure, a cosine shape is used for its finite temporal envelope. In the case of a time dependent ^^^^^^, the pulse amplitude should besuch that ∫ 4^^^^^^^^(^^)^^^^ = ^^.As will be appreciated, other CNOT methods may be applied. Furthermore, other multi-qubit gates, such as a Toffoli gate between or multi-qubit CNOT or other gates may be implemented involving two or more qubit-hosting structures 6 or quantum system 1. For such systems and operations, be it for single- or multi-qubit systems, the various parameters describing the coupling strengths and control signals required for various dynamical configurations need to be characterized. Characterization method The present Inventors have developed a method for characterizing a quantum system 1, whether it be comprised of a single qubit-hosting structure or multiple qubit- hosting structures such as those outlined in the above embodiments, and which crucially operates as a non-linear dynamical system. Figure 6 shows a flow diagram of the method steps for characterizing a quantum system 1 which is operable as a non-linear dynamical system, for instance as described above in relation to Figures 1-4. In a first step 601, the one or more signal generators 11,13 are used to physically prepare the quantum system 1 in a first dynamical configuration. The evolution of a continuously measured quantum system is modelled by the stochastic master equation (SME) formalism. The SME describes the evolution of the quantum state and the corresponding signal ^^^^measured by the observer at time ^^. Specifically, an observer describes the state of an open quantum system 1 by a mixed quantum state, the density matrix ^^^^at time ^^. Mathematically, we associate with the system a Hilbert space of dimension (possibly infinite), and the quantum stateis a Hermitian positive semi-definite operator with unit trace: ≥0, ^^^^[^^] = 1}. Different laws of evolution describe the trajectory of the quantum state{^^}^^∈[0,^^]in state space In the Markovian scenario there are two situations: if the observer cannot observe the quantum systems composing the environment, the state trajectory is described by the Lindblad master equation (ME). Conversely, if the observer measures these external degrees of freedom, the state trajectory is governed by the stochastic master equation (SME). The ME describes the evolution of the average state ^̅^^^ρt when the observer does not monitor the environment, for example for a purely dissipative process or for unread measurements. The evolution is deterministic, it is described by the linearordinary differential equation (ODE): ^^^̅^^^ / ^^^^ = ℒ^^(^̅^^^) = −^^[^^^^, ^̅^^^] + where ℒ^^ is the Liouvillian superoperator, ^^^^ is the Hamiltonian of the system, =^^^^^^† − (1 / 2)^^†^^^^ − (1 / 2)^^^^†^^ is the standard dissipator and {^^1,^^, … , ^^^^,^^} is a collectionof jump operators. For simplicity, we consider in the following a single loss channel with time- independent loss operator ^^. In some situations, this loss channel is fully or partially measured by the observer. The acquired information can be used to update the system state, which gives a better (less mixed on average) description of the state trajectory than the ME. This is precisely the purpose of the SME. The SME describes the evolution of the state ^^^^when the observer continuously measures the loss channel with a detector. The evolution is non-deterministic, it isdescribed by the non-linear stochastic differential equation (SDE): ^^^^^^ = ℒ^^(^^^^)^^^^ + is a superoperator which depends on a stochastic process ^^^^^^:the measured signal. The first part of the equation is the deterministic Lindblad update. The second part models the stochastic back-action of the measurement, which is taken into account continuously, at each time. The detector output is a continuous-time signal defined by the rate of change ofthe stochastic process ^^^^^^ over time: ^^^^ = ^^^^^^ / ^^^^. The path followed by the quantumstate over time is entirely determined by this measured signal: each time the observer performs a new experiment (labelled ^^), they measure a particular realisation of the stochastic signal {^^^^^(^)} , which corresponds to a unique trajectory in state ^^∈ space. These trajectories are called quantum trajectories, and the state ^^^^is said to be conditioned on the information measured by the observer up to time ^^. There are two main types of SMEs: the jump SME and the diffusive SME. The difference lies in the type of detector: if the output is binary (0 or 1), the evolution is described by the jump SME, and otherwise if the output takes a continuous range of values, the evolution is described by the diffusive SME. For example, in quantum optics, the jump SME models photodetection, while the diffusive SME models homodyne or heterodyne detection. Mathematically, they differ by the stochastic process ^^^^^^involved, and the form of the measurement back-action superoperator For the jump SME, most of the time the detector detects nothing and outputs 0, and sometimes, upon detection, the detector clicks, and outputs 1. This stochasticprocess is modelled by the point process ^^^^^^ = ^^^^^^ with law: ℙ[^^^^^^ = 0] = 1 − here ^^ ≥ 0 is the dark countrate (taking into account false clicks), and 0 < ^^ ≤ 1 is the detector efficiency (taking intoaccount missed clicks). The measured signal ^^^^ = ^^^^^^ / ^^^^ is the rate of change of thecounting process ^^ = ∫^^ ^^0^^^^^^, which counts the number of jumps occurring in the timeinterval [0, ^^). The measurement back-action is defined by: ℳ^^(^^, ^^^^) =− ^^)(^^^^ − (^^ + ^^^^^^[^^^^^^†])^^^^), as described for instance in Rouchon, Pierre. "A tutorialintroduction to quantum stochastic master equations based on the qubit / photon system." Annual Reviews in Control 54 (2022): 252-261. For the diffusive SME, the detector output is real-valued, and continuous in time. This stochastic process is given by the Itô process ^^^^^^defined by: ^^^^^^=√^^^^^^[(^^ + ^^†)^^^^]^^^^ + ^^^^^^, where again 0 < ^^ ≤ 1 is the detector efficiency (taking intoaccount measurement imperfections), and ^^^^is a Wiener process taking independent Gaussian distributed increment. The measurement back-action is defined by: described forinstance in Jacobs, Kurt, and Daniel A. Steck. "A straightforward introduction to continuous quantum measurement." Contemporary Physics 47.5 (2006): 279-303. In the jump SME case, the quantum trajectory is discontinuous: the state evolves continuously in state space as long as the detector detects nothing, but upon detection it undergoes a sudden jump. In the diffusive SME case, the state evolves continuously in state space, following a Brownian motion-like trajectory. The unconditioned trajectory described by the ME is recovered by averaging over all possible realisations of the stochastic process driving the SME (or equivalently, over all possible quantumtrajectories or measured signals): ^̅^^^ = ^^[^^^^]. Here ^^ denotes the statistical average over^^^^or ^^^^. Note that although the quantum trajectories described by the jump SME and the diffusive SME are very different in nature, the average state trajectory does not depend on the stochastic process averaged over (jump or diffusive). As will be appreciated, the above continuous-time formulations for the ME and SME can be readily converted into a discrete-time formulation suitable for numerically simulation, for instance via using the well-known Kraus operators formulation. As such, the evolution of quantum system 1 can be modelled by a SME, and specifically the initial state, the form the Hamiltonian and the various jump operators entering the dissipator. The skilled person would be well aware of how to provide such a model for a particular quantum system 1, as described by a few exemplary embodiments below. Importantly, the model of a quantum system 1 is typically known apart from the values of ^^ parameters corresponding to various static or dynamical terms, which canbe gathered in a vector ^^^^^^^^^^ ∈ ℝ^^. These parameters may appear in different parts of the model, such as the initial state, Hamiltonian, jump operators, dark counts or efficiencies, or filter functions as described in more detail below. The first dynamical configuration prepared in step 601 is thus describable via a first stochastic master equation having a first plurality of physical parameters ^^ (^^) ^^^^^^^^, where the superscript in parentheses merely indicates that this refers to the first dynamical configuration. However, the continuous signal ^^^^in relation to the SME is never available to the operator of quantum system 1, it is merely a mathematical object, and not a tangible, measurable quantity. It can be loosely thought of as a series of Dirac delta distributions for the jump SME, and as white noise with a trend for the diffusive SME, and this continuous-time signal is referred to as the sharp signal herein. The actual measured signal is processed by an acquisition chain consisting of amplifiers, filters and digitisation components (e.g. an analogue-to-digital converter, ADC, or a photocounter) generally comprised by detector 8. This chain converts thecontinuous-time sharp signal ^^^^ into a discrete-time signal {^^0, … , ^^^^}. Each value ^^^^ isdefined by integrating ^^^^against the transfer function ^^^^of the acquisition chain for the^^-th time bin: ^^^^ = ∫ ^^^^(^^)^^^^^^^^. The only quantity available to the experimentalist is thisdiscretised measurement record {^^0, … , ^^^^}. It is particularly important to take this filteringinto account when the digitisation time is not negligible compared to the system timescales, as is often the case for superconducting circuits, for example. The digitisation is usually performed by averaging the signal over a duration Δ^^ (a time bin), which is longer than the bandwidth of the various acquisition chain components. In embodiments, the filter function for the ^^-th time bin can then be approximated by a rectangular window ^^^^of duration Δ^^. In the case of the jump SME, the detector typically gives the number of clickevents over a given time interval, and the filter function ^^^^ is the indicatorfunction defined by ^^Ω(^^) = 1 if ^^ ∈ Ω and ^^Ω(^^) = 0 otherwise, wherein […) defines thehal-open interval such that ^^ ∈ [0,1) means 0 ≤ ^^ < 1. The resulting signal takes discretevalues ^^^^ ∈ ℕ: ^^^^ = where the counting process =∫^^2^^1 ^^^^^^ counts the number of jumps occurring in the time interval [^^1, ^^2). In the limit wherethere is at most one jump per time bin Δ^^, the signal is binary ^^^^ ∈ {0,1}, it is just asequence of 0s and 1s: 00100010… , and the filtering accounts for the inevitably finite time resolution of the detector. In the case of the diffusive SME, the signal is typically averaged over someduration Δ^^, the filter function is ^^^^ = ^^[^^Δ^^,(^^+1)Δ^^) with ^^ the gain of the acquisitionchain. The resulting signal is continuous-valued ^^(^^+1)Δ^^^^ ∈ ℝ: ^^^^ = (^^ / Δ^^) ∫^^Δ^^^^^^^^^^. In both case, we call this digitised signal the binned signal. From a practical point of view, this discrete-time signal ^^^^is the only quantity available to the operator of quantum system 1. This digitisation process is illustrated in Figure 7 for the diffusive SME, wherein the filtering and digitisation of the sharp diffusive signal ^^^^against the rectangular windowfilter functions ^^^^ in Figure 7 results in a discrete-time binned signal {^^0, … , Returning to Figure 6, in a second step 603, the detector 8 is thus used to physically measure a time-continuous signal from the quantum system 1 in the first dynamical configuration, which provides a first measured time-discrete signal. The statistics of this time-discrete signal can be estimated experimentally by performing the same experiment several times and averaging over the different realisations. These statistics include the mean of the signal ^^[^^ ] at time bin ^^, the two-^^point correlation function ^^[^^ ^^ ] between bins ^^ and ^^ , and more generally any n-^^ ^^ 1 21 2point correlation functions ^^[^^ … ^^ ] for time bins ^^ ≤ ⋯ ≤ ^^ , where ^^ denotes the^^ ^^ 1 ^^1 ^^statistical average over the stochastic process driving the SME. In particular, the present Inventors have recognised that a stochastic process is characterizable by its correlation functions. Thus, as shown in a next step 605 of Figure 6, the physically preparing 601 and physically measuring steps 603 are repeated a plurality of times to provide a plurality of first measured signals for the first configuration, and one or more first correlation function(s) are calculated from the plurality of first measured continuous signals. Specifically, estimating the correlation functions of the measured signal is straightforward experimentally. For each realisation of the experiment (labelled ^^), the operator has access to a new discrete-time digitised signal … , ^^(^^) ^^}. The ^^-point correlation function is given by: ^^[^^ … ^^^^^^] = lim with ^^ the statistical average over an infinite number of realisations of the experiment. In practice, the correlation functions are estimated with some statistical error due to thefinite number of realisations. We denote ^̂^[^^^^1 … ^^^^^^] the experimental estimate of theexpectation value: Figure 8 illustrates the computation for a three-point correlation function: The present Inventors have also recognised that, in embodiments when starting from the quantum system 1 in a steady-state, a single long trajectory can be measured, instead of repeating the experiment multiple times. The measured signal can then be divided into ^^^^^^^^chunks, and the average is computed by considering each chunk as a different realisation of the experiment. Certain technical considerations must be taken into account for this approach to be rigorous: (i) the steady state must be ergodic under continuous measurement, and (ii) there must be sufficient statistical independence between the different chunks (which is generally ensured by the fact that correlation functions decay to zero at long times). Generally, it is expected that each parameter has a unique impact on the evolution of the state, and therefore on the statistics of the measured signals. To simultaneously estimate the ^^ parameters, first some correlation functions of interest are estimated / measured experimentally (i.e. calculated from physically measured data), and then the correlation functions are numerically simulated from the model of the system (namely the first SME) and fit to these observations, for example with a least-squares method. The result is a vector of parameters ^^ (^^) that minimises the difference between^^^^^^the experimental estimate and the theoretical model. The expectation is that the estimated parameters ^^ (^^) (^^) ^^^^^^are close to the true parameters ^^^^^^^^^^. However, the present Inventors have recognised that the parameters that we wish to estimate may be degenerate or non-identifiable: different combinations of parameters result in the same correlation functions. For instance, consider a quantum system 1 being a lossy harmonic oscillator withHamiltonian ^^ = ^^^^†^^ and jump operator ^^ = √^^ ^^, where ^^ is the oscillator annihilationoperator, is the oscillator frequency, and ^^ is the single-photon loss rate. The loss channel is monitored with efficiency ^^ by homodyne detection along the X quadrature.Suppose the system starts from a coherent state ^^0 = |α0^ ^α0| and our goal is to findthe values of the four parameters ^^ = (^^, ^^, ^^, ^^0). Then it is straightforward to see thatthe parameters ^^ and ^^0are fundamentally non-identifiable from the measured signal.Indeed, the state remains a coherent state during the evolution: ^^^^ = |αt^ ^αt| with αt =^^0^^−^^ / 2^^^^−^^^^^^, because the measurement backaction is null at all time: ℳ(^^^^ , ^^^^^^) = 0.The measured signal allows to uniquely identify ^^ and ^^, but not ^^ and ^^0, as they only appear as a product√^^^^0in the measured signal. The present Inventors have therefore discovered that such a degeneracy or non- identifiability may be addressed by evaluating the correlation functions of the quantum system 1 in another, different, dynamical configuration. Taking the simple example of the lossy harmonic oscillator, placing the quantum system 1 into a second, different, dynamical configuration could for instance be to add a known linear drive, make the oscillator anharmonic, and / or simply prepare quantum system with a different value of the initial state complex amplitude ^^0linked by a known proportionality coefficient (for instance ^^^^0, where ^^ is known). Thus, in a third step 607, the one or more signal generators 11,13 are used to physically prepare the quantum system 1 in a second dynamical configuration different to the first dynamical configuration. Crucially, the second dynamical configuration is describable via a second SME having a second plurality of physical parameters ^^ (^^) ^^^^^^^^, wherein a subset of the second plurality of physical parameters ^^ (^^) ^^^^^^^^is a shared subset with a subset of the first plurality of physical parameters ^^ (^^) ^^^^^^^^, wherein the shared subset of physical parameters characterizes the quantum system 1. In a step 609, the detector 8 is again used to physically measure a continuous signal from the quantum system 1 in the second dynamical configuration, in the same manner as described above in relation to the first dynamical configuration, which provides a second measured continuous signal. For the avoidance of doubt, herein the “first” correlation function(s) should not be interpreted as “one-point” correlation functions, but rather is / are the correlation function(s) of the first dynamical configuration. Thus, these may include the one-point correlation function, but also or alternatively might include higher-order point correlation functions. In a similar vein, the “second” correlation function(s) should not be interpreted as “two-point” correlation functions, but rather is / are the correlation function(s) of the second dynamical configuration (which for instance may include the one-point correlation function of the second dynamical configuration, or indeed the two-point correlation function or any other correlation function). In a step 611, the physically preparing 607 and physically measuring steps 609 are repeated a plurality of times to provide a plurality of second measured signals, and one or more second correlation function(s) are calculated from the plurality of second measured continuous signals, in the same manner as described above in relation to the first dynamical configuration. Again, when starting from the quantum system 1 in a steady-state for the second dynamical configuration, a single long trajectory (i.e. a single acquisition) can be measured and split up into subsections, instead of repeating the experiment multiple times. In a step 613, a plurality of the one or more first correlation function(s) are numerically simulated from the first SME, wherein at least one physical parameter of the first plurality of physical parameters has a different value when numerically simulating each of the simulated one or more first correlation function(s). This numerically enables fitting the one or more first correlation function(s) in a later step. In a step 615, a plurality of the one or more second correlation function(s) are numerically simulated from the second SME, wherein at least one physical parameter of the second plurality of physical parameters has a different value when numerically simulating each of the simulated one or more second correlation function(s). Again, this numerically enables fitting the one or more second correlation function(s) in a later step. Thus, in a step 617, the simulated one or more first (respectively second) correlation functions are fitted to the calculated one or more first (respectively second) correlation functions. For the avoidance of doubt, here “calculated” refers to those correlation functions calculated from physically measured data. By fitting, it will be generally understood that the vector of ^^ parameters ^^ entering the equations used for simulating the one or more first (respectively second) correlation functions are modified to produce particular simulations of the one or more first (respectively second) correlation functions. Any known fitting procedure may then be applied, which compares the simulations to the calculated correlation functions, modifies the vector of parameters and iterates until the simulation has satisfactorily converged to the calculated correlation functions. The vector of parameters may thus be that minimises the difference between the experimental estimate and the theoretical model. Accordingly, values can be assigned to the shared subset of physical parameters which characterizes quantum system 1 and which are shared between the first set of physical parameters ^^ (^^) and the second set of physical parameters ^^ (^^) ^^^^^^^^ ^^^^^^^^. In embodiments, the minimisation may be solved using the Levenberg– Marquardt algorithm using the “least_squares()” function from the SciPy library at Pauli Virtanen et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods, 17 (2020). In embodiments, access to the gradient can greatly improve the fit time and accuracy when estimating many parameters at once. This can be achieved by using modern automatic differentiation libraries to integrate the system of ODEs. For the examples described herein, the present Inventors used the “Tsit5 ODE solver” of the Diffrax library from https: / / github.com / patrick-kidger / diffrax to compute the correlation functions. However, any standard ODE solvers library or methods therein is suitable, especially if the gradient is not needed. The present Inventors have recognized that, in an optional particularly advantageous step 619, steps 601-615 may be repeated a further number of times, and N times in total such that N sets of one or more correlation functions are physically measured and are numerically simulated. Thus in an optional step 621, each of the N different simulated correlation function(s) may be fit with the corresponding N sets of different correlation function(s) which were physically measured so as to assign values to the shared subset of physical parameters which is shared by all of the corresponding (^^) (^^) (^^), ^^ , … , ^^N sets of physical parameters (i.e. N vectors vector of parameters ^^ ).^^^^^^^^ ^^^^^^ ^^^^^^In general, the present Inventors have discovered that the greater number of dynamical configurations (and the more random the differences between each, e.g. the more random the drives in the above example of a driven lossy harmonic oscillator), then the less likely the overfitting procedure. In particular, it is beneficial to prepare, measure, and fit at least 3 different dynamical configurations, and preferably greater than or equal to 10. Figure 9 illustrates the method for extracting physical parameters by fitting correlation function. Approximative or exact formulae may be used to compute these simulated quantities directly from the SME and the initial state of the quantum system 1, as discussed in more detail below. In particular, the present Inventors have recognised particularly numerically efficient methods for numerically simulating correlation functions for fitting physical parameters. An explicit formula for simulating correlation functions from the sharp signal fordifferent times ^^1 < ⋯ < ^^^^ for a constant Liouvillian ℒ^^ = ℒ (i.e. time-constant) can begiven by where ^^^^(^^)is the correlation superoperator for an operator ^^, and is defined either by:(i) ^^^^(^^) = ^^^^ + ^^^^^^^^† for the jump SME, or (ii) ^^^^(^^) = √^^(^^^^ + ^^^^†) for the diffuse SME, wherein the notations take the same forms as defined above. Namely, operator ^^0is the density matrix describing the initial state of quantum system 1, ^^ is the dark count rate, ^^ is the detector efficiency, and ^^ is loss or jump operator (which can be time- independent or take different forms at different times). Said otherwise, the correlation functions can be simulated by taking the trace of a peculiar operator, which is obtained by evolving the initial state with the usual Lindblad evolution but interspersed with the application of the correlation superoperator at each correlation time ^^ .^^The present Inventors have recognised that, if the digitisation time Δ^^ is very small compared to the fastest timescale of the dynamics of quantum system 1, then the correlation functions of the binned signal can be approximated with the sharp signalformula described above ^^. Namely ^^[^^^^1 … ^^^^^^] ≈ ^^[^^^^1′ , wherein the binned signal ^^^^ is replaced with the sharp signal^^ in the center of the time bin ^^′^^^′^ ^^ = ^^Δ^^ + Δ^^ / 2. Importantly, this approximation is notvalid for coincident bins. Moreover, particular care must be taken to ensure that Δ^^ is sufficiently small for the approximation to be correct. A good way to control the validity of the approximation is to compute the exact formula and verify that it matches. The sharp signal correlation functions can be computed by evolving the initial state with the ME using any regular solver known in the art, and applying the correlation superoperator at each correlation time. The numerical cost is the same as solving a single ME, up to the final correlation time ^^ .^^For filtered signals – the correlation function of the filtered signal ^^[^^ … ^^ ] can^^ ^^1 ^^be simulated by explicitly computing the ^^-dimensional integral of ^^[^^ … against corresponding filter functions ^^[^^ … ^^^^ ,^^ ^^1 whereas above when the filter functions have on-overlapping supported. However, this expression quickly becomes prohibitely expensive to compute for large Hilbert space dimensions, and correlation functions involving more than two points. As detailed in Guilmin, Pierre, Pierre Rouchon, and Antoine Tilloy. “Correlation functions for realistic continuous quantum measurement.” IFAC-PapersOnLine 56.2 (2023): 5164-5170, there is a faster and more 5 practical way to compute these quantities by solving modified Lindblad master equations, which we now recall, which the present Inventors have recognised may beneficially be applied in the present characterisation method so as to enable faster simulations required for fitting. Explicitly, we10 ^^ where ^^∞is the so-called generating density matrix at long time, where the generating density matrix at time ^^= ^^ [exp (∫^^ ^^ is given by ^^^^0^^^^′^^^^^^′ ) ^^^^] and ^^^^^^′ = ^^^^′^^^^′. Th ^^e generating density matrix at time ^^, ^^^^, obeys the so-called generalised quantum master equation ^^^^ ^^ / ^^^^ = ℒ^^(^^ ^^), where ℒ ^^ ^^ ^^ ^^ ^^is another superoperator and is ^^= ℒ + (^^ ^^^^ − 1)^^ for t^^he jump SME, or (ii) ℒ = ℒ + ^^ ^^ +15defined either by: (i) ℒ^^ ^^ ^^ ^^ ^^^^ ^^2 ^^ / 2 ℐ for the diffuse SME with ℐ the identity superoperator, wherein like-notations take^^the same forms as defined above. 1,…,^^Thus, to compute ^^ , one can (manually) forward differentiate the ordinary∞differential equation Recursively distributing the derivatives results in an augmented system of coupled linear ODEs (called sensitivity equations), which involve all unordered combinations of partial derivatives with respect to a subset of {^^ , … , ^^ }. Each ODE1 ^^describes the evolution of a fictitious state under the regular Lindblad evolution, with the 25 introduction of additional coupling source terms. The system is solved from time 0 to time ^^, where ^^ is chosen to be greater than any time in the support of Figures 10A-10E give the system of ODEs to solve (in matrix form) to compute correlation functions up to order three for a single signal, wherein for compactness, we drop the time index and mark null entries with a dot. Figure 10A shows the ODE to solve 30 for calculating one-point correlation functions for both the jump and the diffusive SME. For instance, for the one-point correlation function ^^ ^^ this ODEs system is solved from[ ]^^1time ^^ = 0 to time ^^ = ^^ (greater than any time in the support of ^^) by evolving two fictiousstates with initial conditions (^^ , 0), with the one-point correlation function of^^=0filtered signal ^^ being then given by the trace of ^^1^^1 ^^=^^ at final time ^^ = ^^: ^^[^^^^1] = ^^^^[^^1^^]. Figure 10B shows the ODE to solve for calculating two-point correlation functions for the jump SME. Figure 10C shows the ODE to solve for calculating two-point correlation functions for the diffusive SME. Figure 10D shows the ODE to solve for calculating three-point correlation functions for the jump SME. Figure 10E shows the ODE to solve for calculating three-point correlation functions for the diffusive SME. These can be solved in a similar manner as described above in relation to Figure 10A. In particular, the system of ODEs can be viewed as a single large linear ODE, which can be integrated using commonly available ODE solvers, such as high-order Runge-Kutta schemes. Note that for large Hilbert space dimensions, it is not recommended to explicitly store or diagonalize the full matrix as shown in Figures 10A-10E, for both memory and runtime purposes. A better solution is to apply the generator to the state at each time step, which costs only ^^ in memory in time (where is the Hilbert space dimension). For low-order correlation functions, the numerical cost is the same as solving a few MEs up to time ^^. As discussed, binned signals are a special case of filtered signals where the filter functions are rectangular windows. In particular, if the Liouvillian is time-independent, the generator of the large ODE is piecewise constant in time. Thus, it can be integrated by successively exponentiating the generator over each corresponding time interval. For large Hilbert space dimensions, it is not necessary to compute the exponential explicitly, but simply its action on an operator. This can be done efficiently using e.g. Krylov subspace methods. For low-order correlation functions, the numerical cost is the same as solving a few MEs with a time-independent Liouvillian up to time ^^. Moreover, the present Inventors have recognised that simulating correlation functions from binned signals may be particularly numerically advantageous as, for non- coincident time bins, all cross terms involving more than one filtering function (e.g. ^^ ^^1 2shown in Figure 10C) are null. Moreover, the computation can often be efficiently vectorized when computing multiple correlation functions. For example to compute 212[ ] ^^ ^^ ^^ for various time bins 1 ≤ ^^ ≤ ^^, the states ^^and ^^ remain null before time0 ^^[ ] ^^ = ^^Δ^^. The quantities ^^ ^^ ^^ can then be computed in a vectorized fashion by (i)^^ 0 ^^1 evolving the states (^^ , ^^ ) on [0, ^^ ) with the one-point system, then (ii) evolving the^^ ^^ 111( ) resulting states ^^ , ^^ with the Liouvillian on [^^ , ^^ ), and saving the result ^^ , ^^ (at )^^ 1 ^^ ^^^^ ^^^^ ^^intermediate times ^^ , and finally (iii) batch-compute the evolution of the ^^ states^^(^^ , ^^1, ^^2, ^^12) on [^^ , ^^ ) using the two-point system, initialised with the saved results^^ ^^ ^^ ^^ ^^ ^^ ^^+1(^^ , ^^1^^^^ ^^^^ , 0,0) from the previous step.As will be appreciated, the above formulae may be generalised to compute correlation functions between signals from different detectors, that can be of jump and / ordiffusive type. In the above, ^^ is now a set of test functions, one for each detector ^^ =number of measured signals. Splitting detectors between of jump-type ^^ ∈ ^^ and those of diffusive-type ^^ ∈ ^^ , the new generator of^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^ ^^ ^^ / ^^^^ = ℒthe ODE ^^^^ (^^ ) described above is:^^ ^^ ^^ where like notations take the same definitions as before. Numerical examples As discussed above, by preparing and measuring the quantum system 1 in two or more different dynamical configurations, certain parameters can be fit which would otherwise be non-identifiable. For instance, another example of non-identifiability can be found in the system of Berdou, Camille, et al. "One hundred second bit-flip time in a two-photon dissipative oscillator." PRX Quantum 4.2 (2023): 020350, which is a bistable dynamical quantum system comprised of a structure similar to the qubit-hosting structure 6 of Figure 2, but having a heterodyne detector coupled to the first mode a (specifically a transmission line weakly coupled to the “memory” mode resonator). 2 In particular, the present Inventors have recognised that, for a fixed value of ^^ , it is practically impossible to independently estimate ^^ and ^^ from one-, two-, and three-2 2point correlation functions of this system due to the antagonistic combination of (i) the measured signal being a symmetric stochastic process such that only even-order correlation functions are non-null (from a particular symmetry of the system arising from the action of the two-photon dissipation and one-photon loss, as discussed below) and (ii) different combinations of ^^ and ^^ resulting in almost indistinguishable two-point2 2correlation functions due to the intractably high number of trajectories needed to achieve sufficiently low error bars. Here, the dynamical configuration is formed of opposite phase coherent states, when the system has reached its steady state, which are dynamically stabilised in the cavity field (oscillatory mode a) of a cavity (i.e. a resonator portion), wherein the stabilisation is achieved via engineering two-photon dissipation between the cavity andthe environment modelled by jump operator ^^2 = √^^2 (^^2 − ^^22), with two-photon dissipation rate ^^2, complex amplitude ^^2which can be controlled by an external drive,and single-photon loss from the cavity at a rate ^^1 modelled by jump operator ^^1 = ^^,wherein ^^ denotes the annihilation operator of the cavity field. In this reduced model ofthe system which describes just the dynamics of the cavity, the Hamiltonian is ^^ = 0.Figure 11 shows two-point correlation functions ^^[^^0^^^^]for three different combinations of ^^2and ^^2. The parameters ^^2and ^^2cannot be estimated from the two- point correlation functions alone, as the different curves are superimposed so as to be entirely indistinguishable (i.e. statistically degenerate). Although this degeneracy can be lifted by evaluating, e.g., the four-point correlation functions, accurately estimating high-order correlation functions becomes quickly prohibitively expensive in measurement time (due to the much greater standard error of the mean for a fixed number of measurements compared to lower-order correlation functions). Instead, by changing the dynamical configuration, different two-point correlation functions can be resolved, as shown in Figure 12. Note that the “experimental” data used here was in fact modelled via simulating 105stochastic trajectories of a corresponding diffusive SME modelled to mimic experimental conditions, wherein the upper panelshows results for configuration "^^" with ^^2 = 7.0 and the lower panel shows results forconfiguration "^^" with ^^2 = 7.14. Similar “mimicked” experimental data provided fromappropriate models have been used in Figures 11 and 13-14, and are merely provided for clarity in lieu of experimental data acquired in the lab. In particular, for each trajectory (labelled "^^"), we average the measured signalsover 31 time bins of duration Δ^^ = 1 / (2^^), to obtain two time-binned signals, one for eachconfiguration: We estimate the two-point correlationfunctions of the signals for non-coincident time bins ^^[^^0^^^^] (1 ≤ ^^ ≤ 30) by averagingover the 105trajectories. We then jointly fit the estimated correlation functions ^^[^^0^^,^^^^^^^^,^^] with a least-squares method. For this example, we fit the binned signal correlation functions using the sharp signal formula approximation described above. Thus, the present Inventors have recognised that the present invention is particularly suited to characterizing a quantum system 1 comprising one or more cat- qubit-hosting structure(s) 6, as discussed above in relation to Figures 2-4. Each of the one or more cat-qubit-hosting structures 6 comprise at least one resonant portion 9 configured to have a first mode a and a second mode b, and a signal generator 13,316 configured to input control signals to drive the second mode b such that, in combination with the non-linear element 7 and other control signal(s) input to the at least one resonant portion 9 and / or the non-linear element 7, the quantum system 1 is physically prepared in a first and subsequently a second dynamical configuration. For instance, the first and / or dynamical configurations may be a cat-state or a coherent state stabilized or prepared in one or more of the cat-qubit-hosting structures. What is important is that the dynamical configurations are different to each other. Moreover, the present Inventors have recognised that, if the dynamical configuration is prepared using two-photon dissipation and one-photon loss only, all odd- order correlation functions are null. Thus, the first dynamical configuration may comprise preparing a cat-state in the cat-qubit-hosting structure 6 and physically measuring a first continuous signal from the cat-qubit-hosting structure 6. The method may thus comprise calculating only even-order correlation functions from the first dynamical configuration. It may be useful to then break the symmetry which is what causes the null odd- order correlation functions, to also leverage the information in these odd-order correlation functions. Thus, subsequently, the method may comprise adding a symmetry-breaking drive to the quantum system 1, such as a drive on the first mode a, to thereby prepare the quantum system in a second dynamical configuration and physically measure a second continuous signal from the cat-qubit-hosting structure 6. The method may thus comprise calculating both even- and odd-order correlation functions from the second dynamical configuration. The cancellation of odd-order correlation functions can be generalized to any continuously measured system that satisfies conditions (i) that some traceless symmetrysuperoperator ^^ commutes with the system’s Liouvillian [ℒ, ^^] = 0 (such as paritysuperoperator in the case of the above cat-qubit systems where ^^ denotes the annihilation operator of the physical oscillatory mode a), and (ii) itanti-commutes with the signal correlation superoperator {^^^^, ^^} = 0.Returning to quantum systems 1 suitable for stabilizing cat-qubits, as discussed above in relation to Figures 2-4, it is particularly advantageous to integrate the detector 8 with the signal generator used to drive the second mode b (e.g. source 13 in Figures 2-3, and source 316 in Figure 4). This strategy recognises that information which will inherently leak anyway through the transmission line used to transmit signals from source 13,316, and thus is designed to detect and use this information, thereby avoiding the need to place a detector elsewhere (which would introduce another leakage channel) and conveniently reducing the number of lines going to non-linear superconducting quantum circuit 3, 306. Each of the one or more cat-qubit-hosting structures 6 therefore further comprise a detector 8 coupled to the at least one resonant portion 9 and configured to measure continuous signals output from the second mode b. For instance, in embodiments wherein the at least one resonant portion 9 comprises a first 29,320 and a second 31,322 resonant portion configured respectively to substantially host the first mode a and the second mode b, then the detector 8 may be coupled to the second resonant portion 31,322. The present method thus synergistically makes use of the continuous signal inherently output through the transmission line of the signal generator used to drive the second mode b via detecting it using detector 8 couple thereto and fitting correlation functions as described above. Importantly, the present Inventors have recognised that fitting the correlation functions when starting from a system state which is not the steady state gives additional information on the parameters, from fitting the time dynamics of the system. Figure 13A shows one-point correlation functions of X- and P-quadrature binned signals in three different dynamical configurations (denoted by the three different columns). The system is a quantum system 1 for stabilizing a cat-qubit as described above in relation to Figure 2. Here, continuous signal in the form of fluorescence from the second “buffer” mode b output through coupler 23 and filter 25 is detected by detector 8. The top row shows theoretical data ^^^^[^^^^] (i.e. for “true” parameters ^^^^^^^^^^which have been theoretically modelled here), the middle row shows data estimated “measured” data by averaging over many trajectories ^̂^^^[^^^^], and the bottom row shows the fit ^^^^^^^^^^[^^^^] in light shades of the estimated “measured” data ^̂^^^[^^^^] in dark shades (same plot as the middle row). Specifically, Figure 13A shows numerical proof-of-concept, wherein the dynamics of quantum system 1 of Figure 2 are described by: (i) Hamiltonian ^^ which is the idealised ATS Hamiltonian at the so-called saddle point as discussed above plus some additional terms which enter due to the quantum system 1 being non- ideal in reality (due to, e.g. fabrication imperfections and higher-order terms not taken into account when deriving ^^^^^^^^above); and (ii) jump operators describing the loss channels from the first mode a and the second mode b, respectively ^^^^and ^^^^. Namely: where ^^2is the two-photon exchange rate as described above, ^^^^ / ^^are complex drives respectively on the first / second “memory” / ”buffer” modes a / b (with real parts ^^^^^^ / ^^and imaginary parts ^^ , ^^ is the so-called self-Kerr on the first mode a, ^^6is the so-called order-6 self-Kerr on the first mode a, ^^^^is the single-photon loss rate from the first mode a, and ^^^^is the single-photon loss rate from the second mode b. The second mode fluorescence is continuously measured by detector 8 by heterodyne detection along the X- and P-quadratures with efficiency ^^, resulting in two signals ^^^^^ = √^^^^ / 2 ^^^^[^^(^^†^ ^^ − ^^)^^^^]^^^^ + ^^^^^^^^, where ^^^^^^^^and ^^^^^^^^are independent Wiener processes. As will be appreciated, all the parameters can be fitted by the present invention, assuming the amplitude and phase of the drives to be configurable in that they can be switched on or off, scaled by a known factor, and rotated by a known angle. Here, ^^ is chosen real, such that we want to reconstruct the real vector with 10 parameters2^^ = (^^2, ^^ ^^, ^^^^, ^^ ^^, ^^^^ ^^ ^^ ^^ ^^ , ^^, ^^6, ^^^^ , ^^^^ , ^^ ) (for the avoidance of doubt, here ^^ is a vector ofparameters, and not the dark count rate described above in relation to the jump SME). The amplification and digitisation chain of detector 8 converts the continuous- time analogue sharp signal ^^ ^^ / ^^ ^^into the discrete-time digitised signal, which here we take as binned signal ^^ ^^ / ^^ ^^with the ^^-th transfer function ^^^^approximated by rectangular window ^^ ^^ / ^^ ^^ Here the window is of duration Δ^^ = 4 ns with the gain^^ of the acquisition chain being set to unity for this example (and which can be estimated independently by any number of known methods, such as by fitting the autocorrelationof the signals when the drives are switched off and ^^ = 0).2This heterodyne measurement is modelled by splitting the jump operator ^^^^into two parts with halved rates ^^ / 2, and considering the two signals as the result of^^homodyne measurement of each loss channel ^^^^ = √^^^^ / 2 ^^ and ^^ ^^^^ ^^ = √^^^^ / 2 (−^^^^),which gives the same Lindblad dissipator ^^ .^^^^Figure 13A fits one-point correlation functions (i.e. time-average) of the binned ^^ ^^signals ^^ ^^ and ^^ ^^ at different time bins ^^. The one-point correlation function of the[ ] [ ]^^ ^^filtered signal can be found by integrating the system of two coupled ODEs as discussed above in relation to Figure 10A, solved from time 0 to ^^ greater than any time in thesupport of the filter transfer function ^^, with initial conditions (^^^ , ^^ 1^=0 ^^=0 ) = (^^0, 0), andthe one-point correlation function is then ^^[^^^^]. In this case, for ^^ we integrate thesystem of ODEs of Figure 10A by setting rectangular transfer function ^^ =(^^ / Δ^^)^^^^ / ^ [^^Δ^^,(^^+1)Δ^^] with end time ^^ = (^^ + 1)Δ^^ and loss channel ^^ = ^^^ ^^depending on the quadrature. As discussed above, alternatively, one could directly integrate the time- refined average of the sharp signal against the filter function. Yet another alternative, as discussed above, is to approximate the binned signal with the sharp signal counterpartin the middle of the integration integral ^^[^^^ ] ≈ ^^[^^ ] = ^^^^[^^ ^^ ^^^^ℒ^ ^^^^ ^^ (^^0)] with ^^^^ = ^^Δ^^ +Δ^^ / 2 (which leads to relative systematic errors of 0.5% as compared to fitting with the formula according to Figure 10A for this particular example). In the example of Figure 13A, the fluorescence from second mode b is continuously measured for a duration of 200 ns, starting from vacuum, in: (i) a first dynamical configuration wherein only the second “buffer” mode b is driven; (ii) a second dynamical configuration wherein only the first “memory” mode is drive; and (iii) a third dynamical configuration wherein both modes are driven at the same time. To model the experimental data (which would be physically measured in practice) 250,000 stochastic trajectories were modelled for each drive configuration, for a total signal duration of 3 x 250,000 x 200 ns = 150 ms. The trajectories are sampled using the normalized first-orderRouchon method with step size ^^^^ = Δ^^ / 10 = 0.4 ns. For efficiency, all trajectories aresimulated simultaneously by stacking them in a single large array. We run the simulation on a GPU (an NVIDIA L40S card) using the dynamiqs library (https: / / github.com / dynamiqs / dynamiqs). The total runtime is approximately 2 minutes. We reiterate that this step is not needed experimentally, as the data is readily available from the experiment via steps 601-605 and 607-611 as described above in relation to Figure 6. The system of two coupled ODEs as discussed above in relation to Figure 10A is used to fit the data ^̂^^^[^^^^], which is performed using the Levenberg-Marquardt algorithm from the SciPy library. The fitted vector of parameters ^^^^^^^^converges well to the true vector of parameters ^^, as seen in the bottom row of Figure 13A wherein ^̂^^^[^^^^] and ^^^^[^^^^] overlap, even when starting from a relatively far initial guess ^^^^^^^^^^^^. Figure 13B shows a table which summarises the result which thus characterizes quantum system 1 with the fitted parameters, their fitted value (given with a factor / 2^^ where applicable), the estimate from the proposed fitting procedure, and the relativeerror. The chosen ^^2 and ^^^^ correspond to a rate ^^2 / 2^^ = 3.60 MHz and an averagememory photon number ^̅^ = 1.3. The relative error for ^^^^ is not indicated, as thisparameter is statistically degenerate and not correctly estimated in this particular example, but which could be correctly estimated by performing the method in yet further dynamical configurations. Figure 14 shows how the one-point correlation of the modelled sharp signal in the example of Figure 13 varies with each parameter, with the X / P quadratures indicated.For each row of Figure 14 the amplitude of a single parameter ^^^^ is varied (0.7^^^^, ^^^^, 1.3^^^^)respectively indicated by light, medium, and dark shades, where the specific parameter is indicated in the y-axis label). The above example is for a single cat-qubit system, however it could readily be extended to a quantum system 1 comprising a plurality of qubit-hosting structures 6 such as in Figure 4. For instance, for a two cat-qubit system the vector of parameters to fitcould be ^^ = (^^^^^^^^^^, ^^^^^^^^^^, ^^^^^^), where ^^^^^^ is a set of two-qubit parameters describingthe CNOT gate between two cat qubits as described above, and ^^^^^^^^^^ / ^^^^^^^^could each be the set of 10 parameters described above for each of the single-cat qubit system Moreover, as will be understood, further parameters may be included depending on how the Hamiltonian and the loss jump operators describing the quantum system 1 are modelled. In particular, for the cat-qubit system including an ATS as the non-linear element 7 as described above in relation to Figure 2, the ATS may be flux-biased at apoint other than at the 0 − ^^ flux working point (also known as the “saddle point”). As willbe appreciated, in the language of the common magnetic flux ^^^^and the differentialmagnetic flux ^^ , a saddle point is a pair of fluxes satisfying (^^^^ , ^^Δ) = in of the flux quantum. For instance, flux threaded through one of the loops may be non-zero and different to ^^ flux threaded through the other of the loops, in units of the flux quantum. By preparing at least one of the dynamical configurations of quantum system 1 whereinthe ATS is flux biased at a point other than the 0 − ^^ flux working point, the Hamiltonianno longer has its “sin-sin” form as described above, such that other terms enter the modelled Hamiltonian (and therefore other physical parameters to be fitted) which areotherwise particularly small and / or insensitive at the 0 − ^^ flux working point. The presentInventors have recognised that this may lead to a richer probing of the parameter space of quantum system 1, and thus to a more accurate characterization.

Claims

Claims 1. A method of characterizing a quantum system, wherein the quantum system comprises one or more qubit-hosting structure(s), each of the one or more qubit-hosting structure(s) comprising: (I) at least one portion configured to have at least one physical oscillatory mode, (II) one or more signal generator(s) coupled to the at least one portion and configured to input control signals to physically prepare the quantum system in a dynamical configuration, and (III) a detector coupled to the at least one resonant portion and configured to measure continuous signals output from the at least one oscillatory mode to provide measurement signals from the quantum system; the method comprising: (i) physically preparing, with the one or more signal generator(s), the quantum system in a first dynamical configuration, wherein the first dynamical configuration is describable via a first stochastic master equation having a first plurality of physical parameters; (ii) physically measuring, with the detector, a continuous signal from the quantum system in the first dynamical configuration to provide a first measured continuous signal; (iii) repeating steps (i) and (ii) a plurality of times to provide a plurality of first measured continuous signals, and calculating one or more first correlation function(s) from the plurality of first measured continuous signals; (iv) physically preparing, with the one or more signal generator(s), the quantum system in a second dynamical configuration which is different to the first physical configuration, wherein the second dynamical configuration is describable via a second stochastic master equation having a second plurality of physical parameters, and wherein the quantum system is physically prepared in the second dynamical configuration such that a subset of the second plurality of physical parameters is a shared subset with a subset of the first plurality of physical parameters; (v) physically measuring, with the detector, a continuous signal from the quantum system in the second dynamical configuration to provide a second measured continuous signal; (vi) repeating steps (iv) and (v) a plurality of times to provide a plurality of second measured continuous signals, and calculating one or more second correlation function(s) from the plurality of second measured continuous signals;(vii) assigning values to the shared subset of physical parameters which characterizes the quantum system by jointly numerically fitting: (a) a first simulation of the one or more first correlation function(s), numerically simulated from the first stochastic master equation, to the calculated one or more first correlation function(s), and (b) a second simulation of the one or more second correlation function(s), numerically simulated from the second stochastic master equation, to the one or more second correlation function(s).

2. The method of claim 1, the method comprising: physically preparing, with the one or more signal generator(s), the quantum system in a third or further different dynamical configurations, wherein each of the third or further different dynamical configurations is describable via a third or further stochastic master equation having a third or further plurality of physical parameters, and wherein the quantum system is physically prepared in the third or further different dynamical configuration such that a subset of the third or further second plurality of physical parameters is part of the shared subset of physical parameters; physically measuring, with the detector, a continuous signal from the quantum system in the third or further different dynamical configurations to provide a third or further measured continuous signal; repeating said steps of physically preparing and physically measuring the third or further different dynamical configurations a plurality of times to provide a plurality of the third or further measured continuous signals, and calculating one or more third or further correlation function(s) from the plurality of the third or further measured continuous signals; and wherein assigning values to the shared subset of physical parameters which characterizes the quantum system further comprises numerically fitting a third or further simulation, numerically simulated from the third or further stochastic master equation, to the one or more third or further correlation function(s).

3. The method of any preceding claim, wherein the one or more correlation function(s) is comprised of: one-point correlation functions, two-point correlation functions, and / or three-point correlation functions.

4. The method of any preceding claim, wherein at least one of the physically prepared dynamical configurations is a steady state of the quantum system, and whereinrepeating said steps of physically preparing and physically measuring the quantum system in the steady state a plurality of times comprises: physically maintaining, with the one or more signal generator(s), the quantum system in the steady state for a single trajectory; physically measuring, with the detector, a continuous signal from the quantum system in the steady state to provide a steady state measured continuous signal; and dividing the steady state measured continuous signal into a plurality of sub- portions, and calculating the corresponding one or more correlation function(s) from the plurality of sub-portions.

5. The method of any preceding claim, wherein at least one of the physically prepared dynamical configurations is a steady state of the quantum system such that the stochastic master equation describing said at least one of the dynamical configurations has a time-constant Liouvillian ℒ; the method comprising: defining a plurality of ^^ time-bins of width Δ^^ from the measured continuous signal from said at least one of the dynamical configurations, wherein the width Δ^^ is a digitisation time of the acquisition chain of the detector; wherein numerically fitting the simulations of the one or more correlation function(s) of the steady state, the one or more correlation function(s) being ^^-point correlation function(s), comprises calculating the trace of an operator obtained by numerically evolving a numerically simulated initial state of the quantum system ^^0with a Lindblad evolution interspersed with the application of a correlation superoperator at the center ^^^′^ of each ^^-th time-bin from ^^1′to ^^^′^:wherein ^^ is the jump operator of the stochastic master equation describing said at least one of the dynamical configurations, ^^′^^ = ^^Δ^^ + Δ^^ / 2, and ^^^^(^^) is thecorrelation superoperator for an operator ^^ and defined either by: (i) ^^^^(^^) = ^^^^ + ^^^^^^^^†if the detector has binary output with dark count rate ^^ ≥ 0 and detector efficiency and 0 < ^^ ≤ 1, such that the stochasticmaster equation describing the steady state is of the jump type; orthe detector has continuous output with detector efficiency and 0 < ^^ ≤ 1, such that the stochastic master equationdescribing the steady state is of the diffuse type.

6. The method of any preceding claim, wherein at least one of the physically prepared dynamical configurations is a steady state of the quantum system such that the stochastic master equation describing said at least one of the dynamical configurations has a time-constant Liouvillian ℒ; the method comprises: defining a plurality of ^^ time-points from the measured continuous signal from said at least one of the dynamical configurations; wherein numerically fitting the simulations of the one or more correlation function(s) of the steady state, the one or more correlation function(s) being ^^-pointcorrelation function(s), comprises computing the ^^-dimensional integral of ^^[^^^^1 …against a transfer function ^^^^(^^^^) of the acquisition chain of the detector for the ^^-thtime-point,whereinis obtained by calculating the trace of an operator obtained by numerically evolving a numerically simulated initial state of the quantum system ^^0with a Lindblad evolution interspersed with the application of a correlation superoperator at each time-point fromwherein ^^ is the jump operator of the stochastic master equation describing said at least one of the dynamical configurations, and ^^^^(^^)is the correlation superoperator for an operator ^^ and defined either by: (i) ^^^^(^^) = ^^^^ + ^^^^^^^^†if the detector has binary output with dark count rate ^^ ≥ 0 and detector efficiency and 0 < ^^ ≤ 1, such that the stochasticmaster equation describing the steady state is of the jump type; orthe detector has continuous output with detector efficiency and 0 < ^^ ≤ 1, such that the stochastic master equationdescribing the steady state is of the diffuse type.

7. The method of any preceding claim, wherein numerically fitting at least one of the simulations of the one or more correlation function(s) of a given one of the dynamical configurations, the one or more correlation function(s) being ^^-point correlation function(s), comprises: defining a plurality of ^^ time-points from the corresponding measured continuous signal from the given dynamical configuration; assigning a transfer function ^^^^of the acquisition chain of the detector at time- point ^^;defining ^^ =+ ⋯ + α^^^^^^ with {^^1, … , ^^^^} real numbers;deriving a system of 2^^coupled ordinary differential equations by forward differentiating an ordinary differential equation ^^^^ ^^ / ^^^^ = ℒ^^ ^^ ^^such that the system of 2^^coupled ordinary differential equations is describable by:wherein ^^ ^^ is a generating density matrix at time-point ^^, and ℒ ^^ ^^ ^^is a superoperator defined either by: (i) ℒ ^^= ℒ + (^^ ^^^^ ^^ ^^ − 1)^^^^ if the stochastic master equation describing thegiven dynamical configuration is of the jump type, ^^^^(^^)is a correlation superoperator for an operator ^^ defined ^^^^(^^) = ^^^^ + ^^^^^^^^†if the detector has binary output with dark count rate ^^ ≥ 0 and detector efficiency and 0 < ^^ ≤ 1,(ii) ℒ ^^= ℒ + ^^ ^^ + ^^2 / 2 ℐ if the stochastic master equation describing^^ ^^ ^^ ^^^^the given dynamical configuration is of the diffusive type, with ℐ the identify ( ) operator, and ^^ ^^ is a correlation superoperator for an operator ^^ defined^^if the detector has continuous output with detectorefficiency and 0 < ^^ ≤ 1;^^ solving the system of 2 coupled ordinary differential equations from time ^^ = 0to time ^^, where ^^ is chosen to be greater than any time in the support of ^^ , … , ^^ , to1 ^^1,…,^^thereby compute ^^ ; and^^1,…,^^calculating the trace of ^^ to provide the simulations of the one or more^^correlation function(s)8. The method of claim 7, wherein transfer function ^^ of the acquisition chain of^^the detector at time-point ^^ is approximated by a rectangular function of width Δ^^, wherein the width Δ^^ is a digitisation time of the acquisition chain of the detector.

9. The method of any preceding claim, wherein, for each of the one or more qubit- hosting structure(s), the at least one portion is a resonant portion configured to have at least one resonant physical oscillatory mode, each of the one or more qubit-hosting structure(s) further comprises a non-linear element coupled to the at least one resonantportion so as to non-linearly couple to the at least one resonant physical oscillatory mode, and wherein the one or more signal generator(s) are coupled to the at least one resonant portion and / or to the non-linear element, and the detector is coupled to the at least one resonant portion or to the non-linear element and configured to measure continuous signals output from the at least one resonant physical oscillatory mode to provide measurement signals from the quantum system.

10. The method of claim 9, wherein quantum system is for hosting cat qubit; wherein at least one of the one or more qubit-hosting structure(s) is a cat- qubit-hosting structure; wherein the at least one resonant portion is configured to have a first resonant physical oscillatory mode and a second resonant physical oscillatory mode; wherein the non-linear element coupled to the at least one resonant portion is configured to non-linearly couple the first resonant physical oscillatory mode and the second resonant physical oscillatory mode such that, under input control signals to the non-linear element, a two-to-one boson exchange is engineered between the first resonant physical oscillatory mode and the second resonant physical oscillatory mode; wherein a signal generator of the one or more signal generator(s) is coupled to the at least one resonant portion, said signal generator being configured to drive the second resonant physical oscillatory mode; and wherein the detector is coupled to the at least one resonant portion and configured to measure continuous signals output from the second resonant physical oscillatory mode.

11. The method of claim 10, wherein physically preparing, with the one or more signal generator(s), the quantum system in the first dynamical configuration comprises physically preparing a cat qubit hosted in the first resonant physical oscillatory mode; wherein step (iii) comprises calculating only one or more non-null correlation function(s) from the first dynamical configuration which are even-order correlation functions; and wherein, step (vii), numerically fitting the first simulation of the one or more first correlation function(s) comprises numerically simulating only one or more even-order correlation function(s).

12. The method of claims 10 or 11, wherein physically preparing different ones of the dynamical configurations comprises inputting, with the one or more signal generator(s), respectively different control signals to the quantum system, wherein the different control signals are configured to drive the first resonant physical oscillatory mode and / or second resonant physical oscillatory mode.

13. The method of any of claims 10-12, wherein the non-linear element is an asymmetrically threaded SQUID (ATS) comprised of a first superconducting loop and a second superconducting loop, wherein a signal generator of the one or more signal generator(s) is coupled to the ATS and configured to thread the first and second superconducting loops together with a common magnetic flux ^^^^and a differential magnetic flux ^^Δ; and wherein physically preparing different ones of the dynamical configurations comprises inputting, with said signal generator, respectively different pairs of common ^^^^and differential ^^Δmagnetic fluxes (^^^^, ^^Δ), and wherein at least one of the pairs isnot at a saddle point14. The method of any preceding claim, wherein the quantum system comprises a plurality of qubit-hosting structures, wherein each of the qubit-hosting structures is respectively coupled to at least one other of the qubit-hosting structures by a multi-qubit linear coupling element or by a multi-qubit non-linear coupling element; wherein the method comprises: in at least one of physically prepared dynamical configurations, physically performing, with the one or more signal generator(s), a multi-qubit gate operation between a qubit hosted in one of the qubit-hosting structures and at least one other qubit hosted in at least one other of the qubit-hosting structures; wherein at least one physical parameter of the shared subset of physical parameters is a physical parameter which characterizes the multi-qubit gate operation.

15. A quantum characterization unit comprising a quantum system comprising:one or more qubit-hosting structure(s), each of the one or more qubit- hosting structure(s) comprising: (I) at least one portion configured to have at least one physical oscillatory mode, (II) one or more signal generator(s) coupled to the at least one portion and configured to input control signals to physically prepare the quantum system in a dynamical configuration, and (III) a detector coupled to the at least one resonant portion and configured to measure continuous signals output from the at least one oscillatory mode to provide measurement signals from the quantum system; and a command circuit coupled to the quantum system, the command circuit being configured to: (i) physically prepare, with the one or more signal generator(s), the quantum system in a first dynamical configuration, wherein the first dynamical configuration is describable via a first stochastic master equation having a first plurality of physical parameters; (ii) physically measure, with the detector, a continuous signal from the quantum system in the first dynamical configuration to provide a first measured continuous signal; (iii) repeat steps (i) and (ii) a plurality of times to provide a plurality of first measured continuous signals, and calculate one or more first correlation function(s) from the plurality of first measured continuous signals; (iv) physically prepare, with the one or more signal generator(s), the quantum system in a second dynamical configuration which is different to the first physical configuration, wherein the second dynamical configuration is describable via a second stochastic master equation having a second plurality of physical parameters, and wherein the quantum system is physically prepared in the second dynamical configuration such that a subset of the second plurality of physical parameters is a shared subset with a subset of the first plurality of physical parameters; (v) physically measure, with the detector, a continuous signal from the quantum system in the second dynamical configuration to provide a second measured continuous signal;(vi) repeat steps (iv) and (v) a plurality of times to provide a plurality of second measured continuous signals, and calculating one or more second correlation function(s) from the plurality of second measured continuous signals; (vii) assign values to the shared subset of physical parameters which characterizes the quantum system by jointly numerically fitting: (a) a first simulation of the one or more first correlation function(s), numerically simulated from the first stochastic master equation, to the calculated one or more first correlation function(s), and (b) a second simulation of the one or more second correlation function(s), numerically simulated from the second stochastic master equation, to the one or more second correlation function(s).

Images