Method and unit for characterizing or controlling a quantum system

Abstract

A method of characterizing or controlling a quantum system (1), wherein the dynamics of the quantum system is describable via a stochastic master equation having a plurality of physical parameters, and wherein the quantum system comprises: (i) one or more quantum- state-hosting structure(s) (6); (ii) a detector (8) coupled to at least one of the quantum-state- hosting structure(s) (6) and configured to measure a continuous signal output from the quantum-state-hosting structure(s) (6) to provide a measurement signal from the quantum system (1), wherein the detector (8) is configured to digitize the continuous signal over time bins of duration Δt; and (iii) one or more signal generator(s) (11,13) coupled to the at least one portion and configured to input control signal(s) to physically control the dynamics of the quantum system. The method comprises: (I) physically measuring, with the detector, a continuous signal from the quantum system to provide a plurality of digitized discrete-time binned signals {I 1,...,I k}, …, up to time t = kΔt, wherein k is an integer denoting the k-th time bin of the detector; (II) computing a binned state of the quantum system at time (F1), wherein the binned state (F2) is the average of a plurality of solutions of the stochastic master equation at time t = kΔt each for different possible continuous signals wherein, for each of the plurality of solutions P k Δ t , the value of the possible continuous signal integrated over a given time bin is substantially equal to the discrete-time binned signal for said given time bin; and (III) using the binned state (F2) to: assign a value to at least one of the plurality of physical parameters which characterizes the quantum system; and / or assign a value to at least one of the plurality of physical parameters which parametrizes the input control signal(s) to physically control the dynamics of the quantum system; and / or optimize a readout discriminator of a quantum state readout protocol. There is provided a quantum characterisation or control unit, a computer program product, and a computer-readable medium.

Description

[0001] METHOD AND UNIT FOR CHARACTERIZING OR CONTROLLING A QUANTUM SYSTEM

[0002] CROSS-REFERENCE TO RELATED APPLICATION

[0003] This application claims priority from and the benefit of European patent application No. EP 25305139.5, which was filed on 31 January 2025. The entire contents of this application are incorporated herein by reference.

[0004] FIELD OF THE INVENTION

[0005] The invention concerns a method of characterizing or controlling a quantum system measured continuously, and in particular using a discrete-time stochastic master equation to reconstruct the state of the continuously measured quantum system from a digitised (time-averaged) measurement record, as well as a quantum characterization or control unit for doing same. The results of the characterization or control may be used for subsequent optimal control of the quantum system.

[0006] BACKGROUND

[0007] Characterising and controlling the dynamics of a quantum system is a fundamental task in experimental quantum physics. Accurate characterisation and control is vital for successful exploitation of the quantum system, such as a quantum processor for quantum computation, or a quantum device for quantum sensing or quantum communication applications.

[0008] For instance, crafting the control pulses to be sent to manipulate the quantum system typically requires accurate characterisation of the dynamics of the quantum system in order to provide the most optimal control pulse parameters.

[0009] Moreover, successful operation of a quantum processor requires regular recalibration 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.

[0010] Characterising the dynamics of a quantum system requires comparing experimental data measured from the quantum system to a theoretical model of the quantum system. Typically, continuously measured quantum systems are modelled bythe theory of stochastic master equation (SME), which gives an expression for the system state ptat a given time t as a non-trivial function of the continuously measured signal before that time.

[0011] However, in practice, this state itself is not accessible to the observer, as the measured signal is filtered and digitised by an amplification chain of the detection apparatus. Said otherwise, the observer only has access to a discrete-time measurement record {I1, ..., IN}, wherein each discrete value Ikis defined by integrating the continuoustime signal Itagainst the impulse (e.g. filter) response function fkof the acquisition chain for the Zc-th time bin Ik= ∫fk(t)Itdt.

[0012]

[0013] Thus, there is a problem when using this experimental data (which is inherently a discrete-time signal) to reconstruct and compare with the system state pt(which is inherently calculated assuming a continuous-time signal). That is, there is a practical gap between what is observed (the digitised signal Ik), and what is predicted by theory (how the state ptevolves conditioned on the continuous-time signal It).

[0014] This disparity is exasperated when some dynamics of the quantum system varies on time-scales on the order of, or even smaller than, the duration of the time bins At of the acquisition chain. This is the case, for instance, in superconducting circuits, wherein the typical time bin of analogue-to-digital converters or photocounters used therein is on the order of a few nanoseconds, but wherein certain physical dynamics of the superconducting circuit oscillate or vary at even smaller time scales than such time bins.

[0015] This invention aims to address one or more of the aforementioned problems to provide an improved method of characterizing a quantum system.

[0016] SUMMARY

[0017] In a first aspect, there is provided a method of characterizing or controlling a quantum system, wherein the dynamics of the quantum system is describable via a stochastic master equation having a plurality of physical parameters, and wherein the quantum system comprises: (i) one or more quantum-state-hosting structure(s); (ii) a detector coupled to at least one of the quantum-state-hosting structure(s) and configured to measure a continuous signal output from the quantum-state-hosting structure(s) to provide a measurement signal from the quantum system, wherein the detector is configured to digitize the continuous signal over time bins of duration At; and (iii) one or more signal generator(s) coupled to the at least one of the quantum-state-hosting structure(s) and configured to input control signal(s) to physically control the dynamicsof the quantum system. The method comprises: (I) physically measuring, with the detector, a continuous signal from the quantum system to provide a plurality of digitized discrete-time binned signals {I1,...,Ik} up to time t = k t, wherein k is an integer denoting the / c-th time bin of the detector; (II) computing a binned state of the quantum system at time t = k t, pk= E[ρkΔt| I1,..., Ik], wherein the binned state pkis the average of a plurality of solutions ρkΔtof the stochastic master equation at time t = k t each for different possible continuous signals wherein, for each of the plurality of solutions ρkΔt, the value of the possible continuous signal integrated over a given time bin is substantially equal to the discrete-time binned signal for said given time bin; and (III) using the binned state pkto: assign a value to at least one of the plurality of physical parameters which characterizes the quantum system; and / or assign a value to at least one of the plurality of physical parameters which parametrizes the input control signal(s) to physically control the dynamics of the quantum system; and / or optimize a readout discriminator of a quantum state readout protocol.

[0018] Thus, the present inventive method neither averages over all of the possible solutions pkAtof the stochastic master equation when ignoring the plurality of digitized discrete-time binned signals (which amounts to solving the usual Lindblad master equation for the deterministic evolution of the state), nor directly solves the continuoustime stochastic master equation by approximating the plurality of digitized discrete-time binned signals as a continuous measurement record. Rather, the present invention conditions averaging over the possible solutions pkAtof the stochastic master equation on the integral of the continuous signal being what was observed. Said otherwise, the binned state is the average of the states solving the stochastic master equation for all possible realisations of the continuous signals only where the integrated signal value is the one observed for each time bin of the plurality of digitized discrete-time binned signals. That is, the average is only over states pkAtthat are consistent with the observations. This estimate is the best Bayesian estimate of the state from the observer’s point of view, and thus by using this binned state pkto characterize or optimize control of the quantum system, the present inventors have discovered a more precise and acute characterization of the physical parameters or optimization of the control.

[0019] As will be appreciated, operation (III)(a) may be achieved via numerically fitting the binned state, as a function of at least one of the plurality of physical parameters which characterizes the quantum system, to experimental data, such as experimental data comprising a quantum state tomography (e.g. Wigner tomography or Husimitomography), to thereby find the value of said at least one of the plurality of physical parameters which best fits the experimental data.

[0020] Operation (III)(b) may comprise numerically reconstructing the state of the quantum system ρtat any time from t = 0 to t = kΔt, or numerically predicting a future state of the quantum system ρt>kΔtat a future time t > k t.

[0021] Operation (III)(b) may, for instance, comprise optimizing a loss function which depends on at least one target operating attribute to be reached by the quantum system under the input control signal(s) from the one or more signal generator(s) and / or at least one target operating attribute to be reached by the one or more signal generator(s) whilst inputting control signal(s) to physically control the dynamics of the quantum system.

[0022] Each of the one or more quantum-state-hosting structure(s) may comprise a portion configured to host a physical degree of freedom, such as a physical oscillatory mode, and wherein the detector may be coupled to at least one of the quantum-state-hosting structure(s) and configured to measure continuous signals output from the physical degree of freedom to provide measurement signals from the quantum system.

[0023] In an embodiment, the detector has binary output with dark count rate 9 > 0 and detector efficiency and 0 < p < 1, such that the stochastic master equation describing the dynamics of the quantum system is of the jump type, and wherein operation (II) comprises computing the binned state pkby: (i) assigning an initial binned state p0= p0at time t = 0 describing the initial state of the quantum system; and (ii) calculating the binned state of the quantum system at time t = kAt via the iterative formula pk= wherein %kis a completely positive map and is given by %k= dpe‘p / fcexp(k^1)At(£ + (e~ip- 1)CL), wherein exp(kI1)4t(0) denotes the time-

[0024]

[0025] ordered exponential of the superoperator 0, L is the Liouvillian of the stochastic master equation describing the quantum system, and CLis the correlation superoperator defined for an operator 0 by CL(O) = θO + ηLOL†, with L being the measured jump operator of the stochastic master equation describing the dynamics of the quantum system.

[0026] The action of the completely positive map %kon the binned state p̄k-1at the k-th time bin, Kk(p̄k-1), may be computed using a numerical quadrature, and preferably using the Gauss-Legendre quadrature.

[0027] Other numerical quadratures may be used, such as other Gaussian quadratures. In an embodiment, the detector has an output taking a continuous range of values with detector efficiency and 0 < p < 1, such that the stochastic master equation describing the dynamics of the quantum system is of the diffusive type, and whereinoperation (II) comprises computing the binned state pkby: (i) assigning an initial binned state p0= p0at time t = 0 describing the initial state of the quantum system; and (ii) calculating the binned state of the quantum system at time t = kAt via the iterative formula pk= wherein %kis a completely positive map and is given by %k=

[0028]

[0029] (1 / 2π)∫dpeipI_k - (Δt / 2)p²expkΔt(k-1)Δt(L - ipCL), wherein expkΔt(k-1)Δt(O) denotes the time-ordered exponential of the superoperator O, L is the Liouvillian of the stochastic master equation describing the quantum system, and CLis the correlation superoperator defined for an operator 0 by CL(O) = √η(LO + OL†), with L being the measured jump operator of the stochastic master equation describing the dynamics of the quantum system.

[0030] The action of the completely positive map %kon the binned state pk-at the k- th time bin, Kk(p̄k-1), may be computed using a numerical quadrature, and preferably using the Gauss-Hermite quadrature.

