Simultaneous estimation of parameters and the state of an optical parametric oscillator system

In this paper, we consider the filtering problem of an optical parametric oscillator (OPO). The OPO pump power may fluctuate due to environmental disturbances, resulting in uncertainty in the system modeling. Thus, both the state and the unknown parameter may need to be estimated simultaneously. We formulate this problem using a state-space representation of the OPO dynamics. Under the assumption of Gaussianity and proper constraints, the dual Kalman filter method and the joint extended Kalman filter method are employed to simultaneously estimate the system state and the pump power. Numerical examples demonstrate the effectiveness of the proposed algorithms.


I. INTRODUCTION
Quantum estimation is at the heart of many research areas including quantum control, quantum computation and quantum metrology [1]- [3].In quantum state estimation, the aim is to estimate the state of a quantum system given measurement data.Various studies have been presented on both static state estimation and the tracking of a dynamical quantum state [1], [2].For a quantum state estimation problem, we usually assume that the system is well modeled.However, disturbances due to environmental fluctuations or experimental settings may lead to inaccuracies in system modeling [4].Thus, various studies have been done on parameter estimation which aims to estimate unknown parameters in the system modeling from obtained information [5]- [10].Further studies considered that both the state and parameters in a system model should be estimated simultaneously, motivated by either the need to estimate the system state robustly or to estimate the parameters of interest.Simultaneous state-parameter estimation has potential applications in modeling, identification and prediction [11]- [16].
In quantum research, the simultaneous state-parameter estimation problem may also have potential applications in the detection of a classical field by using a quantum sensor and in robust quantum state estimation [17]- [23].A series of works on the state-parameter estimation of a quantum system have already been presented [18]- [23].For example, in [18], [19], the authors modeled the unknown parameter by using a quantum analog system.The works in [21]- [23] employed the quantum-classical Bayesian inference method to solve faulttolerant quantum estimation problems.
In this paper, we employ classical filtering theory for the state-parameter estimation of a quantum system to provide more flexibility on the choice of filtering methods.We consider a nondegenerate optical parametric oscillator (OPO) system which is one of the most interesting and widely used devices in quantum optics [24].The pump power, a key parameter of an OPO system, may fluctuate due to two influences: unwanted disturbances from environmental noise; or by design the pump power is subject to an external signal of interest [4], [17].In classical filtering theory, several algorithms have been proposed for simultaneous state-parameter estimation problems.For example, the dual Kalman filter (dual-KF) method was first proposed in [11] for linear systems and then developed for nonlinear systems in [12], [13].The joint extended Kalman filter (joint-EKF) [14], [15] combines the quantum state and the unknown pump power into a single joint vector.Thus, the combined system becomes nonlinear and the extended Kalman filter (EKF) is used to linearize the system.Based on existing works [14], [15], [25], the joint-EKF is usually expected to be on average more economic than the dual-KF while it may suffer to the divergence problem.Moreover, the joint-EKF is often more sensitive to factors that increase the estimation error due to the approximation procedure imposed by the linearization of the EKF.
Note that classical filtering theory can not be applied to a quantum system directly since the quantum conditional expectation can not be defined properly in a classical probability space due to the uncertainty principle.However, we can find an equivalent classical problem under the assumption that the system is linear Gaussian and with proper constraints on the classical analog system.Thus, both the dual-KF and the joint-EKF can be employed to update estimates of both the state and the unknown parameter.Our main contribution is to formulate the state-space representation for an OPO system, map the quantum filtering problem to its classical analog and demonstrate the efficacy of the two employed methods for the OPO system with a dynamic parameter.Numerical results show that both the dual-KF and the joint-EKF can achieve state-parameter estimation with better performance compared to the case where no filter algorithm is applied to the unknown parameter.
This paper is organized as follows.In Section II, we formulate the system dynamics of an OPO system using state-space representation.In Section III, we first present analysis on a quantum filter and its classical analog.The dynamics of the time-varying pump power are given.Then, the dual-KF method and the joint-EKF method are employed for the simultaneous state-parameter estimation.In Section IV, we investigate the performance of the employed algorithms.Section V concludes the paper.

