Quantum Hamiltonian Identification with Classical Colored Measurement Noise

TL;DR

This paper introduces a method for identifying quantum Hamiltonians in closed systems affected by classical colored measurement noise, using an augmented system model and eigenstate realization to accurately estimate unknown parameters.

Contribution

It presents a novel Hamiltonian identification approach that accounts for classical colored measurement noise via an augmented system model and transfer-function techniques.

Findings

Effective identification of quantum Hamiltonians with colored noise

Verification on a two-qubit system demonstrates accuracy

Method improves robustness against measurement disturbances

Abstract

In this paper, we present a Hamiltonian identification method for a closed quantum system whose time trace observables are measured with colored measurement noise. The dynamics of the quantum system are described by a Liouville equation which can be converted to a coherence vector representation. Since the measurement process is disturbed by classical colored noise, we introduce an augmented system model to describe the total dynamics, where the classical colored noise is parameterized. Based on the augmented system model as well as the measurement data, we can find a realization of the quantum system with unknown parameters by employing an Eigenstate Realization Algorithm. The unknown parameters can be identified using a transfer-function-based technique. An example of a two-qubit system with colored measurement noise is demonstrated to verify the effectiveness of our method.

Tables

Table 1. TABLE I: Comparison of the solutions obtained by the ZS and our methods

	${𝝎}_{\mathrm{𝟏}}\left(\text{GHz}\right)$	${𝝎}_{\mathrm{𝟐}}\left(\text{GHz}\right)$	${𝜹}_{\mathrm{𝟏}}\left(\text{GHz}\right)$
real values	$1.3$	$2.4$	$4.3$
ZS method in [13]	$0.8746$	$2.0601$	$4.1211$
our method	$1.3$	$2.4$	$4.3$

Equations

$$\mathbf{H}=\sum _{m=1}^{N^2-1}{a}_m{\mathbf{X}}_m,$$

$$[{\mathbf{X}}_j,{\mathbf{X}}_k]=\sum _{l=1}^{N^2-1}{C}_{jkl}{\mathbf{X}}_l,$$

$$\dot{\mathbf{\rho }}(t)=-i[\mathbf{H},\mathbf{\rho }(t)].$$

$$\begin{array}{rl}& {\dot{x}}_j(t)=i\sum _{l=1}^{N^2-1}(\sum _{m=1}^{N^2-1}{a}_m{C}_{mjl})x_l(t),\\ & j=1,2,\cdots ,N^2-1.\end{array}$$

$$\bm{y}(t)=\left[\begin{array}[]{cccc}\bm{y}_{1}(t),&\bm{y}_{2}(t),&\cdots,&\bm{y}_{L}(t)\\ \end{array}\right]^{T}\\$$

$$\begin{array}{rl}& \dot{\mathbf{x}}(t)=\mathbf{Ax}(t),\mathbf{x}(0)={\mathbf{x}}_0\\ & \mathbf{y}(t)=\mathbf{Cx}(t),\\ & \mathbf{A}\in {\mathbb{R}}^{K\times K},\mathbf{C}\in {\mathbb{R}}^{L\times K},\end{array}$$

$$S(\omega )={\int }_{-\infty }^{\infty }{e}^{-i\omega \tau }R(\tau )\mathrm{d}\tau .$$

$$S(\omega )=\Gamma (s){\Gamma }^T(-s){|}_{s=i\omega },$$

$$\Gamma (s)=\frac{{\beta }_1{s}^{n-1}+{\beta }_2{s}^{n-2}+\cdots +{\beta }_{n-1}s+{\beta }_n}{s^n+{\alpha }_1{s}^{n-1}+\cdots +{\alpha }_{n-1}s+{\alpha }_n}$$

$$\begin{array}{rl}& \dot{\mathbf{\xi }}(t)=\mathbf{E}\mathbf{\xi }(t)+\mathbf{F}\eta (t),\\ & v(t)=\mathbf{G}\mathbf{\xi }(t),\end{array}$$

$$\begin{split}&\bm{E}=\left[\begin{array}[]{ccccc}0&1&0&\cdots&0\\ 0&0&1&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&1\\ -\alpha_{n}&-\alpha_{n-1}&-\alpha_{n-2}&\cdots&-\alpha_{1}\end{array}\right],\\ &\bm{F}=\left[0\quad 0\quad\cdots\quad 0\quad 1\right]^{T},\\ &\bm{G}=[\beta_{n}\quad\beta_{n-1}\quad\beta_{n-2}\quad\cdots\quad\beta_{1}],\end{split}$$

$$\begin{array}{rl}& \dot{\bar{\mathbf{\xi }}}(t)=\mathbf{E}\bar{\mathbf{\xi }}(t),\bar{\mathbf{\xi }}(0)={\bar{\mathbf{\xi }}}_0,\\ & \bar{v}(t)=\mathbf{G}\bar{\mathbf{\xi }}(t),\end{array}$$