[0031] Other numerical quadratures may be used, such as other Gaussian quadratures. It will be appreciated that the form of the completely positive map %kfor quantum systems described by stochastic master equations of both the jump and diffuse types can be elegantly expressed as Kk

[0032]

[0033] = dpe‘p / fcexp(kI1)4t(£~ipik), where £jis a generating Liouvillian defined for a function j = -ip^kas: (i) Ljt= Lt+ (ej_t- 1)CLif the stochastic master equation describing the dynamics of the quantum system is of the jump type, where CLis the correlation superoperator defined for an operator 0 by CL(O) = θO + ηLOL†if the detector has binary output with dark count rate 9 > 0 and detector efficiency and 0 < p < 1; or (ii) Ltj= Lt+ jtCL+ jt2 / 2 핀 if the stochastic master equation describing dynamics of the quantum system is of the diffusive type, with J the identify superoperator, and CLis a correlation superoperator defined for an operator 0 by CL(O) = √η(LO + OL†) if the detector has continuous output with detector efficiency

[0034]

[0035] and 0 < p < 1. In both cases, L is the measured jump operator of the stochastic master equation describing the dynamics of the quantum system, and £tis the standard Liouvillian of the stochastic master equation describing the dynamics of the quantum system.

[0036] For the avoidance of doubt, the time-ordered exponential

[0037]

[0038] expkΔt(k-1)Δt(O) of a superoperator O is between the earlier time (k - 1)Δt and the later time kΔt > (k - 1)Δt.

[0039] In an embodiment, operation (II) comprises computing the binned state of the quantum system at time t = k t, pk, perturbatively.The detector may have an output taking a continuous range of values with detector efficiency and 0 < p < 1, such that the stochastic master equation describing the dynamics of the quantum system is of the diffusive type, and wherein computing the binned state of the quantum system perturbatively may comprise: (i) assigning an initial binned state p0= p0at time t = 0 describing the initial state of the quantum system; (ii) calculating the binned state of the quantum system at time t = kAt via the iterative

[0040]

[0041] formula p̄k= Kk(p̄k-1) / Tr[Kk(p̄k-1)], wherein Kkis a completely positive map and is given by using: (a) Kk= (e-ik / 2 / √(2πΔt))[풥 + √Δt ikCL]; or (b) Kk= (e-ik / 2 / √(2πΔt))[풥 + √Δt1ikCL+ √Δt2(L - (1-ik²) / 2 CL²)]; or (c) Kk= (e-ik / 2 / √(2πΔt))[풥 + √Δt ikCL+ √Δt2(L - (1-ik²) / 2 CL²) + ...]; or (d) Kk= ..., wherein £ is

[0042]

[0043] the Liouvillian of the stochastic master equation describing the quantum system, and CLis the correlation superoperator defined for an operator 0 by CL(O

[0044]

[0045] ) = y / p(LO + OL1), with L being the measured jump operator of the stochastic master equation describing the dynamics of the quantum system, and ik= Ik / √Δt.

[0046] In an embodiment, the quantum system comprises a plurality of detectors each coupled to at least one of the quantum-state-hosting structure(s) and configured to measure continuous signals output from the quantum-state-hosting structure(s) to provide a plurality of measurement signal from the quantum system, wherein the detector is configured to digitize the continuous signal over time bins of duration At, the method comprising: in operation (I), physically measuring, with the detectors, a plurality of continuous signals from the quantum system to provide a plurality N of sets of digitized discrete-time binned signals {In l,..., In>k} up to time t = / cAt, wherein k is an integer denoting the Zc-th time bin of each detector, and wherein index n e [1; / V] identifies each detector; and in operation (II) computing the binned state pkby: (i) assigning an initial binned state p0= p0at time t = 0; and (ii) calculating the binned state of the quantum system at time t = k t via the iterative formula pk= wherein %kis a

[0047]

[0048] completely positive map and is given by Kk= (1 / 2π)N∫...∫dp1...dpN(∏n=1NeipI_{n,k}) expkΔt(k-1)Δt(Lj), wherein ∫...∫dp1...dpNdenotes the multiple integral across variable pn, wherein

[0049]

[0050] expkΔt(k-1)Δt(O) denotes the time-orderedexponential of the superoperator 0, Ljis a generating Liouvillian superoperator depending on the binary or continuous outputs of each of the detectors.

[0051] In an embodiment, wherein the dynamics of the quantum system is configured to be in a steady state, the method may comprise: in operation (I) physically measuring, with the detector, a plurality of sequential continuous signals from the quantum system to provide a plurality N of sets of digitized discrete-time binned signals {In,1,...,In,k} up to time t = k t, wherein k is an integer denoting the Zc-th time bin of each set, and wherein index n e [1; N] identifies each set; and in operation (II) computing the binned state pkby: (i) assigning an initial binned state p0= p0at time t = 0; and (ii) calculating the binned state of the quantum system at time t = k t via the iterative formula pk= wherein %kis a completely positive map and is given by TrpCfctpfc.i)]

[0052]

[0053] k ~

[0054]

[0055] (1 / 2π)N∫...∫dp1...dpN(∏n=1NeipI_{n,k}) expkΔt(k-1)Δt(Lj), wherein ∫...∫dp1...dpNdenotes the multiple integral across variable pn, wherein

[0056]

[0057] expkΔt(k-1)Δt(O) denotes the time-ordered exponential of the superoperator 0, Ljis a generating Liouvillian superoperator depending on the binary or continuous output of the detector.

[0058] For the avoidance of doubt, herein calculating a quantity via an iterative formula comprises (i) receiving an initial state; (ii) defining a computational rule expressed as a mathematical formula which defines the relationship between consecutive elements in a sequence starting from the initial state, the formula being configured to update a subsequent state based on a preceding state; (iii) executing (e.g. using a conventional computer) a first application of the computational rule to the initial state; (iv) iteratively or recursively applying the computational rule to the result of the preceding application, such that each subsequent computation is based, at least in part, on one or more previously computed states, thereby generating a sequence of computed output states; and (v) terminating further applications of the computational rule after a threshold (e.g. pre-determined maximum) number of iterations. Thus, in this context, the initial state may be set as p0= p0at time t = 0. Then, the map

[0059]

[0060] may be applied to p̄0, such that the next state at the k = 1 time bin is thus calculated via p̄1= K1(p̄0) / Tr[K1(p̄0)]. The map K2may

[0061]

[0062] then be applied to p, such that the next state p2is calculated via p2= and so TrpCzCpi)]’ on iteratively up to the / c-th time bin to provide pk.

[0063] The action of the completely positive map %kon the binned state pk-at the k-th time bin, Kk(p̄k-1), may be computed by evaluating the integrand for a plurality ofvalues of p, wherein at least some of the integrands are evaluated simultaneously using a vectorized-array computation.

[0064] That is, computing Kk(p̄k-1) exactly (to machine precision) typically requires evaluating the integrand for a few tens of values of p wherein, for each value of p, the cost is the same as solving the Lindblad master equation. Thus, by evaluating at least some of the integrands simultaneously using a vectorized-array computation, the computation is more efficient in time. Preferably, all of the integrands for the plurality of values of p are computed simultaneously, thereby reducing the cost of the iterating scheme to solving a single Lindblad master equation for the same system.

[0065] In embodiments, wherein optionally: (A) the plurality of physical parameters comprises a first set {θ1} of one or more physical parameter(s) in the Liouvillian L(θ1) and / or the measurement back-action

[0066]

[0067] of the stochastic master equation describing the quantum system and / or a second set {02} of one or more physical parameters of the initial state p0(e2) at time t = 0 describing the initial state of the quantum system such that the binned state pkis a function of the first {0i} and / or second {02} sets, p̄k(θ1,2), and wherein operation (I I l)(a) comprises: numerically fitting the binned state p̄k(θ1,2) to the plurality of digitized discrete-time binned signals {I1,...,Ik} to assign a value to at least one physical parameter of the first {θ1} and / or second {e2} sets; and / or (B) the plurality of physical parameters comprises a third set of one or more physical parameters {03} which parametrizes the input control signal(s) of the one or signal generator(s) such that the binned state pkis a function of the third set {e3}, p̄k(θ3), and wherein operation (lll)(b) comprises: defining a loss function as a function of the fidelity between a target quantum state ρtargetand the binned state pp̄k(θ3) at time t = kAt, F(ρtargetkΔt, p̄k(θ)); determining a local minimum {0min) of the loss function within the space of physical control parameters of the input control signals; and using the local minimum {0min) as the set of physical parameters of the input control signal(s) in a subsequent control operation of the quantum system.

[0068] Operation (III) may comprise calculating the probability of the discrete-time binned signal for the Zc-th time bin using the binned state pk_! at the ( / c - l)-th time bin using the formula P[Ik| p̄k-1] = Tr[Kk(p̄k-1)].

[0069] As will be appreciated, the probability P[Ik| p̄k-1] of the next digitized discretetime signal to be measured Ikmay be used to determine the most optimal subsequent time bin for the next acquisition which maximizes the information gained from the experiment, or may be used to optimize a loss function and thus determine the physicalcontrol parameters of the signal generator(s) to be used which should achieve the optimized loss function.

[0070] There is also provided, in a second aspect, a quantum characterisation or control unit comprising: a quantum system which comprises: (i) one or more quantum-state-hosting structure(s); (ii) a detector coupled to at least one of the quantum-state-hosting structure(s) and configured to measure continuous signals output from the one or more quantum-state-hosting structure(s) to provide measurement signals from the quantum system, wherein the detector uses time bins of duration At; and (iii) one or more signal generator(s) coupled to one or more quantum-state-hosting structure(s) and configured to input control signals to physically control the dynamics of the quantum system. The command circuit is configured to: (a) control the one or more signal generator(s) so as to physically control the dynamics of the quantum system, and (b) control the detector so as to provide measurement signals from the quantum system; wherein the command circuit is configured to perform the method according to the first aspect and any of its embodiments.

[0071] For the avoidance of doubt, the time bin At is time duration, which is the digitisation time of the acquisition chain of the detector. Formally, the digitized discretetime binned signals {I1,..., Ik} can be defined by integrating the continuous-time signal It(a theoretical quantity not accessible to the observer) against the impulse response fk(t) of the acquisition chain of the for the / c-th time bin Ik= ∫fk(t)Itdt. The

[0072]

[0073] transfer function fk(t) of the acquisition chain of the detector at time-point t = k t may be approximated by a rectangular function of width At, such that the dynamics of the binned quantum state is Markovian. Of course, other transfer functions fk(t) may be known, depending on the specific form of the detector, and preferably such other transfer functions do not contribute to non-Markovian dynamics of the quantum system.