II. STATE-SPACE REPRESENTATION OF AN OPO SYSTEM WITH HOMODYNE MEASUREMENT
We consider an OPO system consisting of a cavity and a nonlinear medium where the nonlinear effects caused by the medium can be enhanced by the cavity [24].In this section, we describe the dynamics of the OPO system using a state-space representation.
The cavity has two input-output channels, which is a common configuration in experiment [24].The first channel with decay rate γ 1 goes to the measurement.Here, a beamsplitter is added to represent the inevitable measurement loss and noise (see Fig. 1).The second channel with decay rate γ 2 is often used to model the loss inside the cavity.The total decay rate is γ = (γ 1 + γ 2 ).
A nonlinear medium can be characterized by the following Hamiltonian where â is the annihilation operator of the cavity and b is the annihilation operator of the pump beam.The coefficients χ is the second order susceptibilities of the crystal medium [24].Therefore, the dynamics of the nonlinear medium is Fig. 1: Schematic of an OPO system.The cavity contains 4 mirrors and a nonlinear medium (the gray square inside the cavity).The red line is abstract representation of the cavity mode.There are two inputs â1 and â2 and the corresponding outputs âout and âlc .The output âout of the cavity is fed into a beamspiltter which yields two beams âm and âlb .â0 is vacuum noise.
Under proper assumptions, one can replace the operator b in (2) by the c-number such that where β = b and • indicates the quantum expectation [28].
It can be seen that depends on both the nonlinear medium and the pump power [24].Therefore, the Hamiltonian of the OPO system can be given as where ω p is the angular frequency of the pump laser and ω r is the cavity resonance angular frequency [24].We assume that the pump and cavity are tuned so that ω p = 2ω r [26].We move to a rotating frame and the dynamics of â are as follows ( Here, we use the corresponding differential form to represent the noises (e.g., d Â1 = â1 dt).In this paper, we are only interested in the amplitude of the pump power.Thus, can be regarded as a real number.
The system consisting of the OPO and the beamsplitter can be regarded as a 3-input-3-output system with corresponding inputs â0 , â1 , â2 and outputs âm , âlb , âlc .The input-output relation of the system is âout where T ∈ [0, 1] is the transmittance of the beam splitter, which corresponds to the measurement efficiency.Let q and p denote the quadratures Then, we have the following output quadratures where Assume that the following Homodyne measurement is applied to the output âm , where θ m is the phase of the Homodyne measurement.Let x = ( q, p) T denote the system state.The dynamics of x and the measurement ŷ can be described by the following state-space equations, where The quantum covariance of two operator vectors ô1 and ô2 is defined as Denote the correlation matrices as Let V denote the covariance matrix V = Cov( x, x).Then, the Heisenberg uncertainty principle [27] gives Since the commutation relation of quadratures gives [ q, p] = ih, the uncertainty principle (14) reads which can be rewritten as where ] is called the symplectic matrix [2].In our case, the symplectic matrix is which yields

III. FILTER FOR BOTH THE STATE AND THE TIME-VARYING PUMP POWER
In this section, we consider the situation where the pump power is a random process.To achieve the simultaneous state-parameter estimation, we first provide standard quantum filter for the system state.Then, we combine the state filter and the classical parameter filter after an explanation of the existence of a classical analog to our quantum filtering problem.
Given the dynamics of the OPO system in (10), a filtered quantum state conditioned on a measurement record can be obtained by using quantum filtering theory [28].However, the estimation of the classical pump power and the quantum state can not be unified due to the difference between quantum and classical mechanics.Fortunately, there exists a classical analog of the quantum filtering equations given that the quantum system is linear Gaussian [2], [29], [30].We first provide a treatment on the equivalence of a quantum conditioned state and its classical analog.Then, the dual-KF method and joint-EKF method are employed to solve the simultaneous stateparameter filtering problem based on the quantum-classical equivalence.