$${\tilde{\mathbf{y}}}_i(t)={\mathbf{y}}_i(t)+\bar{v}(t),$$

$$\begin{split}\bm{\tilde{y}}(t)=&\left[\begin{array}[]{c}\bm{y}_{1}(t)+\bar{v}(t)\\ \bm{y}_{2}(t)+\bar{v}(t)\\ \vdots\\ \bm{y}_{L}(t)+\bar{v}(t)\end{array}\right]=\bm{Cx}(t)+\left[\begin{array}[]{c}1\\ 1\\ \vdots\\ 1\end{array}\right]_{L\times 1}\bar{v}(t).\end{split}$$

$$\begin{array}{rl}& \dot{\check{\mathbf{x}}}(t)=\check{\mathbf{A}}\check{\mathbf{x}}(t),\check{\mathbf{x}}(0)={\check{\mathbf{x}}}_0,\\ & \tilde{\mathbf{y}}(t)=\check{\mathbf{C}}\check{\mathbf{x}}(t),\end{array}$$

$$\check{\mathbf{A}}$$

$$\check{\mathbf{C}}$$

$${\mathbf{G}}_L$$

$$\check{\mathbf{x}}(k+1)$$

$$\tilde{\mathbf{y}}(k)$$

$$\tilde{\mathbf{y}}(k)=\check{\mathbf{C}}{\check{\mathbf{A}}}_d^k{\check{\mathbf{x}}}_0.$$

$$\small\begin{split}&\bm{H}_{rs}(k)=\\ &\left[\begin{array}[]{cccc}\bm{\tilde{y}}(k)&\bm{\tilde{y}}(k+t_{1})&\cdots&\bm{\tilde{y}}(k+t_{s-1})\\ \bm{\tilde{y}}(k+j_{1})&\bm{\tilde{y}}(k+j_{1}+t_{2})&\cdots&\bm{\tilde{y}}(k+j_{1}+t_{s-1})\\ \vdots&\vdots&&\vdots\\ \bm{\tilde{y}}(k+j_{r-1})&\bm{\tilde{y}}(k+j_{r-1}+t_{2})&\cdots&\bm{\tilde{y}}(k+j_{r-1}+t_{s-1})\end{array}\right]\end{split}$$

$$\bm{H}_{rs}(0)=\bm{P}\left[\begin{array}[]{cc}\bm{D}&\bm{0}\\ \bm{0}&\bm{0}\end{array}\right]\bm{Q}^{T}=[\bm{P}_{1}\quad\bm{P}_{2}]\left[\begin{array}[]{cc}\bm{D}&\bm{0}\\ \bm{0}&\bm{0}\end{array}\right]\left[\begin{array}[]{c}\bm{Q_{1}}^{T}\\ \bm{Q}_{2}^{T}\end{array}\right],$$

$$\begin{array}{rl}& \hat{\mathbf{x}}(k+1)={\hat{\mathbf{A}}}_d\hat{\mathbf{x}}(k),\hat{\mathbf{x}}(0)={\hat{\mathbf{x}}}_0,\\ & \tilde{\mathbf{y}}(k)=\hat{\mathbf{C}}\hat{\mathbf{x}}(k),\end{array}$$

$$\check{\mathbf{C}}(s{\mathbf{I}}_{\check{n}}-\check{\mathbf{A}}{)}^{-1}{\check{\mathbf{x}}}_0=\hat{\mathbf{C}}(s{\mathbf{I}}_{\check{n}}-\hat{\mathbf{A}}{)}^{-1}{\hat{\mathbf{x}}}_0.$$

$$\begin{split}&P(s)=\det(s\bm{I}_{\check{n}}-\bm{\check{A}})\\ &Q(s)=\det\big{(}s\left[\begin{array}[]{cc}\bm{I}_{\check{n}}&\bm{0}\\ \bm{0}&\bm{0}\end{array}\right]-\left[\begin{array}[]{cc}\bm{\check{A}}&\bm{\check{x}}_{0}\\ \bm{\check{C}}&\bm{0}\end{array}\right]\big{)}.\end{split}$$

$$\mathbf{H}=\sum _{\alpha =1}^2\frac{{\omega }_{\alpha }}{2}{\mathbf{\sigma }}_z^{\alpha }+{\delta }_1({\mathbf{\sigma }}_{+}^1{\mathbf{\sigma }}_{-}^2+{\mathbf{\sigma }}_{-}^1{\mathbf{\sigma }}_{+}^2),$$

$$\bm{\sigma}_{z}=\left[\begin{array}[]{cc}1&0\\ 0&-1\end{array}\right],\bm{\sigma}_{x}=\left[\begin{array}[]{cc}0&1\\ 1&0\end{array}\right],\bm{\sigma}_{y}=\left[\begin{array}[]{cc}0&-i\\ i&0\end{array}\right],$$