[0074] However, the present invention can also handle transfer functions that contribute to non-Markovian dynamics. For instance, the generalisation to condition on multiple signals can handle filtering functions that overlap over different time-bins and introduce some superficial memory effect.

[0075] The present invention is particularly applicable to characterizing or controlling non-linear open quantum systems. Examples of such include superconducting circuits which may comprise transmon qubits, fluxonium qubits, GKP qubits, etc.

[0076] That is, in embodiments of the first and second aspects, for each of the one or more quantum-state-hosting structure(s) is a qubit-hosting structure(s), the at least one portion is a resonant portion configured to have at least one resonant physical oscillatorymode, 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 nonlinear 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.

[0077] The present invention is furthermore particularly applicable to so-called cat qubits. Thus, in the first and second aspects, the at least one of the quantum-state-hosting structure(s) may be a cat-qubit-hosting structure wherein the physical oscillatory mode hosts a cat qubit. In particular, the cat-qubit-hosting structure may be realized in a superconducting circuit architecture.

[0078] Preferably, the command circuit is configured to control the one or more signal generator(s) so as to physically stabilize and manipulate a dissipative cat qubit.

[0079] Thus, in some embodiments the quantum system is for hosting cat qubit, the at least one resonant portion may be configured to have a first resonant physical oscillatory mode and a second resonant physical oscillatory mode. The non-linear element coupled to the at least one resonant portion may be 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. The signal generator of the one or more signal generator(s) may be coupled to the at least one resonant portion, said signal generator being configured to drive the second resonant physical oscillatory mode. The detector may be coupled to the at least one resonant portion and configured to measure continuous signals output from the second resonant physical oscillatory mode.

[0080] Specifically, stabilizing a cat qubit may comprise physically preparing a cat qubit hosted in the first resonant physical oscillatory mode.

[0081] The quantum characterisation or control unit may be part of a quantum processor. The second aspect may comprise any of the features and embodiments described above in relation to the first aspect.

[0082] There is also provided a computer program product comprising instructions which, when executed by the quantum characterisation or control unit of the secondaspect and any of its embodiments, cause the quantum characterisation or control unit to perform the method of the first aspect and any of its embodiments.

[0083] There is also provided a computer-readable medium having stored thereon the computer program product.

[0084] The quantum characterisation or control unit may be part of a quantum processor.

[0085] BRIEF DESCRIPTION OF THE DRAWINGS

[0086] Various embodiments of the invention will now be described, by way of example only, and with reference to the accompanying drawings in which:

[0087] Figure 1 shows a diagram of a generic quantum system to be characterized according to an embodiment of the invention;

[0088] 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;

[0089] Figure 3 shows a flow diagram of method steps for characterizing or controlling a quantum system, according to an embodiment of the invention;

[0090] Figure 4 shows how the filtering and digitisation of a time-continuous diffusive signal is modelled, according to an embodiment of the invention;

[0091] Figure 5 shows a method for extracting physical parameters characterizing a quantum system, according to an embodiment of the invention;

[0092] Figure 6 shows a method of controlling a quantum system, according to an embodiment;

[0093] Figures 7-9 show numerically simulated results of the method according to an embodiment of the present invention applied to a discrete-level qubit system; and Figures 10 and 11 show numerically simulated results of the method according to an embodiment of the present invention applied to a dissipative cat qubit system.

[0094] DETAILED DESCRIPTION

[0095] Quantum system to be characterized and / or controlled

[0096] Figure 1 shows an exemplary quantum system 1 which may be arranged to stabilize or host one or more qubits. The quantum system 1 therefore comprises one or more qubit-hosting structures 6, each of which is configured to host or stabilize a qubittherein. 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.

[0097] 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 a and a second physical oscillatory mode b. 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.

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

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

[0100] Preparing the quantum system 1 in a dynamical configuration (i.e. setting an initial state of the quantum system 1) 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.

[0101] 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 discussedmore 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.

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

[0103] 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 nonlinear couplers 4 where appropriate, as indicated via the dashed lines in Figure 1.

[0104] 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 or control unit. The quantum characterization or control unit may itself be part of a quantum processor, a quantum communication system, or a quantum detection or sensing system, as the case may be according to the specific technical use the quantum system is put to.

[0105] Thus it will be appreciated that the quantum system 1 may comprise one or more quantum-state-hosting structures 6 (as opposed for specifically qubit-hosting structure), each of which is configured to host or stabilize a quantum state for quantum computation, quantum communication or quantum sensing (e.g. qutrits, qudits, quantum harmonic oscillators, quantum registers, etc.), or indeed any other industrial application exploiting the quantum properties of the quantum system 1. For instance, each quantum-state-hosting structure may comprise a portion configured to host a physical degree of freedom (which may or may not be oscillatory in the typical sense - such as electrons in a metal).The method of the present invention may be applied to a host of different quantum systems, provided a time-continuous signal output from the quantum system is observable by a detector.

[0106] For instance, the quantum system could be a driven anharmonic oscillator under heterodyne detection. Such a quantum system could for instance be realised by a transmon coupled to a superconducting resonator with one or more electromagnetic drives applied to the transmon and / or the superconducting resonator.

[0107] 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 a tapped-ion qubit (with the computational states being specific electronic states of the ion) placed inside a high-quality and high-finesse optical cavity such that the ion interacts with a quantized electromagnetic field mode of the cavity’s field, wherein a laser is used to drive the ion and a photodetector outside configured to detect spontaneous emission from the ion.

[0108] The present Inventors have discovered that the present invention is particularly suited to characterising and controlling so-called cat-qubits, which typically face difficulties in swiftly and accurately extracting physical parameters which characterize the quantum system, and which moreover often exhibit dynamics on the order of or on shorter time periods than the shortest time bins possible for ADC detector systems (e.g less than or on the order of a few nanoseconds). A two-legged cat qubit is defined as a two-dimensional manifold spanned by the so-called cat states \

[0109]

[0110] C±) which are superpositions of two coherent states |

[0111]

[0112] a) and |— a>: \C±) = N±(\a) ± |-a)), where N±= 1 / ^2(1 ± e-2l“l2).

[0113]

[0114] 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 an exponential way with the “size” - i.e. the average number of photons n̄ = |α|2- of the Schrodinger 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.

[0115] The cat qubits can be stabilized or confined by the following exemplary schemes: (A) A parametric dissipative stabilization, with jump operator L2= √κ2(a2- α2), where K2is the two-photon dissipation rate, a is the photon annihilation operator of the memory mode a and a is a complex number defining the cat qubit. This jump operator can be realized by coupling a lossy buffer mode b with dissipation rate Kb, 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

[0116]

[0117] =tg2(a2_a2) / j‘1' + h. c., where b is the photon annihilation operator of the buffer mode b and g2is the two-photon coupling rate, by applying to the four-wave mixing nonlinear element a pump at frequency \2fa-fb\ and a drive of the buffer mode b at frequency fbprovided g2< Kb.

