{'fraction_non_alphanumeric': 0.034024779425567864, 'fraction_numerical': 0.01060634503472874, 'mean_word_length': 4.372919818456883, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective Strategies for Identifying Model Parameters for Open Quantum Systems

The problem of identifiability of model parameters for open quantum systems is considered by investigating two-level dephasing systems. We discuss under which conditions full information about the Hamiltonian and dephasing parameters can be obtained. Using simulated experiments several different strategies for extracting model parameters from limited and noisy data are compared.


I. INTRODUCTION
Control and optimization of quantum systems have been recognized as important issues for many years [1] and control theory for quantum systems has been developed since the 1980s [2][3][4]. There has been considerable recent progress in both theory and experiment [5,6]. However, despite this progress, there are still many challenges. Most quantum control schemes rely on open-loop control design based on mathematical models of the system to be controlled. However, accurate models are not often not available, especially for manufactured quantum systems such as artificial quantum dot atoms or molecules. Therefore, system identification [7] is a crucial prerequisite for quantum control.
In the quantum information domain, procedures for characterization of quantum dynamical maps are often known as quantum-process tomography (QPT) [8][9][10] and many schemes have been proposed to identify the unitary (or completely positive) processes, for example, standard quantum-process tomography (SQPT) [11][12][13], ancilla-assisted process tomography (AAPT) [14][15][16] and direct characterization of quantum dynamics (DCQD) [17]. However, if control of the system's dynamics is the objective, what we really need to characterize is not a global process but the generators of the dynamical evolution such as the Hamiltonian and dissipation operators. The problem of Hamiltonian tomography (HT), though less well-understood, has also begun to be considered recently by a few authors [18][19][20][21]. Although QPT and HT differ in various regards, both try to infer information about the quantum dynamics from experiments performed on systems, and both can be studied from the point of view of system identification with broad tasks including (1) experimental design and data gathering, (2) choice of model sets and model calculation, and (3) model validation.
Recently the quantum system identification problem has been briefly explored from cybernetical point of view, and underlining the important role of experimental design [22]. In this article we follow this line of inquiry. Throughout the paper, we make the following basic as-sumptions: (1) the quantum system can be repeatedly initialized in a (limited) set of known states; (2) that we can let the system evolve for a desired time t; and (3) that some projective measurements can be performed on the quantum system. The main question we are interested in in this context is how the choice of the initialization and measurement affect the amount of information we can acquire about the dynamics of the system. Given any a limited range of options for the experimental design, e.g., a range of measurements we could perform, different choices for the initial states, or different control Hamiltonians, how to choose the best experimental design, and what are the theoretical limitations? Finally, we are interested in efficient ways to extracting the relevant information from noisy experimental data.
The paper is organized as follows: In Sec. II we discuss the model and basic design assumptions. Sec III deals with the general question of model identifiability in various settings, and in Sec IV we compare several different stategies for parameter estimation from a limited set of noisy data from simulated experiments see how they measure up.

II. MODEL AND DESIGN ASSUMPTIONS
To keep the analysis tractable we consider a simple model of a qubit subject to a Hamiltonian H and a system-bath interaction modelled by a single Lindblad operator V , i.e., with system dynamics governed by the master equation where the Lindbladian dissipation term is given by We shall further simplify the problem by assuming that V is a Hermitian operator representing a depolarizing channel or pure phase relaxation in some basis. Without loss of generality we can choose the basis so that V is diagonal, in fact we can choose V = The control Hamiltonian can be expanded with respect to the Pauli basis {σ x , σ y , σ z } with possibly time-dependent coefficients ω α (t). It is convenient to consider a real representation of the system. Following the approach in [23] we expand ρ with respect to the standard Pauli basis for the 2 × 2 Hermitian matrices where the coefficients are v α = Tr(ρσ α ). Similarly expanding the dynamical operators allows us to recast Eq. (1) in following Bloch equation Using this simple model for illustration we subsequently consider the experimental design from three aspects: (1) initialization procedures, (2) measurement choice and (3) Hamiltonian design.
(1) Initialization. We assume the ability to prepare the system in some initial state with respect to the basis {|0 , |1 }, which coincides with the eigenbasis of V . We can formally represent the initialization procedure by the operator Π(θ I ), which is the projector onto the state |ψ I , with I indicating initialization. With these restrictions the design of the initialization procedure is reduced to the selection of parameter θ I . Note that we assume that we can only prepare one fixed initial state, not a full set of basis states.
(2) Measurement. We assume the ability to perform a two-outcome projective measurement where the measurement basis states can be written as so that the choice of the measurement can be reduced to suitable choice of the parameter θ M , and we shall indicate this by writing M (θ M ).
(3) Hamiltonian. In practice we may or may not have the freedom to choose the type of Hamiltonian but it will be instructive to consider the identification problem for the following three cases: (a) ω z (t) ≡ ω z > 0 and ω x (t) = ω y (t) ≡ 0.
In this context the experimental design problem can be reduced to the problem of choosing suitable values θ I and θ M to identify the system parameters γ and ω * in each of the three cases above. The problem considered here is similar to that considered in [21], in particular we still only allow a single initial state and single fixed measurement, but unlike in [21] the initial state and the measurement are not assumed to commute with the dephasing operator.

III. MODEL IDENTIFIABILITY
We now consider the identification problem for the dephasing qubit system, concentrating on experimental design issues for simultaneously discriminating dephasing parameter and Hamiltonian parameter of two-level dephasing systems.
In this special case H commutes with the dephasing operator V . Solving the equation with the initial state (7) gives and applying the binary-outcome projective measure- from which we can see that we can obtain information about the system parameters ω z and γ if and only if sin θ I = 0 and sin θ M = 0, i.e., if neither the initial state preparation Π(θ I ) nor the measurement M commutes with H and V . The measurement traces also yield information about θ I and θ M , i.e., we can determine the relative angles between the initialization and measurement axis and the fixed Hamiltonian/dephasing axis, if they are not known a-priori. We also see that the visibility is maximized if sin θ I sin θ M = 1, which will be the case if the initialization and measurement axis are orthogonal to the joint Hamiltonian and dephasing axis.
These results make physical sense. As [H, V ] = 0, if the initial state preparation Π(θ I ) commutes with H and V then the initial state is a stationary state of the dynamics and the measurement outcome is constant in time c ± = 1 2 (1 ± 1). If sin θ I = 0 then the initial state is not stationary and the state follows a spiral path towards the joint Hamiltonian and dephasing axis but as both H and V are proportional to σ z , Tr[σ z ρ(t)] is a conserved quantity of the dynamics. If sin θ M = 0 then the measurement commutes with σ z , and as Tr[σ z ρ(t)] is a conserved quantity, we again obtain no information about the dynamics from which to identify the system parameters.
In this case, the equation (6) is reduced into the following equation and the solution of Eq. (13) for the initial state (7) is where we set ω If ω 2 x < γ 2 /4 then ω will be purely imaginary and the sine and cosine terms above are replaced by the respective hyperbolic functions. If ω 2 x = γ 2 /4, the expression ω −1 sin( ωt) must be analytically continued.
Eq. (14) implies that we can obtain the estimated value of the following probabilities depends on both ω x and γ this shows that we can obtain full information about the system parameters if and only if cos θ I = 0 and cos θ M = 0. If cos θ I = 0 or cos θ M = 0 then we can still identify γ but not ω x . Again this makes sense.
cos θ I = 0 for θ I = π 2 , i.e., if the initial state is an eigenstate of the Hamiltonian. Since [H, V ] = 0 in this case, eigenstates of H are not stationary. However, since the Hamiltonian and dephasing axis are orthogonal, the initial state remains in a plane orthogonal to the dephasing axis, the z = 0 plane in our case, following the path x(t) = e −γt . Thus we have [H, ρ(t)] = 0 for all times, and we can therefore not obtain any information about the Hamiltonian parameter ω x , but we can still obtain information about the dephasing parameter γ. If the Hamiltonian and dephasing axis were not orthogonal then we would expect to be able to identify both the Hamiltonian and dephasing parameters even if the initial state was an eigenstate of H as in this case it would not remain an eigenstate of H under the evolution.
If cos θ I = 0 but cos θ M = 0 then the measurement commutes with the Hamiltonian. Transforming to the Heisenberg picture, , and again one can show that M ± (t) remains orthogonal to the dephasing axis and Tr[M ± (t)ρ 0 ] = Tr[M ± ρ(t)] is independent of the Hamiltonian H, explaining why we cannot obtain any information about H in this case.
In this case equation (6) is reduced into the following equation whose solution for the initial state (7) is where ω y = ω 2 y − γ 2 4 and Φ y This implies that we can obtain the estimated value of the probabilities where the coefficient functions are As before, if ω 2 y < γ 2 /4 then ω will be purely imaginary and the sine and cosine terms above turn into their respective hyperbolic sine and cosine equivalents, and if ω 2 y = γ 2 /4, the expression ω −1 sin( ωt) must be analytically continued.
In this case it is quite interesting to notice that it is impossible to find such θ I and θ M that cos(θ I − θ M ) = sin(θ I − θ M ) = cos(θ I + θ M ) = 0 and thus that we can identify both model parameters for any choice of the initial state and measurement. This also makes sense because regardless of the choice of θ I and θ M the initial

IV. PARAMETER ESTIMATION
In this section we explore how to estimate dephasing and Hamiltonian parameters from limited noisy measurement data. In an actual experiment we can only estimate the probabilities p ± (t) at a finite number of times t k by repeatedly initializing the system in some fixed state ρ 0 , and letting it evolve for time t k before performing the projective measurement M . Each single repetition of the experiment yields a binary outcome +1 or −1 and we can estimate the probability p ± (t k ) by repeating the experiment N e times and computing the relative frequencies of the respective measurement outcomes ±1, e.g., p ± (t k ) = N± Ne .
To model noisy experimental data we could generate p ± (t) by adding a zero-mean White-Gaussian noise signal g t to p ± (t), i.e., p ± (t) = p ± (t) + g t . By the Law of Large Numbers and Iterated-logarithm Law [24] this gives a Gaussian distribution p ± (t) with mean p ± (t) and variance σ 2 ∼ log log Ne 2Ne for N e → ∞. For N e large this should be a good error model, but for small N e it may not accurately capture the nature of the projection noise, which follows a Poisson distribution. To more accurately model noisy experimental data when the number of measurement repetitions is relatively small we can simulate the actual experiment by generating N e random numbers r n between 0 and 1, drawn from a uniform distribution, and setting p + = N 0 /N e , where N 0 is the number of r n ≤ p + .

A. Fourier and time-series analysis
One commonly used technique to find frequency components in a noisy time-domain signal is the Fourier transform or its discrete version, the discrete or Fast Fourier Transform (FFT). In principle, this allows us to estimate the Hamiltonian parameters (frequencies) ω z , ω x andω y [20], and the damping rate γ can be estimated from the Lorentzian broadening of the Fourier peak [21]. When the signal is sparse, i.e., we only have a relatively small number of sample points, and noisy, however, this approach becomes problematic. Fig. 3 shows that we can still identify a peak in the spectrum for sparse noisy signals but the accuracy of the estimate is limited, and accurately estimating the dephasing rate γ from the broadening of the peak, given the distortion, is very challenging.
Alternatively, once an estimated value for ω z has been obtained, γ can be deduced in other ways. In the simplest case A, where the evolution is given by Eq.
Givenp + (t), if θ I , θ M and ω z are known, we can therefore in principle calculate z(t) and thusγ. However, in practice the problem is that p + (t) and all the other parameters are not known precisely. In this case we could estimate γ using the Least Square Method [25] and the principle of Sequential Analysis [26,27].
To explore this approach we generated simulated noisy data signals p + (t) with standard deviation σ t according to the procedure above, assuming θ I = θ M = π 2 , ω z = 50, γ = 50 and the sampled period T s = 1.4µs. If the number of measurements N e at each time t is N e = 100 then the standard deviation of estimated error is log log Ne 2Ne = 0.0874. The original noiseless probabilityp + (t), the noisy signal with standard deviation 0.0874 and the power spectral density are shown in Fig.4. Again the frequency ω z can be easily determined by the peak value in Fig. 4(right). We attempt to determineγ based on p + (t) using simulated time-series. For a single time series p + (t) the distribution ofγ around the true value of 50 is shown in Fig. 5(left). The estimation is unsatisfactory even if the signal is sampled densely at a high time resolution. The results can be noticeably improved by averaging over multiple time-series of p + (t). If the mean value of the γ-estimates over several time series is used as the final estimate for γ, the estimated error ofγ can be remarkably reduced as shown in Fig. 5 (right), which indicates the mean value ofγ with 1 − 40 time-series respectively.
The number of simulated time-series of p + (t) required is determined by the desired estimation accuracy. To achieve a mean value of the distribution of ca. 50 with standard deviation less than 1.0 simulations suggest that at least 5 time-series with N e = 100 are necessary, i.e., about 500 measurements have to be performed at each sample time t to ensure that the standard deviation of the estimated error is less than 1.0. Instead of five time series with N e = 100 we could choose four times series with N e = 125 or two with N e = 250 to reach the nearly same estimation accuracy, as shown in Fig. 6. In the actual experiment one can obtain the mean-value sequence ofγ by continually updating p(t) at each sampled time t. If the mean-value sequence is found to converge to a certain fixed value and its standard deviation satisfies the requirement, the experiment measurement should be stopped.
To estimate γ by finite measurement data one will have to identify a suitable time t d when to terminate the experiment, ideally when the true value of e −γt d is almost 0. A p + (t) is almost equal to p − (t) for t ≥ t d , the signal after this time will be effectively a pure noise signal. Therefore increasing the signal length beyond a certain critical time t k0 will not improve the estimation accuracy.
Furthermore, we consider the effect of different ω z and γ on the estimation accuracy. If the true value of γ is less than the value of 50 above, the estimation of ω z and γ will be more accurate (as shown in Fig. 7) as the exponential e −γt decays more slowly. If the value of ω z is reduced, i.e., the oscillation period is increased while its envelope remains unchanged the accuracy of the estimation is not influenced (as shown in Fig. 8).

B. Bayesian Analysis of Sparse Signals
The previous section shows that we can in principle simultaneously identify both dephasing parameter and Hamiltonian parameters for simple open systems using Fourier and time-series analysis. Both Fourier and timeseries analysis can deal with very noisy signals but they require rather dense sampling of the signal with is both time-consuming and resource-intensive. Ideally we would like to be able to estimate the parameters by sparse sampling the signal. One promising approach in this regard is Bayesian estimation. The main idea behind the Bayesian approach is to choose the parameters to be determined, here ω and γ, to maximize a certain likelihood function We can still identify a peak in the Fourier spectrum (right) around ω = 1 but the accuracy of this estimate is limited, and estimating γ from the peak broadening would be a challenge. where d is the measured data and the p are the probabilities predicted by the model, which depend on the parameters to be determined, and σ is the error variance.
Following the same approach as in [19], we write the signals as a linear combination a small number of basis functions determined by the functional form of the signals. Here the measurement signalsp ± (t) can be written as a linear combination of two basis functions p ± (t) = α 1 g 1 (t) + α 2 g 2 (t).
The values of the basis functions and the coefficients are given in Table I   all cases, even if the signal is very noisy, the sampling is sparse and the measurement and initialization choices are not ideal to maximize the visibility or signal-to-noise ratio such as θ M = π 4 and θ I = π 3 . The squeezed peaks indicate that the accuracy of the frequency estimation is  much higher than the accuracy of the γ-estimates.
An additional advantage of the Bayesian estimation does not require a-priori knowledge of the initialization or measurement angles θ I and θ M . Rather, the estimation procedure provides values for the coefficients of the basis functions, which are related to the parameters θ I and θ M .

V. CONCLUDING DISCUSSION
We have studied the issue of model identifiability and experiment design for open system dynamics for a dephasing qubit. From the examples in Sec. III we can derive some general insights about the limits of identifiability and the role of experiment design. Unlike for process tomography, where we require the ability to prepare the system in many different input states, and the ability to measure a complete set of observables, we can in general extract information about both the Hamiltonian and the dephasing parameters by repeating a single experiment, initializing the system in single fixed initial state and measuring a fixed basis. There are certain limitations, however. We gain no information about the value of the system parameters if the initial state is a stationary Power spectral density Frequency (Hz) wz=50Hz,r=50 wz=20Hz,r=50 wz=50Hz,r=20 wz=20Hz,r=20 Standard deviation of noise is 0.0874 FIG. 7: Power spectral density of measured signal in different cases state of the system, or if the measurement is a conserved quantity, although knowledge of the stationary states or conserved quantities restricts the dynamics and thus provides indirect information about the system. Even if the initial state is not stationary or the measurement is not a conserved quantity, we may fail to obtain information about the Hamiltonian parameters is if the operators H and V are orthogonal and M commutes with H, for instance. These limitations also apply to higher dimensional systems, although for such systems additional restrictions on the identifiability of model parameters may arise.
If the experiment design is such that the model parameters are identifiable there are various ways to extract the relevant parameters from a set of noisy samples (time series), including Fourier analysis, time series analysis and Bayesian estimation. Although all of these approaches are in principle able to provide the required information, Bayesian estimation appears to be superior to the alternatives, in particular when we are dealing with a limited number of noisy data points. One reason the Bayesian estimation is capable of providing far more accurate estimates for the parameters given the same input data is that it utilizes information about the structure of the signal in the form of the choice of the basis functions we are projecting onto. This allows us to overcome restrictions on the uncertainties of the parameter estimates imposed by Nyquist's law.