$$\bm{\sigma}_{+}=\left[\begin{array}[]{cc}0&1\\ 0&0\end{array}\right],\bm{\sigma}_{-}=\left[\begin{array}[]{cc}0&0\\ 1&0\end{array}\right].$$

$$\displaystyle\left[\begin{array}[]{c}\dot{x}_{1}(t)\\ \dot{x}_{2}(t)\\ \dot{x}_{3}(t)\\ \dot{x}_{4}(t)\end{array}\right]$$

$$y(t)$$

$$S(\omega )=\frac{10^{12}{\omega }^2+4\times 10^{26}}{{\omega }^4-3.999\times 10^{13}{\omega }^2+4\times 10^{26}}.$$

Full text

I Introduction

In recent years, great progress has been achieved on quantum technologies, such as quantum computation [1], quantum communication [2], and quantum metrology [3]. Relying on the precise models of relevant quantum systems, these quantum-mechanics-based techniques can achieve better performance than classical-mechanics-based counterparts. However, under some circumstances, some parameters in these models may not be known, which result in degraded performance. For acquiring these parameters, a fundamental step is identification. Classical system identification theory uses some methods, such as wavelet cross spectrum analysis, least squares methods, or maximum-likelihood estimators [4], to estimate the system parameters utilizing the input and output data of the system.

Inspired by classical system identification theory, research on quantum parameter identification is undergoing developments, among which Hamiltonian identification has been widely studied. Since the Hamiltonian of a quantum system determines the evolution of the quantum state, quantum Hamiltonian identification becomes an important research area, where the identifiability for checking the existence and uniqueness of solutions of quantum Hamiltonian identification is a core issue. An identifiability condition was firstly established for a closed quantum system driven by a laser field, where the populations of all states are measured [5]. By relaxing the above measurement condition; i.e., measuring limited observables, an identifiability condition for the same problem was found in [6]. Alternatively, for a closed many-body quantum system, Sone and Cappellaro  [7] analyzed the identifiability for the Hamiltonian identification under a single quantum probe and tested the estimation performance in the presence of Gaussian noises. Also, the identifiability for linear passive quantum systems coupling to quantum bosonic fields was developed in [8].

Studies on identifiability for quantum Hamiltonian identification provide theoretical analysis for designing identification algorithms which can be divided into two categories; i.e., time-domain and frequency-domain approaches. In the time-domain approaches, a class of observer-based methods was proposed. Kosut and Rabitz [9] used an invariant asymptotic state observer to estimate parameters in Hamiltonian with a gradient algorithm. Afterwards, an adaptive observer of exponential convergence was proposed by Bonnabel and Mirrahimi [10] to directly estimate the parameters with Gaussian measurement noises and control. Moreover, quantum tomography technique was applied to quantum Hamiltonian identification. In [11], a two-step identification algorithm was presented for a closed quantum system, which is based on the framework of quantum process tomography. Jagadish and Shaji [12] also proposed an algorithm to identify the coupling Hamiltonian between a qubit and its environments using measured data from quantum process tomography. Although the identification algorithm is usually effective, the tomography process is time-consuming. In addition, a system-realization-based method was proposed to identify unknown parameters in the Hamiltonian of a linear passive quantum system [8].

On the other hand, the frequency-domain approaches were explored for quantum Hamiltonian identification. Zhang and Sarovar [13] estimated parameters in the Hamiltonian of spin systems using equivalent transfer functions where time trace observables are measured. This method was later extended to open quantum systems [14]. Fourier analysis was also applied to quantum Hamiltonian identification. Cole et. al. [15] adopted Fourier transform of the measurement of one observable to identify Hamiltonian of a closed two-level quantum system. Schirmer et. al. [16] also estimated the Hamiltonian of a two-level quantum system based on Fourier analysis and Burgarth et. al. [17] provided the solution for the $N$-level case. In the previous works, the measurement process is assumed to be ideal or disturbed by Gaussian noise, which is convenient for the analysis and design of identification algorithms. However, in practice, we use classical devices in the measurement process which may carry classical colored noise [18, 19]. Ignoring the impact of colored noise may degrade the performance of the identification algorithms. Nevertheless, it is still an open problem on quantum Hamiltonian identification under the condition that the measurement process is disturbed by classical colored noise.

In this paper, we propose an augmented system method to identify parameters in the Hamiltonian of a closed quantum system where the data of measured time trace observables carry classical colored noise. We consider a closed finite-level quantum system which can be represented in a coherence vector representation. Correspondingly, when the observable of the system is specified, the dynamics of the system can be described by a reduced equation for the coherence vector. To combine the colored noise into the model for identification, a spectral factorization method is utilized such that an augmented system model for quantum Hamiltonian identification can be obtained. Moreover, with the data of the time trace observables, an Eigenstate Realization Algorithm is employed to find a realization of the augmented system. Equalling the transfer function of both the original system with unknown parameters and the realization generated by the measurement data, we can obtain a set of nonlinear algebraic equations for the unknown parameters, which is difficult to be solved analytically. Numerically, these equations can be solved using a PHCpack [20]. Finally, we provide an example of a two-qubit system with measurement process disturbed by classical colored noise.