[0118] (B) A Kerr Hamiltonian = y (at2- a2) (a2- a2^, where K is the amplitude

[0119]

[0120] of the Kerr Hamiltonian, a is the photon annihilation operator, and |a|2is the mean photon number.

[0121] (C) A detuned Kerr Hamiltonian (at2- a2) (a2- a2^ - Aa^a, where K

[0122]

[0123] is the amplitude of the Kerr Hamiltonian, a is the photon annihilation operator, a is a complex number defining the cat qubit, and A is the detuning factor.

[0124] (D) A two-photon exchange (TPE) HamiltonianH / fl= g2(ci2~ a2)o’+ + h.c., where g2is the complex two-photon coupling rate, a is the photon annihilation operator, a 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).

[0125] (E) A dissipative squeezing stabilization, with jump operator Lsc= (^(cosh(r)a + sinh(r)e‘0a‘1)2- a2^), where KSCis the squeezed two-photon dissipation rate, a is the photon annihilation operator of the memory mode a, a is a complex number defining the cat qubit, r and 9 are the modulus and argument of the complex squeezing parameter = reie. This jump operator can be realized by coupling a lossy buffer mode b with dissipation rate Kb, and a four-wave mixing non-linear element - typically a Josephson junction or an ATS -, to the memory mode a and by engineering

[0126]

[0127] the Hamiltonian = gsc(cosh(r)a + sinh(r)e‘6'cit)2- a2) / / 1' + h. c., where b is the

[0128]

[0129] photon annihilation operator of the buffer mode b and gscis the squeezed two-photon coupling rate, with several pumps at frequencies \2fa- fb\, fband 2fa+ fb, and a drive of the buffer mode b at frequency fbprovided gsc< Kb.

[0130] (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 - is stabilized by engineering the Hamiltonian = g2(a2+ cta - a2)bt+ h. c., where g2is the amplitude of the two-photon coupling rate as described above achieved via a pumpat frequency \2fa- fb\, a is the annihilation operator of the memory mode a, A is a complex number which phase and amplitude result from the amplitude of longitudinal coupling produced by a pump at frequency fb, a is a complex number resulting from a drive of the buffer mode b at frequency fband b is the annihilation operator of the buffer mode b; a comparison between the moon cat qubit and the squeezed cat qubit could be established by expressing A as a function of the complex squeezing parameter = reieas follows: A = 2tanh (r).

[0131] (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 appropriate voltage bias, provides the required parametric interaction at the frequency \2fa- fb\ necessary to achieve dissipative stabilization. In particular, the two-photon coupling rate g2is therefore not limited by the amplitude of the two-photon pump: g2= (EJ / 4)φa2φb, where Ej is the Josephson energy of the one or more Josephson junctions, <pais the zero-point fluctuation of the phase of the memory mode a, and <pbis the zero-point fluctuation of the phase of the buffer mode b.

[0132] (H) A resonant dissipative stabilization, for which the Applicant filed the European patent application EP 21306965.1, with jump operator L2= √κ2(a2- α2), where K2is the two-photon dissipation rate, a is the photon annihilation operator of the memory mode a and a is a complex number defining the cat qubit. This jump operator can be realized by coupling a lossy buffer mode b with dissipation rate Kb, and a three-wave mixing nonlinear element to the memory mode a which engineers the Hamiltonian

[0133]

[0134] = g2(a2-α2)b†+ h. c., where b is the photon annihilation operator of the buffer mode b provided the mode frequencies verify substantially 2fa= fband g2< Kbto which a drive of the buffer mode at frequency fbis added.

[0135] 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 arrangedto make possible three- or four-wave mixing between a first physical oscillatory mode a and a second physical oscillatory mode b.

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

[0137] 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 fa=

[0138]

[0139]

[0140] and the second mode b has a resonant frequency fb=

[0141]

[0142]

[0143] where ωaand ωbare the respective angular frequencies of the first mode a and the second mode b.

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

[0145] The memory mode a has a high-quality factor Qawhile the buffer mode b has a low-quality factor Qb. 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 Q = / / A, wherein f is resonant frequency and A is the spectroscopic linewidth; or (b) via a time-domain measurement wherein a tone is sent in, and a return signal is measured after a pre-determined time, wherein Q = f * T with T 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.

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

[0147] 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 forparametric 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.

[0148] 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 faand fb. 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 nonlinear element 7. Such a participation can be quantified by the zero-point fluctuation of the superconducting phase across the non-linear element 7, noted <pafor the first mode a and <pbfor the second mode b.

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

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

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

[0152] 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), similarto 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 pm3 of a suspended nanostructure.

[0153] In such acoustic resonators, the cat qubit is encoded in the phonons of the resonator.

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

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

[0156] 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).

[0157] 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 arrangedat least to drive the second mode b by delivering radiation at a frequency substantially equal to the second resonant frequency fbto the resonant portion 9.

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

[0159] 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 faand fb, 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.

[0160] 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

[0161]

[0162] crosses and label in Figure 2) being shunted by an inductance (indicated by the ELlabel 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.

[0163] 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 a2of the jump operator

[0164]

[0165] 'norder 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 by parametrically pumping the ATS at the frequency fp= \2fa- fb\, as successfully shown by Lescanne, Raphael, et al. " Exponential suppression of bit-flips in a qubit encoded in an oscillator." Nature Physics 16.5 (2020): 509-513. Advantageously, the pump frequency satisfies fp» g2to make this parametric pumping work as well as possible.

[0166] It is known to the skilled person that, in the case of the non-linear element 7 being an ATS 7, when biased at its flux working point also known as the “saddle point” (0 - n, namely 0 flux threaded through one of the loops and n 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:

[0167] EL 2

[0168] HATS = -2£} sin(<pz(t)) sin(rp) + — (<p - <pA(t)),

[0169]

[0170] where <p = <pa(a + a1') + <pb(b + ) is the total superconducting phase difference across the ATS 7, pais the zero-point fluctuation of the phase of the first mode a across the ATS 7 and rpbis the zero-point fluctuation of the phase of the second mode b across the ATS 7, Ej is the Josephson energy of the side junctions and ELis the inductive energy of a central inductance, <pz(t) corresponds to a common flux modulation of the two loops of the ATS 7 and <pA(t) 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.

[0171] 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 an example, by pumping the common flux at the frequency fp= \2fa-fb\,

[0172]

[0173] = e2phcos (2nfpt), the non-linear resonant part of the Hamiltonian writes in the rotating frame: HATS= (E]e2php2pb / 2') *

[0174]

[0175] (a2bt+ h.c.), which is typically the two-to-one photon exchange Hamiltonian needed to engineer the two-photon stabilization.

[0176] When an external DC magnetic field is set such that a 0 mod 2n magnetic flux threads one of the loop and a n mod 2?r 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-teeconnected 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.

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

[0178] When the electromagnetic source 17 is set at the frequency \2fa- fb\, the nonlinear 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 fb.

[0179] Alternatively, the electromagnetic filter 25 may be configured as a band stop filter at a frequency faand 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.

[0180] Alternatively, it may be configured as a low-pass (respectively high-pass) filter if fa > fb (resP fb > fa)- In other embodiments, the electromagnetic filter 25 can be omitted when 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.

[0181] As previously explained, the second mode b is driven at its resonant frequency fb. This two-photon drive is performed by a electromagnetic source 13 set at frequency fb- 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.

[0182] 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 electromagneticsource 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.

[0183] The electromagnetic source 11 may be arranged to drive the first mode a by delivering electromagnetic radiation at the frequency fato 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 Hzexpressed asHz / h= eza + h. c., where the complex rate ezresults from the amplitude and phase of the drive of the first mode a.

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

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

[0186] 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 fbto the linear electromagnetic network b / a to drive the second mode b.

[0187] 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 <pzor “sigma flux”) and provide (and optionally modulate) the differential flux (herein also referred to as <p&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 fbto 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.

[0188] As will be appreciated, the qubit-hosting structure of Figure 2 is just one specific example, and variations thereof may be conceived for dissipatively stabilizing a cat qubit. For instance, the non-linear element 7, which is an ATS in Figure 2, may be replaced by a three-wave mixing non-linear element 7, which may for instance be formed by at least one loop including a first Josephson junction, a central inductive element and a second Josephson junction which are configured such that when a predetermined current of a constant intensity is applied by a current source (which is a signal generator replacing the electromagnetic sources 17,27 and the electromagnetic network 19 of Figure 2), the resonant frequency fbis substantially equal to twice the resonant frequency fa. Such a three-wave mixing non-linear element and the associated superconducting circuit (and variations thereof) are described in EP21306965.1 filed by the Applicant. Other alternatives of the qubit-hosting structure may be used, such as that comprising a DC voltage-biased Josephson junction, as described in EP 23306839.4 and Aissaoui, Thiziri, et al. " A cat qubit stabilization scheme using a voltage biased Josephson junction." arXiv preprint arXiv:2411.08132 (2024).

[0189] Characterization or control method

[0190] The present Inventors have developed a method for characterizing or controlling a quantum system 1, whether it be comprised of a single quantum-state-hosting structure or multiple quantum-state-hosting structures such as those outlined in the above embodiments, and which crucially operates as a non-linear dynamical system.

[0191] Figure 3 shows a flow diagram of the method steps for characterizing a quantum system 1 which is operable as an open quantum dynamical system, for instance as described above in relation to Figures 1-2. It is noted here that the present invention can be used to characterize or control any open quantum system which is continuously measured. That is, the quantum system 1, although typically (and indeed preferably)non-linear for many applications, need not necessarily be non-linear. Ergo, the quantum system may not comprise a non-linear element 7 in some embodiments.

[0192] 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 Itmeasured by the observer at time t.

[0193] Specifically, an observer describes the state of an open quantum system 1 by a mixed quantum state, the density matrix ptat time t. Mathematically, we associate with the system a Hilbert space

[0194]

[0195] of dimension

[0196]

[0197] (possibly infinite), and the quantum state is a Hermitian positive semi-definite operator with unit trace:

[0198]

[0199] p e = {p e, p > 0, Tr[p] = 1}. Different laws of evolution describe the trajectory of the quantum state {j°}te[o, T] 'nstate space

[0200]

[0201] 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).

[0202] The ME describes the evolution of the average state ptpt 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 linear ordinary differential equation (ODE): dpt / dt = £t(

[0203]

[0204] pt) = -i[Ht,pt] + Xk=i®Lkt(Pt), where £tis the Liouvillian superoperator, Htis the Hamiltonian of the system, 2)L(JO) =

[0205]

[0206] LpL1' - (l / 2)LtLp - is the standard dissipator and

[0207]

[0208] ..., LN tis a collection of jump operators.

[0209] For simplicity, we consider in the following a single loss channel with timeindependent loss operator L. 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.

[0210] The SME describes the evolution of the state ptwhen the observer continuously measures the loss channel with a detector. The evolution is non-deterministic, it is described by the non-linear stochastic differential equation (SDE): dpt= £t(pt)dt +

[0211]

[0212] M^tpt.dYf), where is a superoperator which depends on a stochastic process dYtthe 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.

[0213] The detector output is a continuous-time signal defined by the rate of change ofthe stochastic process dYtover time: It= dYt / dt. The path followed by the quantum state over time is entirely determined by this measured signal: each time the observer performs a new experiment (labelled j), they measure a particular realisation of the stochastic signal

[0214]

[0215] , which corresponds to a unique trajectory { ^] in state

[0216]

[0217]

[0218] space. These trajectories are called quantum trajectories, and the state ptis said to be conditioned on the information measured by the observer up to time t.

[0219] 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 dYtinvolved, and the form of the measurement back-action superoperatorL.

[0220] 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 stochastic process is modelled by the point process dYt= dNtwith law: P[cWt= 0] = 1 - W[dNt= 1], and W[dNt= 1] = (e + pTr[Lptl ])dt, where here 9 > 0 is the dark count rate (taking into account false clicks), and 0 < p < 1 is the detector efficiency (taking into account missed clicks). The measured signal It= dNt / dt is the rate of change of the counting process Nt= f*dNt, which counts the number of jumps occurring in the time interval [0,t). The measurement back-action is defined by: JvCLp,dN) = ~

[0221]