A. Filter for the state
Given the system Hamiltonian in (4) and the measurement in (9), the conditioned dynamics of the system can be obtained by using the standard quantum filtering theory [28].Meanwhile, a quantum state can also be characterized by a Wigner distribution which is a pseudo-probability distribution in the phase space over a classical configuration corresponding to the quantum configuration [2].The Wigner function appears like a joint classical probability distribution in classical cases.A quantum system can be described as a linear Gaussian system if its Wigner function is Gaussian and dynamics are linear.The Wigner function is positive-definite for a Gaussian quantum state while it can be negative for general quantum states.For a linear quantum system, the Gaussianity can be preserved since future states are linear combinations of the initial Gaussian state.Therefore, the dynamics of the system (10) are completely described by the time evolution of the first and second statistical moments of the quadrature coordinates [2], [30]- [32].In our case, the moments of the unconditioned Gaussian state are Since the state-space equations ( 10) are linear, with the assumption that the initial state is Gaussian, the conditioned state remains a Gaussian state with the following moments where dw = ydt −C x c dt is the innovation.Here, the subscript c means the quantity is conditioned on a measurement record.Note that both (10) and ( 20) are isomorphic to those of a classical linear Gaussian system.Then, we can find the following classical system analog to the system (10) [2] where x = (q, p) T is a vector of classical random variables and y is a classical random variable.v is a classical Wiener process.The coefficient matrices are given in (11).
Here, we present the main restrictions applied to the classical system (21) inherited from the quantum origin.The two restrictions are originated from the unitary evolution of the unconditioned system and the Heisenberg uncertainty principle.For the unconditioned system (21), the unitarity places the following fluctuation-dissipation restriction on the drift and diffusion matrices A and D For the conditioned state, we have the following fluctuationobservation relation which preserves the uncertainty relation ( 16) Since we assume that the inputs are vacuum states, which are Gaussian, the covariance for the vacuum noises is h 2 I where I is the identity.The Heisenberg uncertainty relation for the noise vector v gives [33] where

B. Time varying pump power
In optical cases, the pump power may fluctuate due to environmental perturbations or instrumental settings.Here, we assume that is described by the following stochastic process where µ < 0 is the drift coefficient and g is the diffusion coefficient.v is a classical Wiener process.The constant c characterizes the expected value of the stochastic process and we denote it as the tendency constant.Thus, the corresponding state-space equations are The covariance of two vectors o 1 and o 2 of classical random variables is defined as where E(•) denotes a classical expectation.The correlation matrices are D dt = Cov(gdv , gdv ),

C. The dual-KF method
In this section, we briefly revisit the dual-KF method from the Bayes estimation perspective and provide the continuous dual-KF algorithm.Given ( 21) and ( 27), we aim to obtain estimates of both and the state x.The task can be achieved by using the Maximum a Posteriori (MAP) method which aims to maximize the joint conditional density ρ x, |y [13], [25].The dual-KF method separates the conditional density into two terms ρ x, |y = ρ x| ,y ρ |y .
The estimate c is found by maximizing ρ |y while the estimate x c is found by maximizing ρ x| c ,y .The corresponding cost function for the state x is where ] is the optimal prediction.Since the probability density is Gaussian, minimization of the above cost function regarding x can be reached by using the Kalman-Bucy filter (20).
Similarly, the cost function for the estimation of parameter is where y − = E[y|y t−dt 0 , t−dt 0 ] and − c are the optimal predictions.The cost can be minimized using the following Kalman-Bucy filter: where dη = ydt −Cx c dt is the innovation and is the linearization coefficient which can be calculated using (20) while the conditioned x c is used to approximate x.At every time step, the conditioned state x c is updated using the current estimated parameter c while the estimated c is also updated using the current state x c .Here, we summarize the algorithm as follows.
Algorithm 1: (continuous dual-KF) Initialized with: Update for the state: x .Update for the pump power:

D. The joint-EKF method
The joint estimation aims to maximize the following conditional density To generate the above MAP estimates, we define the following new joint state consisting of both x and The state-space representation of the joint state z is where The noise correlation matrices are (39) The generalized system becomes nonlinear and an EKF be applied to generate an approximation of the MAP estimation.We denote this method as the joint-EKF method and the algorithm is as follows.