The remainder of this paper is organized as follows. Section II describes the model for identification of finite-level quantum systems. In Section III, we develop a colored noise realization which is utilized to augment the original system model. The procedure to obtain identified Hamiltonian is presented in Section IV. In Section V, the effectiveness of our algorithm is verified in an example of a two-qubit system with classical colored measurement noise. The conclusions are drawn in Section VI.

II Identification model for finite-level quantum system

In this paper, we consider a closed $N$ dimensional quantum system. The system Hamiltonian $\bm{H}$ satisfies $i\bm{H}\in\mathfrak{su}(N)$, where the Lie algebra $\mathfrak{su}(N)$ can be represented by $N\times N$ traceless skew-Hermitian matrices. The dimension of $\mathfrak{su}(N)$ over $\mathbb{R}$ is $N^{2}-1$. Therefore we can expand the Hamiltonian as

[TABLE]

where $\mathcal{X}=\{i\bm{X}_{m},m=1,2,\cdots,N^{2}-1\}$ is a set of orthogonal bases for $\mathfrak{su}(N)$ [13]. The commutation relations for the elements in $\mathcal{X}$ are

[TABLE]

where the commutator $[\cdot,\cdot]$ is calculated as $[\bm{X},\bm{Y}]=\bm{XY}-\bm{YX}$ for two operators $\bm{X}$ and $\bm{Y}$ and $C_{jkl}$ are antisymmetric constants with respect to the interchange of any pair of indices [21]. This property indicates that $C_{jkl}$ equals to zero with any two identical indices.

The Hamiltonian determines the dynamics of the density matrix $\bm{\rho}(t)$ of the system as

[TABLE]

which is the so-called Liouville equation [22]. The density matrix $\bm{\rho}$ describes the probability distribution of the system states, which is an $N\times N$ Hermitian and positive semi-definite matrix with $\text{tr}(\bm{\rho})=1$. We have assumed $\hbar=1$.

The Liouville equation (3) can be transformed into a coherence vector representation [23] which is convenient for the design of the identification algorithm. In this representation, the state of the system is alternatively described by a coherence vector $x=[x_{1},\cdots,x_{N^{2}-1}]^{T}$ where $x_{j}$ is the expectation value of $\bm{X}_{j}$; i.e., $x_{j}=\text{tr}(\bm{X}_{j}\bm{\rho})$. The corresponding dynamical equation is thus written as

[TABLE]

To observe the quantum system, we can choose $L$ observables ${\bm{O}_{1},\bm{O}_{2},\cdots,\bm{O}_{L}}$ and their expectations can be taken as the outputs of the system; i.e.,

[TABLE]

where $\bm{y}_{i}(t)=\langle\bm{O}_{i}\rangle={\rm tr}[\bm{O}_{i}\bm{\rho}]$. An observable $\bm{O}_{i}$ can be expanded in terms of the bases in $\mathcal{X}$ as $\bm{O}_{i}=\sum_{j}o_{j}^{(i)}\bm{X}_{j}$ where the corresponding bases $\mathcal{M}=\{\bm{X}_{\bm{\mu}_{1}},\bm{X}_{\bm{\mu}_{2}},\cdots,\bm{X}_{\bm{\mu}_{p}}\}$ span a minimal space containing the $L$ observables. Here, $\bm{\mu}_{i}$ denotes the indices for the corresponding bases and the number of $\bm{\mu}_{i}$ is $p$.

With respect to the bases in $\mathcal{M}$, an accessible set of $\mathcal{M}$ can be generated using a filtration process [13]. Denoting $\mathcal{F}_{0}=\mathcal{M}$, an iterative procedure can be calculated $\mathcal{F}_{i}=\mathcal{F}_{i-1}\cup[\mathcal{F}_{i-1},\mathcal{X}]$, where $[\mathcal{F}_{i-1},\mathcal{X}]=\{\bm{X}_{j}|\text{tr}(\bm{X}_{j}^{{\dagger}}[g,h])\neq 0,g\in\mathcal{F}_{i-1},h\in\mathcal{X}\}$, until $\mathcal{F}_{i}$ saturates. Supposing the final set is $\mathcal{F}=\{\bm{X}_{\bm{\mu}_{1}},\bm{X}_{\bm{\mu}_{2}},\cdots,\bm{X}_{\bm{\mu}_{K}}\}$ with a size $K$, the corresponding reduced coherence vector is $\bm{x}=[x_{\bm{\mu}_{1}},x_{\bm{\mu}_{2}},\cdots,x_{\bm{\mu}_{p}},\cdots,x_{\bm{\mu}_{K}}]^{T}$. Therefore, the corresponding reduced dynamical equation for the reduced coherence vector can be written as