[0222] p)(dN - (0 + pTr[LpL' )dt), as described for instance in Rouchon, Pierre. " A tutorial introduction to quantum stochastic master equations based on the qubit / photon system." Annual Reviews in Control 54 (2022): 252-261.

[0223] For the diffusive SME, the detector output is real-valued, and continuous in time. This stochastic process is given by the ltd process dYtdefined by: dYt= / pTr[ L + L^ptdt + dWt, where again 0 < p < 1 is the detector efficiency (taking into account measurement imperfections), and Wtis a Wiener process taking independent Gaussian distributed increment. The measurement back-action is defined by:

[0224]

[0225] JvCL(p, dY) = y / p(Lp + pi1' - Tr[(L + L^)p]p)(dY - yfpTr[(L + Lt)Jo]dt), as described for instance in Jacobs, Kurt, and Daniel A. Steck. " A straightforward introduction to continuous quantum measurement." Contemporary Physics 47.5 (2006): 279-303.

[0226] In the jump SME case, the quantum trajectory is discontinuous: the state evolvescontinuously 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 quantum trajectories or measured signals): pt= E[t]. Here E denotes the statistical average over Ntor Wt. 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).

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

[0228] 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. Importantly, the model of a quantum system 1 is typically known apart from the values of p parameters corresponding to various static or dynamical terms, which can be gathered in a vector Qtruee IF. These parameters may appear in different parts of the model, such as the initial state, Hamiltonian, jump operators, dark counts or efficiencies, physical parameters which parameterizes the control signals from the signal generator(s) 11,13, or even filter functions as described in more detail below. Herein the notation {...} may also be used to denote a set of values or parameters which can be conceptually (that is mathematically or numerically) gathered in a vector.

[0229] The dynamics of quantum system 1 is thus describable via a stochastic master equation having a plurality of physical parameters OtrUe- However, as mentioned in above, the continuous signal Itin 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.

[0230] 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 the continuous-time sharp signal Itinto a discrete-time signal {I0,..., IN}. Each value Ikisdefined by integrating Itagainst the transfer function fkof the acquisition chain for the / c-th time bin: Ik

[0231]

[0232] = f The only quantity available to the experimentalist is this discretised measurement record { / 0,

[0233]

[0234] It is particularly important to take this filtering into account when the digitisation time is not negligible compared to the system timescales, as is often the case for superconducting circuits, for example.

[0235] The digitisation is usually performed by averaging the signal over a duration At (a time bin), which is longer than the bandwidth of the various acquisition chain components. In embodiments, the filter function for the / c-th time bin can then be approximated by a rectangular window fkof duration At.

[0236] In the case of the jump SME, the detector typically gives the number of click events over a given time interval, and the filter function fk=

[0237]

[0238] D[kAt7(k+1)At)is the indicator function defined by Dn(t) = 1 if t e 1 and Dn(t) = 0 otherwise, wherein [...) defines the half-open interval such that x e [0,1) means 0 < x < 1. The resulting signal takes discrete values Ike N: Ik=

[0239]

[0240] Itdt = / V[kAt(k+1)At), where the counting process

[0241]

[0242] =t2 dNt counts the number of jumps occurring in the time interval [t1;t2). In the limit where there is at most one jump per time bin At, the signal is binary Ike {0,1}, it is just a sequence of Os and 1s: 00100010..., and the filtering accounts for the inevitably finite time resolution of the detector.

[0243] In the case of the diffusive SME, the signal is typically averaged over some duration At, the filter function is fk= l[kat,k+1)at) with G the gain of the acquisition

[0244]

[0245] chain. The resulting signal is continuous-valued Ike IR: / k= (G / At) £f+1)At / tdt. We note here that there is a choice in the initial enumeration of k. For instance, if the initial time bin is mathematically designated by k = 0, the integral for Ikgoes from / cAt to (k + l)At. However, if the initial time bin is instead designated by k = 1, then integral for Ikgoes from (k - l)At to / cAt, which is entirely equivalent under translation of k. Hereon in, we will typically take the initial time bin as being designated by k = 1, although any alternative designation of the initial time bin can of course be defined under the appropriate translation in k.

[0246] In both cases, we call this digitised signal the binned signal. From a practical point of view, this discrete-time signal Ikis the only quantity available to the operator of quantum system 1.

[0247] Thus, in a first step 301, the detector 8 coupled to the at least one qubit-hosting structure 6 is used to physically measure the continuous signal from quantum system 1to provide a plurality of digitized discrete-time binned signals

[0248]

[0249] { / 1;up to time t = k t, wherein k is an integer denoting the Zc-th time bin of the detector.

[0250] This digitisation process is illustrated in Figure 4 for the diffusive SME, wherein the filtering and digitisation of the sharp diffusive signal Itagainst the rectangular window filter functions fkin Figure 4 results in a discrete-time binned signal

[0251]

[0252] { / 1;

[0253] Returning to Figure 3, in a second step 303, a new object describing the quantum system is computed, which is referred to as the “binned state”, using the discrete-time binned signal {I1,...,Ik}}.

[0254] Specifically, the binned state pkis computed at time t = k t, defined by:

[0255] Pk=IE [PfcAt I ’ ■■■> c]'

[0256] As will be appreciated, this formulation is similar to the form of the ME recovered by averaging over all possible realisations of the stochastic process driving the SME described above: pt= E[t]. Indeed, it is the exact same average on the stochastic process (quantum trajectories) pt, except that this averaging is conditioned on the integral of the modelled continuous signal being what was observed. Phrased differently, the averaging is performed only over the states p that are consistent with the observations (i.e. the measurement record of the discrete-time binned signal {1,..., / k}), not on all possible values of p. This estimate is the best Bayesian estimate of the state from the observer’s point of view, it is the quantum state as usually defined by textbooks, but crucially taking into account of only the information actually acquired by the observer.

[0257] As an example, lets consider the binned state p±and corresponding discrete-time binned signal I±for the first time bin k = 1. All that is known is the value of Z1;such that the integral of possible continuous signals (i.e. trajectories from the SME) must equal this value,

[0258]

[0259] ltdt, where the filter function is the indicator function as defined above and all parameters of the acquisition chain have been set to unity. The binned state is thus p1 =E{ / t}o£t£At[pt({zt}o<t<4t)| = / i] = Jim s / =1? where the only p

[0260]

[0261] kept in the sum are those having a signal which satisfies the condition f^Itdt = / i.

[0262] Generally then, the binned state pkis the average of a plurality of solutions pk&tof the stochastic master equation at time t = k t each for different possible continuous signals wherein, for each of the plurality of solutions ρkΔt, the value of the possible continuous signal integrated over a given time bin is substantially equal to the discretetime binned signal for said given time bin.

[0263] The present inventors have further discovered numerically efficient schemes for computing the binned state pk.In particular, the binned state pkat step k knowing { / 1;..., Ik] is determined by an iterative scheme, starting from the initial state p0= p0(which is known, either from a theoretical model or from an initial experimental characterization, such as an initial state tomography), which takes the following iterative form:

[0264] _ _ ^k(Pfc-i)

[0265]

[0266] Pk~ Tr[Xk(pk_i)I

[0267] Here, %kis a completely positive update map, a type of object which is well known in the art. In this case, %konly depends on the measurement result Ikat step k:

[0268] = dpe^exp^1)At

[0269]

[0270] with £~ipikaso-called generating Liouvillian £jdefined for a function j = -ip^keither by: (i) £Jt= £t+ (eJt- 1)CLfor the jump SME, or (ii) £Jt= £t+ jtCL+;2 / 2 J for the diffuse SME with J the identity superoperator, and where £tis the standard Liouvillian entering the SME describing the dynamics of the quantum system. Here, CL( ) is the correlation superoperator for an operator 0, and is defined either by: (i) CL(O) = 90 + rjLOL^ for the jump SME, or (ii) CL( ) = ^(LO + OL+) for the diffuse SME, wherein the notations take the same forms as defined above. As will be appreciated, we can have many jump operators, but may have only one that we measure, which is the measured jump operator L entering the definitions of the correlation superoperator. Of course, there may also be many measured jump operators (e.g. specially in the case of many detectors) - the “measured jump operator” should be interpreted as a single measured jump operator, or a plurality of measured jump operators, depending on the situation.

[0271] Here,

[0272]

[0273] expkΔt(k-1)Δt(O) denotes the time-ordered exponential of the superoperator 0 (i.e. the generating Liouvillian £~iplk, which is a superoperator, to be applied on the operator pk-to evaluatekpk-y).

[0274] The update map %khere is itself not trace-preserving, however dividing by its trace Tr^C fc-i)] ensures that the overall state update pk= is a completely

[0275]

[0276] positive trace-preserving map.

[0277] Namely, operator p0is the density matrix describing the initial state of quantum system 1 (and operators p, p2,..., pketc. are the subsequently iteratively determined binned states at those time bins), 9 is the dark count rate, p is the detector efficiency, and L is loss or jump operator (which can be time-independent or take different forms at different times). The generating Liouviliian, correlation superoperators, and generating density matrix (introduced below), are rigorously described at least in: (i) P. Guilmin P.Rouchon, and A. Tilloy, “Correlation functions for realistic continuous quantum measurement.” IFAC-PapersOnLine 56.2 (2023): 5164-5170; and (ii) P. Guilmin, P. Rouchon, and A. Tilloy, “Parameters estimation by fitting correlation functions of continuous quantum measurement,” arXiv preprint, no. 2410.11955, 2024, doi: 10.48550 / ARXIV.2410.11955.

[0278] Thus, the binned state pkis calculated iteratively by first setting the initial state as p0= p0at time t = 0. The map

[0279]

[0280] is then applied to p0, such that the next state at the k = 1 time bin is thus calculated via p±= The map %2is then be applied to

[0281]

[0282] i such that the next state p2is calculated via p2=

[0283]

[0284] and so on iteratively up to Tr|JC2( i)J

[0285] the / c-th time bin to provide pk.

[0286] Here, we provide a short derivation of this iterative form for determining the binned state, wherein for simplicity we consider only a single loss channel with a timeindependent jump operator L, a single signal from the detector 8, and simple signal filter function as described in examples above fk= Dk.

[0287] The main ingredient is the aforementioned is the so-called generating density matrix at time t is given by ptJ= E [exp and dYt> = It>dt' (obtainable as

[0288]

[0289] shown in appendix S1 of the P. Guilmin et. al. reference (ii) cited above).

[0290] The generating density matrix at time t, ptJ, obeys the so-called generalised quantum master equation dptJ / dt = £tJ(

[0291]

[0292] ptJ), where £Jtis the aforementioned generating Liouvillian, which is a superoperator.

[0293] Formally, the solution of dptJ / dt is given by ptJ= expo(£ )(o), wherein the notation expo(f) denotes the time-ordered exponential of function f between initial time 0 and time t. In other works, the time-ordered exponential is sometimes written using the time-ordering symbol X: exp^ (f) = Texp

[0294]

[0295] fsds^.

[0296] It will then be appreciated that for the constant function j = alkwhere a e C, the generating density matrix is given by p“lfc= E[eaIkpt]. Accordingly, by “massaging” conditional probabilities, we have:

[0297] Pk=EtPfcat I

[0298] Pk = E[Pfcat I Ik> Pk-l\>

[0299] - _ iE[<y(zfc— ik)pkAt i pk-i]

[0300] p

[0301]

[0302] k~ Em -fo i / vd ’

[0303] where the conditional probability formula has been used P[A I B, C] = P[4 n B I C] / P[B I C], We thus condition over a signal that is either discrete-valued for the jump SME: IkeN, such that 8 denotes the Kronecker delta function, or continuous-valued for the diffusive SME: Ike HR, such that 8 denotes the Dirac delta function.

[0304] The conditioning on the signal value Ikcan be connected with the aforementioned generating density matrix using the inverse Fourier transform of the 8 function: 5(x) = f dpeipx, where the bounds of the integral are -n to n for the jump SME (and 8 is the Kronecker delta function) and -co to oo for the diffusive SME (and 8 is the Dirac delta function). In the below, the bounds of the integral are herein temporarily dropped for notational simplicity.

[0305] The numerator of pk=E[^(rJ„fcrffc),p^tlPfc71]> denoted pk, can be evaluated as follows:

[0306] Pk = IE [8(Ik— IklPk&t I Pk-1\

[0307] = ^r]E^ dpe^-^p^t | pk-i]

[0308] = dpeipIkE[e~ipIkpkAtI pk-1]

[0309] = dpeip,kp^t\

[0310]

[0311] where pk^l kis the generating density matrix defined above for j = -ip^kand initial value jOfc-i at time t = (k - l)At. Thus, by using the above solution of dptJ / dt, namely ptJ= expo( f)(jo0), we retrieve the above formulation of pk= J dpe^exp^

[0312]

[0313] 1)At

[0314] The integral to evaluate the application of %kto pk-can be computed accurately using an efficient numerical quadrature, for example with Gaussian quadratures. Computingk(pk-i) exactly (to machine precision) typically requires evaluating the integrand for a few tens of values of p. For each value of p, the numerical cost is the same as solving the ME (Lindblad master equation). Using vectorized-array computations, we can compute multiple values of p simultaneously, which reduces the numerical cost of the iterative scheme %kto that of solving only a single ME for the same quantum system 1.

[0315] As will be appreciated, quadratures are an exemplary method to approximate integrals by summing a few values of the integrand at specific points: f f(x)dx ~wif(.xd where d is the number of quadrature points,

[0316]

[0317] are the quadrature weights, and xt are the quadrature points.

[0318] Explicitly, for the jump SME, the bounds of the integral are -n to n, such that we have: %k=

[0319]

[0320] dpe‘p / fcexp(kX)At(£ + (e~ip- 1)(?L), with like terms taking the samedefinitions as provided above. For example, one can use the so-called Gauss-Legendre quadrature to compute this integral numerically exactly.

[0321] Explicitly, for the diffusive SME, the bounds of the integral are -co to co, such that we have %k= ^

[0322]

[0323] CoadPeipIk~Tp2expk^_1)&t(JL- ipC^, with like terms taking the same definitions as provided above. For example, one can use the so-called Gauss- 2 Hermite quadrature to compute this integral numerically exactly. Note that the e~vfactor of the quadrature is natively present in this map.

[0324] Moreover, in practice, the Liouvillian may be vectorized to compute these time-ordered exponentials. As this may be undesirably costly for large Hilbert space, so in practice one may only compute the actions of superoperators on states, without every materializing the superoperators. Said otherwise, only %k(pk-i) is computed, and never

[0325] The present Inventors have also discovered that the above expressions for Xfc(jOfc-i) may be computed approximately in a perturbative manner.

[0326] For instance, taking the example of the diffuse SME, wherein ^k(Pk-l) — fZo dp6ipIk~~p2exp(fc1)4t(£ - ipCLPk-i), when At is smaller than the other time scales involved in the dynamics (but not infinitely smaller), one can expand the second (time-ordered) exponential and retain only terms up to a certain order in At. Crucially, the binned signal Ikshould be typically of order VAt at least for sufficiently small windows (after which the growth is linear). Moreover, it is noted that integrals of the form Jk= — J dp pkeip k~~pare of order l / (At)k / 2. This means that the term of order pAt in the time-ordered exponential behaves as a term of order VAt and thus should be expanded to higher order.

[0327] Accordingly, as an example, Taylor expanding the appropriate exponential and ordering contributions with their power of At up to second order in At (and identifying and using the relevant Gaussian integral identities), we have:

[0328] Q~ i-k!^ _ -£ %k= [<7 + VAt ikCLyj2nAt, — 2 ( 1 — ik_ +VAt l£-—^c2

[0329] - 3), I +VAFk-J-CL3+^(CL£ + £CL)

[0330] , — 4 (ik- 6ik+ 3, ik- 1 _ 1 +VAt ( - el + (e2£ + 2CL£CL+ £CL2) + -£2

[0331]

[0332] where the notation ik= Ik / f t has been used for simplicity, and the like-terms take their respective forms defined above for the diffuse SME. As will be appreciated, to compute this numerically, there is no need to compute any more expensive integrals, rather each term of the sum may be simply computed. Thus, the binned state can be computed perturbatively, as shown in the above non-limiting example.

[0333] As will be appreciated, the above formulations for %kmay be generalized to condition on multiple signals, or multiple values of the same signal. That is, the above formulae may be generalised to compute the binned state from signals from different detectors (or multiple signals from the same detector, for instance if continuously measuring the quantum system 1 in a steady state, the same signal can be segmented into a plurality of sub-signals), that can be of jump and / or diffusive type.

[0334] Thus, to condition on N signals in the above coming from detectors (dn...,dN'), or indeed a single detector with the signal partitioned into N sub-signals, j is now a set of test functions, one for each signal j = { / wLzeiii ip with N the number of measured signals. Splitting detectors between those of jump-type n e Sjumpand those of diffusive-type n e Sdiffusive, the new generator of the ODE dptJ / dt = LtJ(ptJ) described above is:

[0335] • 2 p = tt+ £ 0'“ -!) + 2, +

[0336]

[0337] n^jump nES diffusive where like notations take the same definitions as before. Thus, the generating density matrix for multiple functions j = {j\, is:

[0338]

[0339] The above form for the binned state is then straightforward to extend to conditioning on multiple signals. For example, for detector d with binned signal ( / 1;..., Ik) and detector d2with binned signal ( / i< - k)> where the change of notation to J here is merely to distinguish the second detector signals, then:

[0340] Pk=IE[PfcAt I (71, ■■■ < / / <)< (Jl> ■■■ Jk)]

[0341] Pk=EtPfcat I Ik’Jk’Pk-1]

[0342] _ _ IE [<5( / fc ~ Ik)8(Jk ~Jk)Pk&t I Pk-11

[0343] P

[0344]

[0345] k~ Wlk- Ik)(Jk -Jk) \ Pk-l] ■

[0346] Here for simplicity, the time bins At and initial times are assumed the same for both detectors. Now using the inverse Fourier transform for each 8 in the numerator of Pk, denoted pkPk ~ IE [5(4 — Ik)8(Jk ~Jk)Pkbt I Pk-li

[0347] / 1 \2

[0348] =(3-) dpdqeipOk lk)eiq<Jk Jk)pkM I Pk-i] \Z7T /

[0349] / 1 \2

[0350] = JJ dpdqeip / fcel( / 7fc]E[eipIkeiqJlipkMI Pk-i]

[0351] / 1 \2

[0352] = (^) ff dpdqeip‘keiqJkpkJ&t.

[0353]

[0354] where pkJ&t=p^plk~iqlkjsthe aforementioned generating density matrix for multiple functions j = -ip^k,-iq^k') and initial value pk-±at the time t = (k - l)At. Finally, the new completely positive map %konly depends on the measurement results Ikand Jkat step k:

[0355] 1 1 \2

[0356] %k =dpdqeipIkeipJkexpk^tff (£q),

[0357]

[0358] with £Jthe above generating Liouvillian for multiple signals for j = (-ip^.-iq^).

[0359] Accordingly, for N detectors (or same detector at different time durations) (d1;...,dN), each with discrete-time binned signals (ll k,

[0360]

[0361] at step k (where the first index numerates the detectors or signals), we have:

[0362] N \

[0363] / 1 \wr r

[0364] Kk= 7- ■■■ dpi -.dpjv

[0365] \2n / J J neip / nfc)exP(k-i)At(^7)<

[0366]

[0367] with £Jthe above generating Liouvillian for multiple signals for j = f-ip^,...,-ipN^kFor example, for detectors d1, d2) and measured signals (Jkk) at step k, for the SME of jump type we have:

[0368] (1 \2r77r71— dp dqeip!ke^kexP^1)At(£ + (e-ip- l)CL+ (e~iq- 1)CL\

[0369]

[0370] ^71 / J-nJ-nwhereas for the SME of diffusive type we have:

[0371] / 1 \2r°° r°° At?.. At2 / \%

[0372]

[0373] k =J _ (yr.dPJ — (yr.elqJk~q eXPtfc-l)At ~[PCLdl~t(lCLd2)-

[0374] The generalisation is useful when multiple signals are measured, to perform simulation with adaptive time stepping (where one needs to sample an intermediate signal value with the correct statistics), or for filtering functions that overlap over different time-bin, and introduce some superficial memory effect.

[0375] The present Inventors have further recognised that the map %kmay also be used to give the probability of acquiring a discrete-time binned signal Ikknowing the previous binned state pk-via:

[0376] P[ I Pk-i] = Tr[%k(pk_!)].This result is particularly useful for numerical simulations and to optimize quantum state readout protocols.

[0377] That is, if we want to simulate a trajectory from scratch, one would need to sample Ikwith the correct statistics. This is enabled by using the above calculated probability of acquiring a discrete-time binned signal Ikknowing the previous binned state pk-±, PUfc I Pk-l - As an example, using the above perturbative expansion of [%fe(jOfe)], one finds:

[0378] e-Z^ / (2At)

[0379] Wk I Pk-i] = (1 + AIk+ BIk+ C / k3+ DI ft,

[0380]

[0381] \2nt\t where A, B, C, and D are simply the traces of the terms of the corresponding order in Ikthe above perturbative expression. Thus, one can simply precompute A, B, C, and D and then sample from this probability distribution, yielding a simulated Ik, which can thus subsequently be used to simulate pk(as if we had physically measured Ik). This may provide a better numerical discretisation of the SME than the existing schemes (which typically chooses a numerical discretisation time and evolves the state with a discretized scheme, for example a Euler scheme, which is only valid for negligibly small discretisation times).

[0382] Regarding optimizing quantum state readout protocols, as an example, conventional readout protocols are optimized based on the assumption that the probability density of the time- integrated signal (acquired via observing the state in the superconducting readout resonator) is explicitly Gaussian in form. This yields a simple formula for the readout fidelity (the so-called signal-to-noise ratio, SNR). However, it has been recognized that this assumption is not always valid, for instance for fast readout protocols, where the readout resonator state is driven strongly.

[0383] Using the binned state, and also the probability of discrete time signal Ikbased on the preceding binned state P[ / kIk_i], the exact probability density of the time-integrated signal (i.e. the discrete-time binned signal) for a given initial state. The present Inventors have discovered that this allows for refining the SNR (which is just an approximation of the fidelity) and thereby finding the optimal physical parameters (i.e. which parametrizes the input control signal(s) from the signal generator(s) 11,13) for a readout protocol. Crucially, this optimization of a readout protocol is valid even when the final signal distribution is not Gaussian, which is particularly relevant as decay errors occurring during the typically long readout durations contributes to the distortion of the distribution to non-Gaussian forms. For instance, so-called T1 decay errors during the readout result in the excited state being mis-identified as the ground state (e.g. for typicaldiscrete-level qubits such as a transmon wherein the excited state is one computational state and the ground state is the other computational state). The time-averaged signal distribution for the excited state is not gaussian, but bimodal with a small part centered within the gaussian which corresponds to the time-averaged ground state signal.

[0384] Returning again to Figure 3, in a third step 305, the binned state pkis used for a subsequent purpose, either to numerically reconstruct the state, to extract values of physical parameters of the quantum system 1, or to determine the physical parameters for subsequent physical operations.

[0385] For instance, the binned state pkmay be used to either: (a) characterize the quantum system 1; and / or (b) determine physical parameter values for input control signal(s) from the signal generator(s) 11,13 to physically control the dynamics of the quantum system 1 (i.e. in subsequent time); and / or (c) optimize a readout discriminator of a quantum state readout protocol of the quantum state 1.

[0386] For instance, generally, it is expected that each physical parameter of the quantum system 1 has a unique impact on the evolution of the state, and therefore on the measured signals. To estimate (possible simultaneously) the p parameters {0}, first the plurality of digitized discrete-time binned signals [I1,...,lkis measured. As the physical parameters enter the Liouvillian £ 6) and / or the measurement back-actionL(0) of the SME, the binned state is inherently also a function of the physical parameters pk(9. Thus, one may numerically fit the binned state pk(0) to the plurality of digitized discrete-time binned signals {I1,...,Ik}} to assign a value to the p parameters {0}, for example with a least-squares method. The result is a vector of parameters 9fitthat minimises the difference between the experimental results and the binned states.

[0387] Figure 5 illustrates this method for extracting physical parameters. Approximative or exact formulae may be used to compute the binned states, as discussed in above.

[0388] Figure 6 shows a method for determining values for physical parameters which parametrizes the input control signal(s) of signal generator(s) 11,13 to physically control the dynamics of the quantum system 1.

[0389] For instance, at step 600, a desired final (target) state is provided. At step 601, the binned state is determined from the measurement record of discrete-time binned signals as described herein (i.e. for an initial operation for optimization to be performed on quantum state 1 having at least approximately known physical parameters or for a known operation to be performed on quantum state 1 having initial physical parameters for optimization).At step 602, a fidelity between the target state and the binned state is calculated, which may generally be used to define a loss function fLoss. The loss function fioss typically depends on one or more target operating attributes Xtargetto be reached by the quantum system 1, but can also depend on some intermediate stage of the dynamics. Typically, the target operating attribute Xtargetis a target quantum state Ptargetora target quantum unitary Utargetto be reached by the quantum system.

[0390] For instance, a possible objective is to prepare a target quantum p^t3et- for instance a Fock state with five photons p^t3et= |5><5| reached by time t = kM. In such a case, the following loss function fLosscan be built: fLoss(

[0391]

[0392] p^t3et<0) =1“ where: F is the fidelity between two quantum states.

[0393] The fidelity is a known measure of mathematical distance between two quantum states, and more particularly between two density operators. The fidelity between two density operators p and a is commonly defined as follows: F(p, o) = (tr J / ptj / p)2. Consequently, the previous loss function fLosscan be seen as a measure how different the binned quantum state pk(9 defined by the p-dimensional parameter vector 9 and the target state

[0394]

[0395] are from each other. By minimizing fLoss, we bring the experimental state closer to the target state.

[0396] As mentioned above, a different objective can also be to prepare a target quantum unitary Utarget. Such a target quantum unitary Utargetcan be defined as a mapping from a set of initial quantum states {p,..., Pm} to a set of respective target quantum states {pf..., pf} (i.e. a quantum gate). This mapping can be written as follows for every /

[0397]

[0398] e [1, m]|: p* = Utargetp\U^arget. In such a case, the following loss function fLosscan be built: f

[0399]

[0400] L(Utarget,9) = m-^=iF(UtargetpiU^arget, U Preparing a target quantum state or a target quantum unitary are common objectives. However, it is known to the skilled person that other objectives can be achieved and that other target operating attributes can be defined.

[0401] Thus, with the loss function defined, a local minimum {9min} of the loss function within the space of physical control parameters of the input control signals may be determined.

[0402] One way may be, for instance, to compare the calculated fidelity with a threshold value, as shown at step 603. If the fidelity is below the threshold value, then the physically controllable parameters of operation (namely the set of physical parameters which parametrizes the input control signal(s)) are modified in SME in step 604. The binnedstate is computed once again 601, and the subsequent fidelity is calculated once again 602, and compared with the threshold value in step 603. The steps 601 -604 are repeated until a calculated fidelity is above the threshold value, in which case in step 605 the parameters of the final iteration, corresponding to physically controllable parameters of the operation, are extracted and used as the set of physical parameters of the input control signal(s) in a subsequent control operation of the quantum system.

[0403] The fidelity may be given as the overlap between the approximated final state and desired final state. The threshold may be between 0.7 (70%) and 1 (100%). For instance, the threshold may be chosen such that the fidelity between the approximated final state and the desired final state is above 0.8, preferably above 0.9, and most preferably above 0.99. Moreover, as will be appreciated, modifying parameters between successive iterations of choices of SME may use any known optimization algorithm, such as gradient descent.

[0404] Finally, the binned state may also be used to optimize the discriminator in readout protocols of the quantum system 1 (e.g. qubit readout protocols when the quantum-state-hosting structure 6 is for hosting qubits).

[0405] For instance, in the example of the above cat qubits, when preparing the cat state in the parity state s = It?1) and performing a measurement (i.e. state readout) operation, the acquired plurality of digitized discrete-time binned signals SL=

[0406]

[0407] { / , during the readout operation will follow a distribution

[0408]

[0409] in the space of signal traces (where i denotes the iteration of the measurement). Thus, we can define a signal trace (the discrete-time binned signals) by the set of k complex values it took at the corresponding k time steps. So the signal trace space is U = C".

[0410] A single-shot signal trace discriminator is a function f of the signal trace space to a binary output set {'+7- '} which takes a single-shot signal trace as an input, and returns a measurement label, in our case ‘+’ or as an output. The discriminator can be seen as a binary partition of the signal trace space.

[0411] For each prepared state (in our case s = It?1)) the associated discriminator error£s ’■= P(f(. S') =7 s'|s) is defined as the probability that the signal trace S, which follows the random distribution

[0412]

[0413] in the space of signal traces, gets discriminated to an output label f (S) which does not correspond to the prepared state s. The mean error of the discriminator / is defined as (e++ e_) / 2.

[0414] For a given measurement process, there is a single discriminator functionoptthat minimizes the mean discriminator error. This optimal discriminator is defined asf r<?>1 =[ / + / ifs'+(s)>s'-(s)

[0415]

[0416] optI ' — ' otherwise

[0417] where gsis the probability density function of a random variable following the probability distribution 5S.

[0418] The binned state can be used to find the optimal discriminator by computing, for a specific choice of system and control parameters, the exact value of the probability density function g+and g_ for the digitized discrete-time signal measured experimentally. The system and control parameters can then be optimized numerically to find the discriminator with the optimal fidelity fopt(S).

[0419] Numerical examples

[0420] As discussed above, by preparing and measuring the quantum system 1, certain parameters can be more precisely and / or accurately fit as opposed to conventional methods relying entirely on the SME as the model with which to compare experimental data.

[0421] For example, let’s consider a qubit with Hamiltonian H = | ox+ 1 oyand a single

[0422]

[0423] jump operator L = o_. We consider that the jump operator is measured with perfect diffusive continuous measurement of L (with g = 1.0), and we simulate the system evolution for t e [0, T] with T = 2.0 and starting from initial state p0= |e><e|.

[0424] We can plot the average Lindblad evolution by showing the expectation value of all Pauli matrices, as shown in Figure 7.

[0425] If we are interested in the probability distribution of the signal averaged over a long time I =

[0426]

[0427] dYt, one option is to simulate many SME trajectories with a finely grained time-step, then averaging the simulated signal between [0, T] for each trajectory, and group the results by their value to estimate the probability distribution with an imperfect histogram. But using the binned state, we can directly compute the exact probability of the signal without simulating any trajectory. In Figure 8, we compare for this qubit system the estimated density probability (for 1000, 10,000 and 100,000 simulated trajectories) with the exact probability computed with the binned state.

[0428] Here, the “theory” curve in Figure 8 has been computed in double precision (float64), i.e. 101 points (x axis) where for each point we compute a degree 31 quadrature takes 0.5ms on a CPU (using vectorized computation).

[0429] Now, we can also plot the theoretical probability distribution for various values of g, as shown in the upper panel of Figure 9. This can be used to numerically fit the valueof the efficiency, by comparing the measured experimental distribution with the theory prediction, e.g. using a simple least-squares method. More generally, any parameter that has an impact on the time-averaged signal probability distribution can be fitted using this method.

[0430] We can also study the probability distribution for various values of the integration time T, shown in the lower panel of Figure 9.

[0431] As another numerical example, let’s consider the reduced system for the dissipative cat qubit (e.g. a system as described above in relation to Figure 2), after adiabatic elimination of the buffer. We consider the so-called “deflate” protocol, where the buffer drive is turned off, and the memory state empties through the buffer. The Hamiltonian is H = 0 and there is a single jump operator L = a2. \Ne consider that the jump operator is measured with imperfect diffusive continuous measurement of L (with p = 0.5), which corresponds to homodyne measurement of the buffer along the X quadrature, and we simulate the system evolution for t e [0, T] with T = 2.0 and starting from initial state a coherent state with a = 3.

[0432] We can plot the average Lindblad evolution by showing the state Wigner at various times, in the upper panel of Figure 10.

[0433] We then propose to compare the existing schemes to reconstruct the system state from a measurement record. To do so, we simulate the SME evolution for a state ptwith a very fine-grained 8t = 10-3, we then integrate the corresponding signal over time bins of duration At = 0.1, and try to reconstruct the system’s state pkfrom this integrated signal, with a Rouchon unnormalized scheme (CP to machine precision but TP to order one, just renormalised with the trace which is the wrong projection, comparable to the Euler scheme), the Rouchon normalized scheme (CPTP up to machine precision), or with the proposed binned state. We repeat this process for 100 trajectories, and compare (on average over these 100 trajectories) the fidelity of the state reconstructed from the integrated signal pk(for each numerical scheme) with the true state pt. There are 20 time bins of duration At e [0, T], so we have 20 states to compare. References for Rouchon schemes can be found in “A tutorial introduction to quantum stochastic master equations based on the qubit / photon system” cited above. The results are shown in Figure 11.

[0434] As can be seen, the fidelity of reconstruction is much higher for the present Invention dubbed “Robinet” (> 0.9) than the two other schemes, especially at short time where the deflate dynamics is fast compared to the digitization timescale.

[0435] We can also compare the average state for each reconstruction. The Rouchon(unnormalized) state averaged over the 100 trajectories is shown in the second panel of Figure 10. As can be seen, the unnormalized Rouchon scheme the states makes no sense, it’s completely false.

[0436] The Rouchon (normalized) state averaged over the 100 trajectories is shown in the third panel of Figure 10. The normalized Rouchon scheme gives a physically meaningful state, but the evolution is quite far from the correct evolution, and reaches a fidelity as low as 0.25 for small times.

[0437] The Robinet state averaged over the 100 trajectories is shown in the lower-most panel (fourth) of Figure 10. The Robinet reconstructed state is relatively close to the true state. It has a lower average purity, but this is expected as we lost some information compared to the true state trajectory.

[0438] As one can see, all the previous existing numerical methods used to reconstruct the system’s state from the measurement record fail when the digitization time is not negligible. The fidelity between the binned state and the “true” system state (for a very finely grained signal) is not one, but this is expected as we lost some information when averaging the signal. The binned state is the best prediction one can make on the system’s state given the digitized measurement signal.

Claims

Claims1. A method of characterizing or controlling a quantum system, wherein the dynamics of the quantum system is describable via a stochastic master equation having a plurality of physical parameters, and wherein the quantum system comprises: (i) one or more quantum-state-hosting structure(s); (ii) a detector coupled to at least one of the quantum-state-hosting structure(s) and configured to measure a continuous signal output from the quantum-state-hosting structure(s) to provide a measurement signal from the quantum system to a command circuit, wherein the detector is configured to digitize the continuous signal over time bins of duration At; and (iii) one or more signal generator(s) coupled to the at least one of the quantum-state-hosting structure(s) and configured to input control signal(s) to physically control the dynamics of the quantum system; the method comprising:(I) physically measuring, with the detector, a continuous signal from the quantum system to provide a plurality of digitized discrete-time binned signals { / 1;..., Ik} up to time t = / cAt, wherein k is an integer denoting the Zc-th time bin of the detector;(II) computing, by the command circuit, a binned state of the quantum system at time t = k t, pk= E[feat| / 1;..., Ik], wherein the binned state pkis the average of a plurality of solutions ρkΔtof the stochastic master equation at time t = k t each for different possible continuous signals wherein, for each of the plurality of solutions ρkΔt, the value of the possible continuous signal integrated over a given time bin is substantially equal to the discrete-time binned signal for said given time bin; and (III) using, by the command circuit, the binned state pkto:(a) assign a value to at least one of the plurality of physical parameters which characterizes the quantum system; and / or(b) assign a value to at least one of the plurality of physical parameters which parametrizes the input control signal(s) to physically control the dynamics of the quantum system by the one or more signal generator(s); and / or(c) optimize a readout discriminator of a quantum state readout protocol.

2. The method of claim 1, wherein the detector has binary output with dark count rate 9 > 0 and detector efficiency and 0 < p < 1, such that the stochastic master equation describing the dynamics of the quantum system is of the jump type, and wherein operation (II) comprises computing the binned state pkby:(i) assigning an initial binned state p0= p0at time t = 0 describing the initial state of the quantum system; and(ii) calculating the binned state of the quantum system at time t = kAt via the iterative formula pk=wherein %kis a completely positive map and is givenby dpe^exp^1)4t(L + {e~^ - l)eL), wherein exp^^CO) denotes the time-ordered exponential of the superoperator 0, L is the Liouvillian of the stochastic master equation describing the quantum system, and CLis the correlation superoperator defined for an operator 0by eL(O) = 90 + pLOL]', with L being the measured jump operator of the stochastic master equation describing the dynamics of the quantum system.

3. The method of claim 2, wherein the action of the completely positive map %kon the binned state pk-1at the / c-th time bin, %k(pk-i), is computed using a numerical quadrature, and preferably using the Gauss-Legendre quadrature.

4. The method of claim 1, wherein the detector has an output taking a continuous range of values with detector efficiency and 0 < p < 1, such that the stochastic master equation describing the dynamics of the quantum system is of the diffusive type, and wherein operation (II) comprises computing the binned state pkby:(i) assigning an initial binned state p0= p0at time t = 0 describing the initial state of the quantum system; and(ii) calculating the binned state of the quantum system at time t = kAt via the iterative formula pk= wherein %kis a completely positive map and is givenby %k= ^J”oodpeip / fc" TP2exp^_t1)At(2;- jpCL), wherein exp^^C0) denotes the time-ordered exponential of the superoperator 0, L is the Liouvillian of the stochastic master equation describing the quantum system, and CLis the correlation superoperator defined for an operator 0 by CL(O) = (LO + O, with L being the measured jump operator of the stochastic master equation describing the dynamics of the quantum system.

5. The method of claim 4, wherein the action of the completely positive map %kon the binned state pk-at the / c-th time bin, Ckpk-1), is computed using a numerical quadrature, and preferably using the Gauss-Hermite quadrature.

6. The method of claim 1, wherein operation (II) comprises computing the binned state of the quantum system at time t = kAt, pk, perturbatively.

7. The method of claim 6, wherein the detector has an output taking a continuous range of values with detector efficiency and 0 < p < 1, such that the stochastic master equation describing the dynamics of the quantum system is of the diffusive type, and wherein computing the binned state of the quantum system perturbatively comprises:(i) assigning an initial binned state p0= p0at time t = 0 describing the initial state of the quantum system;(ii) calculating the binned state of the quantum system at time t = kAt via the iterative formula pk= wherein %kis a completely positive map and is givenby using:— ifc / 2 _ i (a)Kk =£££M [V + VAt or (b) ^ / < = 7^ [J + VAt i.^ + VAt (£ _ ^eL2)];Or& = St W++ ^t2(£ - ^CL2) + yfti3(^21 d+A ' (CLL +reL))]; or(d) = ^2= [J + VAt\kCL+ VAt2(£ + VAt3(l^lcl + (CL£ +ACj ] + VAt4(CL2£ + 2CL£CL+ £CL2) +wherein is the Liouvillian of the stochastic master equation describing the quantum system, and CLis the correlation superoperator defined for an operator 0 byCL O') = LO + O, with L being the measured jump operator of the stochastic master equation describing the dynamics of the quantum system, and ik= Ik / t.

8. The method of claim 1, the quantum system comprises a plurality of detectors each coupled to at least one of the quantum-state-hosting structure(s) and configured to measure continuous signals output from the quantum-state-hosting structure(s) to provide a plurality of measurement signal from the quantum system, wherein the detector is configured to digitize the continuous signal over time bins of duration At, the method comprising:in operation (I), physically measuring, with the detectors, a plurality of continuous signals from the quantum system to provide a plurality N of sets of digitized discrete-time binned signals [In l, -, In,k} up to time t = / cAt, wherein k is an integer denoting the Zc-th time bin of each detector, and wherein index n e [1; / V] identifies each detector; and in operation (II) computing the binned state pkby:(i) assigning an initial binned state p0= p0at time t = 0; and(ii) calculating the binned state of the quantum system at time t = k t via the iterative formula pk= wherein %kis a completely positive map and is givenby Xfe =(^) / ■■■ / dPi ■■■dPw (nn=!e‘p / n'fc) exp^1)At(£7), wherein J... JdP1...dpwdenotes the multiple integral across variable pn, whereinexpkΔt(k-1)Δt(O) denotes the time-ordered exponential of the superoperator0, is a generating Liouvillian superoperator depending on the binary or continuous outputs of each of the detectors.

9. The method of claim 1, wherein the dynamics of the quantum system is configured to be in a steady state, the method comprising:in operation (I) physically measuring, with the detector, a plurality of sequential continuous signals from the quantum system to provide a plurality N of sets of digitized discrete-time binned signals {In l, -, In,k} up to time t = / cAt, wherein k is an integer denoting the Zc-th time bin of each set, and wherein index n e [1; / V] identifies each set; andin operation (II) computing the binned state pkby:(i) assigning an initial binned state p0= p0at time t = 0; and(ii) calculating the binned state of the quantum system at time t = k t via the iterative formula pk= wherein %kis a completely positive map and is givenby Xfe =(^) J -JdPi -dPN (rin=i efp / "fc)exp^1)at(£-'), wherein J... JdP1...dpwdenotes the multiple integral across variable pn, whereinexpkΔt(k-1)Δt(O) denotes the time-ordered exponential of the superoperator 0, £jis a generating Liouvillian superoperator depending on the binary or continuous output of the detector.

10. The method of any of claims 2-9, wherein the action of the completely positive map %kon the binned state pk-at the / c-th time bin, %k(Pk-i), is computed byevaluating the integrand for a plurality of values of p, wherein at least some of the integrands are evaluated simultaneously using a vectorized-array computation.

11. The method of any preceding claim, wherein:(A) the plurality of physical parameters comprises a first set {^i} of one or more physical parameter(s) in the Liouvillian J3J) and / or the measurement back-actionL(ei) of the stochastic master equation describing the quantum system and / or a second set {e2} of one or more physical parameters of the initial state p0(e2) at time t = 0 describing the initial state of the quantum system such that the binned state pkis a function of the first {^i} and / or second {e2} sets, pk(d 2), and wherein operation (II l)(a) comprises:numerically fitting the binned state pk(d12) to the plurality of digitized discrete-time binned signals {I1,...,Ik}} to assign a value to at least one physical parameter of the first {^i} and / or second {e2} sets; and / or(B) the plurality of physical parameters comprises a third set of one or more physical parameters {03} which parametrizes the input control signal(s) of the one or signal generator(s) such that the binned state pkis a function of the third set {e3},fe(03), and wherein operation (lll)(b) comprises:defining a loss function as a function of the fidelity between a target quantum state p^t3etand the binned state pp̄k(θ3) at time t = kAt,determining a local minimum {6min} of the loss function within the space of physical control parameters of the input control signals; andusing the local minimum {6min} as the set of physical parameters of the input control signal(s) in a subsequent control operation of the quantum system.

12. The method of any of claims 2-11, wherein operation (III) comprises calculating the probability of the discrete-time binned signal for the Zc-th time bin using the binned state pfc-i at the (k - l)-th time bin using the formula P[ / kI pk-] = Trl fcC fc-i)].

13. A quantum characterisation or control unit comprising:a quantum system which comprises:(i) one or more quantum-state-hosting structure(s);(ii) a detector coupled to at least one of the quantum-state-hosting structure(s) and configured to measure continuous signals output from the one or more quantum-state-hosting structure(s) to provide measurement signals from the quantum system, wherein the detector uses time bins of duration At; and(iii) one or more signal generator(s) coupled to one or more quantum-state-hosting structure(s) and configured to input control signals to physically control the dynamics of the quantum system; and a command circuit configured to: (a) control the one or more signal generator(s) so as to physically control the dynamics of the quantum system, and (b) control the detector so as to provide measurement signals from the quantum system;wherein the quantum characterisation or control unit is configured to perform the method according to any preceding claim.

14. A computer program product comprising instructions which, when executed by the quantum characterisation or control unit of claim 13, cause the quantum characterisation or control unit to perform the method of claims 1-12.

15. A computer-readable medium having stored thereon the computer program product of claim 14.

Images