IV. NUMERICAL RESULTS
In this section, we first apply the proposed algorithms to a special case where parameters are fixed to values that are close to real experiments.Here, we aim to show the effectiveness of the two methods without comparing their performance.Then we test the performance improvement regarding the measurement efficiency, the tendency constant and diffusion coefficient of the classical signal.

A. Case study
In this section, simulation results of both the dual-KF method and the joint-EKF method are presented.For comparison, we simulate x c using the dual-KF method, the joint-EKF method and the KF method.For the KF method, the Kalman filter is applied to the system state while no filter is applied to and we assume that = c is a constant.x T represents the true evolution of the system state.Here, x T is simulated using (20) given known and a complete measurement record of all three outputs âm , âlb and âlc .See the Appendix for a detailed explanation of the simulation of x T .
To quantify the performance improvement, we define the relative performance improvement (RPI) for as The subscript c indicates either the dual-KF method or the joint-EKF method.The RPI of c evaluates the percentage of the MSE reduced by using the two proposed methods compared with the KF method.We simulate the system for N times and characterize the RPI using its mean value I m c and the standard error of the mean (SEM) S which are defined as follows The relative improvement for the quantum state I x c = (I q c ; I p c ) is defined similarly and the corresponding mean and SEM are calculated in the same way.
The following parameters are selected for the case study.
• θ m = π 12 rad, • γ 1 = 0.95 rad/s, γ 2 = 0.05 rad/s, • c = 0.5, µ = −0.01rad/s, g = 0.028.Note that the escape efficiency γ 1 /γ is based on experiments [34].We use the total decay rate as γ = γ 1 + γ 2 = 1 rad/s for simplicity, which determines the scaling of time evolution.The actual values in the experiments are in order of 2π × 10 × 10 6 rad/s.Fig. 2 presents time evolution of one trial of the simulation.Fig. 2 (a) gives the time evolution of the estimated c .The black line is the true evolution of .The green line is straight which means that is a constant value since there is no filter Fig. 2: One estimation trial of (a) c , (b) q c and (c) p c by using the dual-KF method, the joint-EKF method and the KF method.
applied to for the KF method.The red and blue lines are generated by the dual-KF method and the joint-EKF method, respectively.It can be seen that on average the blue and red lines are closer to the black line compared with the green line, especially when there is a large deviation of the real value of to the tendency constant c.Fig. 2 (b) and (c) give time evolution of q and p of the dual-KF, the joint-EKF and the KF methods, respectively.The red line is generated by using the dual-KF method and the blue line is generated by using the joint-EKF method.The KF method (the green line) does not take any information from measurement to update c .The RPIs for the dual-KF method are 73.4% ( c ), 79.9% (q c ) and 75.1% (p c ) while the RPIs of the joint-EKF method are 48.6% ( c ), 71.8% (q c ) and 69.5% (p c ).It is clear that both the proposed two methods can follow the true evolution of the system state x better than the KF method.dual-KF joint-EKF I m c (%) 48.6 38.5 I m q c (%) 51.1 42.9I m p c (%) 39.5 34.9 TABLE I: Mean RPI of the dual-KF method and the joint-EKF method.The corresponding SEM is at most 0.92%.
While the performance improvement of the proposed methods may not be convincing with only one trial result, the superiority can be confirmed by Table I which gives the mean RPI of the dual-KF and joint-EKF methods calculated over N = 1 × 10 3 trials.The result shows that both the dual-KF method and the joint-EKF method have a performance improvement of over 34.9% while the dual-KF method provides higher RPIs on average.We also consider different initial states for both the vacuum state and the coherent state including (c) Fig. 3: RPI of (a) c , (b) q c and (c) p c with the measurement efficiency T increasing from 0 to 1.
x 0 = (0 0), x 0 = (1.40), x 0 = (0 1.4) and x 0 = (1.41.4).The RPIs are consistent for the four different initial states.Thus, we conclude that both the dual-KF method and the joint-EKF method can provide the better estimates of the system state and the pump power simultaneously than the KF method while the dual-KF method is better on average than the joint-EKF method for given system parameters for the proposed case.The superiority of the dual-KF method can be attributed to the fact that the joint-EKF method is sub-optimal since linearization is used for approximation.