[TABLE]

where $\bm{x}_{0}$ is the initial state, $\bm{A}_{jl}=-i\sum_{m=1}^{N^{2}-1}a_{m}C_{m\bm{\mu}_{j}\bm{\mu}_{l}}$ and $\bm{C}$ is configured such that $\bm{y}(t)$ are expectation values of our measured observables. The equation (5) affords the basic model for the identification problem.

III Augmented Model For Output Disturbed By Classical Colored Noise

In existing Hamiltonian identification studies, ideal case or Gaussian noise disturbing measurement results have been considered. However, in practice, the measurement results may be polluted by classical colored noise arising from measurement devices. Hence, it is necessary to introduce the classical colored noise into the dynamics and thus we obtain a complete model for the purpose of Hamiltonian identification.

III-A Classical Colored Noise Model

In general, classical colored noise can be characterized by a shaped power spectral density (PSD) $S(\omega)$ describing the signal power distribution over all the frequency components [24]. Since the PSD and autocorrelation $R(t)$ form a Fourier transform pair, the PSD can be obtained by

[TABLE]

The spectral factorization theorem tells that a positive, rational and strictly proper PSD $S(\omega)$ can be factorized as

[TABLE]

where $\Gamma(s)$ is a causal transfer function which results from the internal dynamics of noises. Here, $\Gamma(s)$ characterizes a mapping of a white noise input $\eta(s)=\mathcal{L}[\eta(t)]$ to a colored noise output $v(s)=\mathcal{L}[v(t)]$. The operator $\mathcal{L}[\cdot]$ is Laplace transform and $\eta(t)$ and $v(t)$ are the white noise input and the colored noise in the time domain, respectively. Note that the power spectral density of $v(t)$ is $S(\omega)$ which will reduce to a flat one when the output noise is white. Also, for an irrational power spectral density, we can find its rational approximants using Padé approximation [25] and then a transfer function can be obtained by the factorization.

For a given transfer function $\Gamma(s)$, it is easy to construct a corresponding minimal realization [26]. For the single-input-single-output (SISO) transfer function

[TABLE]

which is strictly proper, we can write its realization in a controllable canonical form as

[TABLE]

with

[TABLE]

where $\bm{\xi}\in\mathbb{R}^{n}$ can be considered as the internal state vector of the noise with an initial state $\bm{\xi}(0)=\bm{\xi}_{0}$. It is straightforward to obtain the dynamics of the expectation $\bm{\bar{\xi}}$ of the internal mode $\bm{{\xi}}$ as

[TABLE]

where $\bar{v}(t)$ is the expectation of the colored noise $v(t)$ and the input term vanishes since $\bar{\eta}(t)=0$.

III-B Augmented System Model

We consider that the classical colored noise $v(t)$ is additive and we assume that the $i$-th output $\bm{y}_{i}(t)$ of the quantum system; i.e., the time trace of the $i$-th observable $\bm{O}_{i}$, is disturbed by the expectation value of the classical colored noise. Hence, the polluted output can be expressed as

[TABLE]

where $\bm{y}_{i}(t)$ is the original quantum output for the $i$-th observable and $\tilde{\bm{y}}_{i}(t)$ is the polluted one. Hence, the polluted output of the system can be written as

[TABLE]

Further, taking the expression of $\bar{v}(t)$ in (11) into (13) and denoting a new state vector $\bm{\check{x}}$ as \bm{\check{x}}=\left[\begin{array}[]{cc}\bm{x}(t)^{T}&\bar{\bm{\xi}}(t)^{T}\end{array}\right]^{T}, we can combine the original system with the noise model as

[TABLE]

with

[TABLE]

where $\bm{\check{x}}\in\mathbb{R}^{\check{n}},\bm{\check{A}}\in\mathbb{R}^{\check{n}\times\check{n}},\bm{\check{C}}\in\mathbb{R}^{L\times\check{n}}$ and the initial state of the augmented state is \bm{\check{x}}_{0}=\left[\begin{array}[]{cc}\bm{x}_{0}^{T}&\bar{\bm{\xi}}_{0}^{T}\end{array}\right]^{T}. The order of the augmented model $\check{n}$ satisfies $\check{n}=K+n$ where $K$ and $n$ are the orders of the quantum system model and the colored noise realization, respectively.

Now, we obtain an augmented system (14) of the finite-level quantum system whose output is disturbed by classical colored noise. A similar model can be found in the design of Kalman filter under classical colored measurement noise [27]. Note that both the dynamics of the quantum system and the classical colored noise contribute to the polluted outputs $\bm{\tilde{y}}(t)$. However, it is difficult to distinguish the unpolluted quantum output $\bm{{y}}(t)$ from the noise directly. A possible method is in demand for extracting the quantum information of the original system from the polluted output.

IV Hamiltonian Identification with the Time Trace Observables Polluted by Classical Colored Noise

IV-A Problem Statement

In Section III, we have introduced the dynamics of the coherence vector for the finite-level quantum system and considered the measurement process disturbed by classical colored noise which is represented by a linear system realization for a given spectrum density $S(\omega)$. Consequently, we have obtained a parameterized augmented system model for describing the total dynamics of the coherence vector and the internal modes of the noise. We aim to identify the unknown coefficients in $\{a_{m}\in\mathbb{R},m=1,2,\cdots,M\}$ for a quantum system with the Hamiltonian (1) using the augmented model. Hence, we state our Hamiltonian identification problem as follows.

Given the structure of the augmented system model (14) for the closed quantum system with an initial state $\bm{\check{x}}_{0}$, our identification problem is to estimate the unknown parameters $\{a_{m}\in\mathbb{R},m=1,2,\cdots,M\}$ in the model utilizing the disturbed measurements $\bm{\tilde{y}}(t)$ of the time traces of the observables for the underlying system.

IV-B Measurement Process

To access the information of the quantum system, we measure the time trace observables (3). We assume that we have a large number of copies of the system prepared in an identical initial state. Thus we can measure an observable on many identical systems at a time and then obtain the expectation value of the observable. We call a set of identical systems we measure as an ensemble. We measure in this way because in quantum mechanics quantum measurement for any observable will generally change the system state. Moreover, since we cannot obtain measurement results for the same system at different time instants, we need to make measurements on different ensembles for different instants.

Concretely speaking, we sample the observables with an equal interval $\Delta t$. Denoting the measured value of $\bm{O}_{i}$ of the $j$-th copy at a time instant $k\Delta t$ as $y_{i}^{(j)}(k)$, the measured value of $\bm{O}_{i}$ accompanied by colored noise $v^{(j)}(k)$ can be expressed as $\tilde{y}_{i}^{(j)}(k)=y_{i}^{(j)}(k)+v^{(j)}(k)$. Here, we have written $\{y_{i}(k\Delta t)\}$ as $\{y_{i}(k)\}$ for simplicity. After measuring many copies at different time instants, we can average over the measurement results and obtain (12). The measurement process is shown in Fig. 1.

Note that when we consider multi-observable time traces, commutative observables can be measured simultaneously. However, due to the uncertainty principle [28], measurements of non-commutative observables should be carried out on different ensembles.

IV-C Identification Algorithm for Measurement Data with Classical Colored Noise

To link the discrete measurement data to the system model of the augmented system model (14), we discretize (14) for a given sampling interval $\Delta t$ and thus we have

[TABLE]

with $\bm{\check{x}}(0)=\bm{\check{x}}_{0}$, where $\bm{\check{A}}_{d}=e^{\bm{\check{A}}\Delta t}\in\mathbb{R}^{\check{n}\times\check{n}},$ $\bm{\check{C}}\in\mathbb{R}^{L\times\check{n}}$, and $k$ stands for the $k$-th time step. Hence, the initial state response of (IV-C) can be written as

[TABLE]

It is difficult to solve $\{a_{m}\}$ from (26) since it is transcendental in $\{a_{m}\}$ and the system dimension of $\check{n}$ is unknown. However, these difficulties can be overcome by combining a system realization method [13] and an Eigenstate Realization algorithm (ERA) [29]. Using the ERA, the dimension $\check{n}$ can be determined. In addition, a system realization can be constructed based on the measurement data such that the unknown parameters can be obtained by solving a set of nonlinear equations arising from the the equivalence between the transfer functions of the original system and the realization. It should be mentioned that we calculate the realization using the measurement data disturbed by classical colored noise naturally which results from both the dynamics of quantum system and the noise.

The ERA begins with a generalized Hankel matrix

[TABLE]

which is constructed by the measurement result. Its dimension is $rL\times s$ with two integers $r$ and $s$. To determine $\check{n}$ accurately, it is good to choose a sufficient number of measurement results; i.e., two large integers $r$ and $s$ are preferred.

We choose the measurement results from the initial time; i.e., $k=0$, and then we can have a singular value decomposition of $\bm{H}_{rs}(0)$ as

[TABLE]

where $\bm{P}\in\mathbb{R}^{rL\times rL},\bm{Q}\in\mathbb{R}^{s\times s}$ are unitary matrices, and they are partitioned into $\bm{P}_{1},\bm{P}_{2}$, and $\bm{Q}_{1},\bm{Q}_{2}$ with respect to the dimension of the diagonal square matrix of $\bm{D}$. Since the diagonal elements of $\bm{D}$ are positive singular values of $\bm{H}_{rs}(0)$, its dimension is determined by the number of the singular values.