B. Performance improvement regarding varying parameters
In this section, we present the mean RPIs regarding the following varying parameters: • the measurement efficiency T , • the diffusion coefficient g, • the tendency constant c.The mean RPIs and the corresponding SEMs are calculated over N = 1 × 10 3 trials.
In Fig. 3, T varies from 0 to 1 while the other parameters are given as in Section IV-A.It can be seen that the mean RPI of both dual-KF method (red line) and joint-EKF method (blue line) has a positive relationship with the measurement efficiency.This is reasonable since a higher measurement efficiency means more information about the system is used for the estimation of and a better estimated yields better estimated state vector x.
In Fig. 4, the diffusion coefficient g increases from 0.005 to 0.035 while the other parameters are given as in Section IV-A.With the increasing of g, the mean RPI of c increases from around 0 up to around 48% (dual-KF) and 38% (joint-EKF) while the mean RPI of q c and p c increase from around 0 up to around 51% (dual-KF), 43% (joint-EKF) and 40% (dual-KF), 36% (joint-EKF), respectively.Thus, we see that the proposed algorithms show more improvement when the Fig. 4: RPI of (a) c , (b) q c and (c) p c with the diffusion coefficient g increasing from 0.005 to 0.028.diffusion coefficient is larger.The reason is that a larger noise ratio results in a larger deviation of to its tendency constant c and the superiority of the proposed methods is clearer.An extreme situation is that the performance of proposed methods is the same with the KF method when the diffusion coefficient is zero.In Fig. 5, the tendency constant c increases from 0.3 to 0.7 while g = 0.025 and the other parameters are given as in Section IV-A.With the increasing of c, the mean RPI of c increases from 10% to 33% for dual-KF and 6% to 25% for joint-EKF while the mean RPI of q c and p c increase from 10% to 27% for dual-KF of q c , from 6% to 23% for joint-EKF of q c , from 4% to 23% for dual-KF of p c , from 3% to 20% for joint-EKF of p c , respectively.The results show that there is a positivity relationship between the tendency constant c and the RPI although the increase of RPI is small compared with the results in Fig. 3 and Fig. 4. Thus, we see that the RPI is more determined by the measurement efficiency and the diffusion coefficient than the tendency constant.

V. CONCLUSION
In this paper, we considered the state and parameter estimation problem of an OPO system where the pump power is subject to a stochastic process.We first formulated the system dynamics in the state-space representation.Then, we provided analysis on quantum filtering and its classical analog under the assumption that the quantum system is linear and Gaussian.Details of restrictions on obtaining a classical analog for the quantum system were provided.Thus, the simultaneous estimation problem of both the state and the unknown parameter can be solved by using the dual-KF method and the joint-EKF method.For the dual-KF, the estimation of the state and the pump power is decoupled and two Kalman filters run concurrently.For the joint-EKF method, the state vector is augmented with the unknown pump power.Therefore, the augmented system is nonlinear and an EKF is then employed.The simulation results show that both these methods can achieve good estimation of the system.An average improvement of 30%-50% can be reached compared with Kalman filter method for the studied case.The simulation results also show that the performance of the proposed algorithms increases with increasing of parameters T , g and c.This work can also be extended to the state-parameter estimation of other systems.

APPENDIX
The true evolution of the system state x T is approximated by the conditioned state x c obtained by using a complete measurement record and known .We consider Homodyne measurement to all the three outputs âm , âlb and âlc .Thus, measurement of âm is The complete measurement is ŷ = ( ŷθ m , ŷθ lb , ŷθ lc ).For simulation, we let θ m = θ lb = θ lc .The corresponding measurement matrices are given in (45).Then, the system state can be estimated using (20) with the coefficient matrix given in (45) while the innovation is dw = ŷdt −C x dt.

Fig. 5 :
Fig.5: RPI of (a) c , (b) q c and (c) p c with the tendency constant c increasing from 0.3 to 0.7.