Using the identity matrix $\bm{I}$ with a subscript indicating its dimension, we define matrix $\bm{E}^{T}_{L}=[\bm{I}_{L},\bm{0}_{L},\cdots,\bm{0}_{L}]_{L\times rL}$ and $\bm{e}_{1}$ is the first column of $\bm{I}_{s}$. According to [29], we can establish a numerical realization

[TABLE]

with $\bm{\hat{A}}_{d}=\bm{D}^{-1/2}\bm{P}_{1}^{T}\bm{H}_{rs}(1)\bm{Q}_{1}\bm{D}^{-1/2}$, $\bm{\hat{C}}=\bm{E}^{T}_{L}\bm{P}_{1}\bm{D}^{1/2}$, and $\bm{\hat{x}}_{0}=\bm{D}^{1/2}\bm{Q}_{1}^{T}\bm{e}_{1}$ [29]. It is clear that the order of this realization is exactly the dimension of $\bm{D}$. Further, because of the equivalence between system models (IV-C) and (29), the dimension $\check{n}$ of uncertain model (IV-C) is expected to be the same as that of (29) or $\bm{D}$. Then letting $\bm{\hat{A}}=\log\bm{\hat{A}}_{d}/\Delta t$, the pair $(\bm{\hat{A}},\bm{\hat{C}},\bm{\hat{x}}_{0})$ formulates a continuous-time realization describing dynamics of both quantum system and colored noise.

Thus far, we have completed the process of developing a numerical time-continuous realization from measured data and also determined the system dimension $\check{n}$. Moreover, we have built up an augmented model $(\bm{\check{A}},\bm{\check{C}},\bm{\check{x}}_{0})$ in (14). The corresponding transfer functions from the initial states to the outputs of the two models should be equal [26]; that is,

[TABLE]

The left-hand side contains the parameters to be estimated, while the right-hand side is completely determined by measured data. To solve $\{a_{m}\}$ from this equation, we just need to equal the coefficients of $s$ in all orders of the both sides. Firstly, the left-hand side can be simplified as $Q(s)/P(s)$ [26], where

[TABLE]

In fact, the coefficients of $s$ in different orders in $P(s),Q(s)$ are polynomials of the unknown parameters $\{a_{m}\}$, and the corresponding coefficients in the right hand are numbers obtained through experiment. Next, what we need to do is to solve these polynomial equations for unknown Hamiltonian parameters $\{a_{m}\}$. Since these equations are often nonlinear or in high order, professional numerical tools such as PHCpack [20] can be used to obtain final results. To this end, we have introduced the whole procedure to identify unknown parameters in the Hamiltonian under disturbed measurements with classical colored noise . This procedure can be applied to finite-dimensional closed quantum system with measurements disturbed by classical colored noise as long as the noise can be represented by linear system models though spectral factorization of its PSD.

V An Example for A Two-qubit System

In this section, we consider a two-qubit system whose Hamiltonian is written as

[TABLE]

with Pauli matrices

[TABLE]

and the corresponding ladder operators

[TABLE]

Here, the superscripts in the Hamiltonian label the qubits. For the identification aim, we measure the local observable $\bm{\sigma}_{x}^{1}$ of the first qubit and thus with the observable-induced accessible set and the Hamiltonian, a dynamical equation of the coherence vector of the two-qubit system can be written as

[TABLE]

where $x_{1}(t)=\langle\bm{\sigma}_{x}^{1}(t)\rangle$, $x_{2}(t)=\langle\bm{\sigma}_{y}^{1}(t)\rangle$, $x_{3}(t)=\langle\bm{\sigma}_{z}^{1}\bm{\sigma}_{x}^{2}(t)\rangle$ and $x_{4}(t)=\langle\bm{\sigma}_{z}^{1}\bm{\sigma}_{y}^{2}(t)\rangle$. The initial state is set as [x_{1}(0)\quad x_{2}(0)\quad x_{3}(0)\quad x_{4}(0)]^{T}=\left[\begin{array}[]{cccc}0&1&0&0\\ \end{array}\right]^{T}. To simulate the real dynamics of the quantum system, we set the real parameters as $\omega_{1}=1.3~{}\text{GHz},\omega_{2}=2.4~{}\text{GHz}$ and $\delta_{1}=4.3~{}\text{GHz}$.

As for the classical colored noise $v(t)$ added in measurement process, we assume that its power spectrum density is expressed as

[TABLE]

Factorizing the spectrum $S(\omega)$, we can obtain a transfer function $\Gamma(s)$ as below

[TABLE]

A realization in a controllable canonical form can be found as

[TABLE]

with a two dimensional internal mode $\bm{\xi}(t)$.

We verify the validity of this noise realization (56) by checking its power spectral density. First, imposing a white noise signal $\eta(t)$ on both the transfer function (55) and the noise realization (56) with an arbitrary initial state, we obtain two output colored noise signals. Then we estimate the PSDs of the two signals using Welch’s overlapped segment averaging estimator [30]. The two estimated PSD curves are compared with the theoretical PSD (54) as shown in Fig. 2. It can be seen that the three curves are of a similar tendency. However, we witness a disparity due to the imperfect white noise and the error for estimating the PSD. Therefore, in simulation, we can produce the expectation of colored noises using the realization (56) with an arbitrary initial state.

The parameters to be identified are $\omega_{1},\omega_{2}$ and $\delta_{1}$ as well as the matrices and the state in (11) for the classical colored noise. Before we analyze the measurement data, we may not know the order $n$ for the colored noise model. However, with the noise model, we can still write an augmented system model for the parameter identification as

[TABLE]

with

[TABLE]

where the augmented state vector is $\bm{\check{x}}(t)=[x_{1}(t),x_{2}(t),x_{3}(t),\bar{\bm{\xi}}(t)^{T}]^{T}$. Note that for the classical colored noise, there exist many realizations which are equivalent and related by a similarity transformation. Therefore, we can just fix the part elements of the matrices in the noise realization to reduce the parameters to be identified. In our simulation, we assume the output vector $\bm{G}=[1,1,\cdots,1]_{1\times n}$.

Moreover, with the sampling time $\Delta t=0.1\mu s$ and the final time $T=12\mu s$, the measurement and the real results of the output $\langle\sigma_{x}^{1}\rangle$ are plotted as the black-dot and yellow-dot lines, respectively, in Fig. 3, where the polluted measurement result has a discrepancy of its real value. However, the measurement result $\{\tilde{y}(k)\}$ can still be used to construct the Hankel matrix $\bm{H}_{rs}(0)_{rL\times s}$ where we let $r=20,L=1,s=100$. Consequently, we can obtain the corresponding singular value decomposition, where we plot the singular values of $\bm{H}_{rs}(0)$ in logarithmic scale in Fig. 4. We can easily find a huge gap between the dominant singular values and the other quite small ones. Therefore, we can determine the dimension of the augmented system (57) according to the number of the dominant singular values and thus the order of the colored noise realization can be determined. Hence, we have the dimension of the augmented system $\check{n}=6$ and the order of the noise realization $n=2$.

Following the procedure in Section IV of the ERA, a realization $(\bm{\hat{A}},\bm{\hat{C}},\bm{\hat{x}}(0))$ with the dimension $\check{n}$ can be identical to that of $(\bm{\check{A}},\bm{\check{C}},\bm{\check{x}}_{0})$ (57); i.e.,

[TABLE]

whose both sides are polynomials in $s$. Specifically, the left-hand side of (59) contains the parameters to be identified while the right-hand side of (59) is numerically constructed by the measurement data. Equalling the coefficients of $s$ in the same order on the both sides, we can collect a polynomial equation set. In our example, the nine lowest order polynomial equations containing nine unknown variables are

[TABLE]

where $e_{ij}$ is the element of $\bm{E}$ in the $i$-th row and $j$-th column and $\xi_{0i}$ is the $i$-th element of $\bm{\xi}_{0}$. A set of solutions can be obtained utilizing a PHCpack [20].

In Table I, we compare our identified result with that obtained using Zhang and Sarovar’s method (ZS method) [13] which is not specially designed for the colored noise case. Our method can precisely identify the real values of the parameters in the Hamiltonian if the number of copies is not limited. While the results obtained by the method in [13] are different from the real values.

Moreover, in our method the estimates for the noise realization are $e_{11}=10^{6}(-50.025+410.92i),e_{12}=10^{6}(-49.97+415.41i),e_{21}=10^{6}(49.94-406.43i),e_{22}=10^{6}(49.925-410.92i),\xi_{01}=10^{6}(-4053.8+3692.7i),\xi_{02}=10^{6}(4053.5-410.92i)$. After substituting the solution of our method back into the augmented system model (57) and the solution of ZS method back into the quantum system model (47), the outputs produced by the two identified systems are compared with the real measurements in Fig. 5. It can be seen that our identification results coincide with the real system which shows that our method can improve the identification accuracy with classical colored measurement noise.

VI Conclusion

In this paper, we designed a procedure to identify the Hamiltonian of closed quantum system under classical colored measurement noise and showed its performance on a two-qubit example. In our method, an augmented system model was constructed and ERA was used to eliminate the impact of classical colored measurement noise on the precision of Hamiltonian identification. In principle, our identification procedure is applicable to various colored noise in the measurement process and does not require any prior information (e.g., PSD) about noise. The future research will focus on extending our method to open quantum systems for Hamiltonian identification when there exists classical colored measurement noise.