Estimation of a general time-dependent Hamiltonian for a single qubit

The Hamiltonian of a closed quantum system governs its complete time evolution. While Hamiltonians with time-variation in a single basis can be recovered using a variety of methods, for more general Hamiltonians the presence of non-commuting terms complicates the reconstruction. Here using a single trapped ion, we propose and experimentally demonstrate a method for estimating a time-dependent Hamiltonian of a single qubit. We measure the time evolution of the qubit in a fixed basis as a function of a time-independent offset term added to the Hamiltonian. The initially unknown Hamiltonian arises from transporting an ion through a static laser beam. Hamiltonian estimation allows us to estimate the spatial beam intensity profile and the ion velocity as a function of time. The estimation technique is general enough that it can be applied to other quantum systems, aiding the pursuit of high-operational fidelities in quantum control.

E stimation of the underlying dynamics, which drive the evolution of systems is a key problem in many areas of physics and engineering. This knowledge allows control inputs to be designed, which account for imperfections in the physical implementation. For closed quantum systems, the time dependence of a system is driven by the Hamiltonian through Schrödinger's equation. If the Hamiltonian is static in time, a wide range of techniques have been proposed for recovering the Hamiltonian 1-4 , which have been applied to a variety of systems including chemical processes 5 and quantum dots 6,7 . These methods often involve estimation of the eigenvectors and eigenvalues of the Hamiltonian via spectroscopy, or through pulse-probe techniques for which a Fourier transform of the time evolution gives information about the spectrum.
These methods are not directly applicable to time-dependent Hamiltonians, which are becoming increasingly important as quantum engineering pursues a combination of high-operational fidelities and speeds, often involving fast variation of control fields, which are particularly susceptible to distortion before reaching the quantum device [8][9][10][11][12] . The time-varying case has thus far been studied in cases where the variation is along a single dimension in the Hilbert space, which for the commonly studied spin is a single spatial direction. In the case that the measured fields dominate the evolution (strong field limit), measurement of the system evolution as a function of time suffices for the reconstruction. For fields which are weaker than other available control fields (weak-field limit) the latter can be used to modulate the effect of the signal Hamiltonian on the quantum system [13][14][15] , providing an excellent signal-to-noise ratio. A further complication arises when a time-varying Hamiltonian contains non-commuting terms (for example, time-variation along two spin axes), because the evolution of the quantum system depends not only on their separate influences, but also on products arising from the noncommutativity. For unspecified time-dependent coefficients, no analytical solution to Schrödinger's equation exists 16,17 . In the weak-field limit, strong control fields can be used to separate out the different components using modulation, however, when the Hamiltonian itself is strong (as is the case in fast quantum control) these techniques cannot be applied.
In this article, we propose and demonstrate a method for reconstructing a general time-dependent single qubit Hamiltonian with non-commuting terms. The technique involves observing the evolution of the spin projection on the z-axis, while applying a static offset to one of the terms of the Hamiltonian. By varying the static offset, we build up data sets, which contain sufficient information to extract the full time-dependent Hamiltonian. Parameterizing the two time-dependent terms using basis splines (B-splines), we introduce an iterative fitting technique, which finds the Hamiltonian that best matches the data. We benchmark the reconstruction method experimentally by transporting a single trapped ion through a static laser beam, a technique suited to scaling up trapped-ion quantum information processing 18,19 . We perform two consistency checks on the Hamiltonian estimation using four separate reconstructions. For the first two, we compare two cases which use the same ion velocity profile, but different laser beam positions. For the second consistency check, we use the same laser beam, but change the velocity profile between the two by using different sets of timevarying control potentials. The method produces consistent experimental parameters in both cases, indicating the success of the reconstruction technique. Our method is applicable to spin Hamiltonians of the general formĤ¼ P i f i t ð Þŝ i , where the f i (t) are arbitrary time-dependent functions andŝ i are the Pauli operators.

Results
Hamiltonian estimation method. In our experiments, a Hamiltonian with two non-commuting time-dependent terms arises when we perform quantum logic gates by transporting an ion through a static laser beam 18,19 . In this case, the Hamiltonian describing the interaction between the ion and the laser can be written in an appropriate rotating frame aŝ which includes a time-varying Rabi frequency O(t), and an effective detuning d(t), which is related to the first-order Doppler shift of the laser in the rest frame of the moving ion (see Methods for details).
To reconstruct the Hamiltonian, we make use of two additional capabilities. First, we can switch-off the Hamiltonian at time t off on a timescale much faster than the qubit evolution. Subsequently measuring in theŝ z basis, we can obtainŝ z t off ð Þ h i. On its own, this does not allow us to separate the contributions from O(t) and d(t). To do this, we use a second capability, which is the ability to add a controlled offsetĤ s ¼' d Lŝz =2 to the Hamiltonian, resulting inĤ I t ð Þ þĤ s . The resulting spin measurement is now dependent on both t off and the set value of d L . Repeating the experiment for a range of values of d L but with otherwise identical settings, we obtain an estimate of the expectation value which we denote aŝ s meas z t off ; d L ð Þ . Hamiltonian extraction involves finding the functions d(t) and O(t), which generate spin populationsŝ sim z t off ; d L ð Þ that most closely match the data. We minimize the reduced w 2 cost function.
where n ¼ N À n À 1 is the number of degrees of freedom, with N the number of data points, n the number of fitting parameters and s meas (t off , d L ) the s.e. on the estimated hŝ meas Þi. This is subject to the initial condition C t ¼ 0; d L ð Þ j i ¼0 j i, and the following restrictions, which are imposed by quantum mechanics for all d L .
One challenge in obtaining an estimate for the Hamiltonian is that we must optimize over continuous functions d(t) and O(t).
To address this, we represent d(t) and O(t) with a linear combination of B-spline polynomials, which allow the construction of smooth functions using only a few parameters 20 . Any smooth function S(t) can be written in terms of B-spline polynomials B i,k (t) and a set of weights a i as The polynomial B-spline functions B i,k (t) are of order k with each polynomial centred at a time t i , which is parameterized by the index i. Further details and an example can be found in Methods.
Using the B-spline form for d(t) and O(t), the cost function is minimized by adjusting the weights of the B-spline decomposition. Solving this optimization problem in general is hard, because it is non-linear and non-convex due to the nature of Schrödinger's equation and the use of projective measurements. This produces a non-trivial relation between the weights and the spin populations as discussed in Methods. To overcome this challenge, we have implemented a method which we call 'Extending the Horizon Estimation' (EHE) in analogy to a well-established technique called 'Moving Horizon Estimation' (MHE) 21 .
The key idea is that because our measurement data arises from a causal evolution, we can also estimate the Hamiltonian in a causal way. Instead of optimizing J over the complete time span at once, we first restrict ourselves to a small, initial time span reaching only up to the start of the qubit dynamics 0ot off oT 0 , where we denote T 0 as the time horizon for the first step. Optimizing J over this short time span requires fewer optimization parameters and is simpler than attempting to optimize over the full data set. Once we have solved this small sub-problem, we extend the time horizon to T 1 where T 1 ¼T 0 þ t and re-run the optimization, extrapolating the results of the initial time span into the extended region to provide good starting conditions for the subsequent optimization. This procedure is iterated until the time span extends over the whole data set T n max ¼ max t off ð Þ ð Þ . The method allows us to reduce the number of B-spline functions used to represent d(t) and O(t), and also reduces the amount of data considered in the early stages of the fit, when the least is known about the parameters. This facilitates the use of non-linear minimization routines, which are based on local linearization of the problem and converge faster near the optimum. More details regarding the optimization routine can be found in Methods.
Conceptually EHE is very similar to MHE. The main difference is that in MHE the time span has a fixed length and thus its origin is shifted forward in time along with the horizon. In EHE, the origin stays fixed at the expense of having to increase the time span under consideration. MHE avoids this by introducing a so-called arrival cost to approximate the previous costs incurred before the start of the time span. This keeps the computational burden fixed over time, which is very important as MHE is usually used to estimate the state of a system in real-time, often on severely constrained embedded platforms. Since neither constraint applies to our problem, we decided to extend the horizon rather than finding an approximate arrival cost. This is advantageous since finding the arrival cost in the general case is still an open problem 22 . Due to the similarity between MHE and EHE, we anticipate future improvements by adapting techniques used in MHE to EHE. This might be used to reduce the data-processing required for reconstruction, which for EHE scales as T 2 n max .
Experimental implementation. To test the ability of the method to reliably extract a Hamiltonian from data, we apply it to the Hamiltonian for an trapped-ion qubit during transport through a near-resonant laser beam. Our qubit is encoded in the electronic states of a calcium ion, which is defined by . This transition is well-resolved from all other transitions, and has an optical frequency of o 0 /(2p)C411.0420 THz. The laser beam points at 45°to the transport axis, and has an approximately Gaussian spatial intensity distribution. The time-dependent velocity _ z t ð Þ of the ion is controlled by adiabatic translation of the potential well in which the ion is trapped. This is implemented by applying time-varying potentials to multiple electrodes of a segmented ion trap, which are generated using a multi-channel arbitrary-waveform generator, each output of which is connected to a pair of electrodes via a passive third-order low-pass Butterworth filter. The result is that the ion experiences a time-varying Rabi frequency O(t) and a laser phase which varies with time as where f(z(t)) ¼ k z (z(t))z(t) with k z (z(t)) the laser wave vector projected onto the transport axis at position z(t) and o L the laser frequency. The spatial variation of k z (z(t)) accounts for the curvature of the wavefronts of the Gaussian laser beam. To create a Hamiltonian of the form of equation (1), we work with the differential of the phase, which gives a detuning the laser detuning from resonance. For planar wavefronts, k 0 z z ð Þ¼0, and d(t) corresponds to the familiar expression for the first-order Doppler shift (see Methods for details).
The experimental sequence is depicted in Fig. 1. We start in zone B by cooling all motional modes of the ion to " no3 using a combination of Doppler and electromagnetically induced transparency cooling 23 , and then initialize the internal state by optical pumping into 0 j i. The ion is then transported to zone A, and the laser beam used to implement the Hamiltonian is turned on in zone B. The ion is then transported through this laser beam to zone C. During the passage through the laser beam, we rapidly turn the beam off at time t off and thus stop the qubit dynamics. The ion is then returned to the central zone B to perform state readout, which measures the qubit in the computational basisŝ z (for more details see Methods). The additional HamiltonianĤ s is implemented by offsetting the laser frequency used in the experiment by a detuning d L . For each setting of t off and d L the experiment is repeated 100 times, allowing us to obtain an estimate for the qubit populationsŝ meas We first perform a comparison in which the ion velocity is the same but the beam position is changed. Thus we expect to obtain two different profiles for O(t) but the same velocity profile, which is closely related to d(t). Experimental data is shown in Fig. 2 alongside the results of fitting performed using our iterative method. The beam positions used for each data set differ by B64 mm along the transport axis, but the transport waveform used was identical. It can be seen from the residuals that the estimation is able to find a Hamiltonian, which results in a close match to the data. The estimated coefficients of the Hamiltonian extracted from the two data sets are shown in Fig. 3a,b. To estimate the relevant errors of our reconstruction, we have performed non-parametric resampling with replacement, optimizing for the solution using the same set of B-spline functions as was used for the experimental data to provide a new estimate for the Hamiltonian. This is repeated for a large number of samples, resulting in a distribution for the estimated values of d(t) and O(t) from which we extract statistical properties such as the s.e. The error bounds shown in Fig. 3 correspond to the s.e. on the mean obtained from these distributions (see Methods for further details). It can be seen that the values of d(t) for the two different beam positions have a similar form but a fixed offset for the region where the reconstructions overlap. We believe that this effect arises from the non-planar wavefronts of the laser beam. Inverting the expression for d(t) to obtain the velocity of the ion, we find Using this correction, we find that the two velocity profiles agree if we assume that the ion passes through the centre of the beam at a distance of 2.27 mm before the minimum beam waist, a value which is consistent with experimental uncertainties due to beam propagation and possible mis-positioning of the ion trap with respect to the fixed final focusing lens. The velocity estimates taking account of this effect are shown in Fig. 3c.
Our second comparison involves using two different velocity profiles but with a common beam position. The resolution in both time and detuning were lower in this case than for the data shown in Fig. 2 (see Methods for the data). Figure 4 shows the results of the reconstruction. We observe that the estimated Rabi frequency profiles agree to within the error bars of the reconstruction. One interesting feature of this plot is that the error bars produced from the resampled data sets increase near the peak. We believe that this happens because the sampling time of the data is 0.5 ms, which starts to become comparable to the Rabi frequency (the Nyquist frequency is 1 MHz). To optimize the efficiency of our method, it would be advantageous to run the reconstruction method in parallel with data taking, thus allowing updating of the sampling time and frequency resolution based on the current estimates of parameter values.

Discussion
Our method for directly obtaining a non-commuting timedependent Hamiltonian uses straightforward measurements of the qubit state in a fixed basis as a function of time and a controlled offset to the Hamiltonian. Unlike schemes based on dynamical modulation or continuous strong driving, it avoids the need for control fields which act more strongly on the qubit than the Hamiltonian to be measured. This is a key advantage in quantum technologies where the Hamiltonian of interest is often already close to the limit of system drive strength. A processtomography-based approach would require that for every time step multiple input states be introduced, and a measurement made in multiple bases [24][25][26] . This requires a much greater level of control than the method presented above. An effective modulation of the measurement basis arises in our approach due to the additional detuning d L . Nevertheless, it is also worth noting that tomography provides more information than our method: it makes no assumptions about the dynamics aside from that of a completely positive map while we require coherent dynamics. Extensions to our work are required to provide a rigorous estimation of the efficiency of the method in terms of the precision obtained for a given number of measurements, and to see whether a similar approach could be taken for non-unitary dynamics. We have recently used these methods to improve the control over the ion velocity, which is of direct value in optimizing transport operations in scalable trapped-ion quantum information processing 11,12,27 , and will be essential for realizing multi-qubit transport gates 18 . We expect them to be applicable across a wide range of physical systems where such control is available, including those considered for quantum computation 4,6,7,28-31 .

Methods
Derivation of Hamiltonian. The interaction of a laser beam with frequency o L and wave vector k(z(t)) with a two-level atom with resonant frequency o 0 and time-dependent ion position z(t) ¼ (0, 0, z(t)) can be described in the Schrödinger picture by the Hamiltonian where the Rabi frequency O(z(t)) gives the interaction strength between the laser and the two energy levels. We can define the laser phase at the position of the ion as )z(t) and k z (z(t)) ¼ |k|cos(y(t)) being the projection of the laser beam onto the z-axis along which the ion is transported. Here y(t) is the angle between the wave vector k(z(t)) and the transport axis evaluated at position z(t). Moving to a rotating frame using the unitary transformation U¼e À i F t ð Þ 2 and applying the rotating wave approximation with respect to optical frequencies, we obtain Defining a static detuning d L ¼ o L À o 0 , we obtain which is the expression used in the main text.
B-spline curves and optimization algorithm. The set of polynomial B-spline functions B i,k (t) of order k are recursively defined over the index i over a set of points K ¼ {t 0 , t 1 , .., t n þ k }, which is referred to as the knot vector 20 . Figure 5 gives a visualization of the B-splines B i,k (t) and a B-spline curve. The B-spline construction ensures that any linear combination of the B-splines is continuous and has (k À 2) continuous derivatives. The knot vector K determines how the basis functions are positioned within the interval [t 0 , t n þ k ]. We notice that for our Hamiltonian the spacing of the B-splines is not critical, which we think is due to the smoothness of the variations in our Hamiltonian parameters d(t) and O(t). We therefore used the Matlab function spap2 to automatically choose a suitable knot vector and restricted ourselves to optimizing the coefficients a i . We collect all coefficients a i for d(t) and O(t) and store them in a single vector a.
A detailed algorithmic summary of our implementation of the EHE method is given below.
1. Searching for a starting point: here we reconstruct the Hamiltonian for a first, minimal time horizon such that we can then use this as a starting point to iteratively extend the horizon as described in step 2.
(a) Choose an initial time horizon such that it contains the region where the first discernible qubit dynamics occur. (b) Cut down the number of fitting parameters as much as possible, for example, by using few B-splines of low order. This amounts to choosing empirically a low number of B-splines (and thus the length of a 0 ), which might represent d(t) and O(t) over the given region. (c) Use a non-linear least-squares fitting routine to minimize J by varying the parameters a 0 . In the case that the initial fit is poor or no minimum is found, try new initial conditions, change the number of B-spline functions, or manually adjust the function using prior knowledge of the physical system under consideration.

T n+1
Time This procedure is used to provide a starting point for the optimization over the initially chosen window, which is typically performed with a set of higher order B-splines. From this starting point, we iteratively extend the fitting method to the full data set as follows.
2. Extend the horizon: this step is repeated until the whole time horizon is covered.
It consists of the following sequence, which is illustrated in Fig. 5.
(a) Extend the time horizon by t from T n to T n þ 1 ¼T n þ t. If all these fail, we have to resort to increasing the bound on w 2 red .
3. Post-processing: the following steps are optional and were performed manually in cases where we wished to improve the fit or examine its behaviour.
(a) The optimization over the whole time horizon was re-run using different numbers of B-splines for d(t) and O(t). This was used to check the sensitivity of the fit. (b) The optimization over the whole time horizon was re-run using a starting point based on the previously found optimum plus randomized deviations. This tested the robustness of the final fit.
Wavefront correction. For plane waves, we find that _ f t ð Þ¼k Á v t ð Þ, which is the well-known expression for the first-order Doppler shift. For transport through a real Gaussian beam, the wave vector direction changes with position. Taking this into account, the derivative of f(t) becomes where k 0 z ¼ dk z =dz and _ z t ð Þ is the component of the ion velocity along the z-axis. We extract d(t) using our Hamiltonian estimation procedure, thus to obtain the velocity of the ion we use As the ion moves through the beam it experiences the same magnitude of the wave vector |k| ¼ 2p/l, but the angle y between the ion direction and the wave vector changes. Written as a function of this angle, the velocity becomes where y 0 (z(t)) ¼ dy(z(t))/dz(t). We parameterize our Gaussian beam according to Fig. 6. The phase is given as a function of both the position along the beam axis x and the perpendicular distance from this axis k by 32 where the Gaussian beam parameters include the beam waist W(x), the radius of curvature R(x), the Rayleigh range x R and the Guoy phase shift z(x). These are given by the expressions where W 0 is the minimum beam waist and l the laser wavelength. The ion moves along the z-axis as shown in Fig. 6. In the kx-plane a unit vector e l (k, x) perpendicular to the wavefronts is given by and the unit vector e v pointing along the direction of transport is given by The angle y(x) between the wave and position vector is then given by the dot product which can be written in terms of the full set of parameters above as where in our experiments a ¼ 3p/4. Using equations (12) and (18), we examined the value of x cl required for the velocity to match for our two beam positions. We find that they agree for x cl ¼ À 2.27 mm, which is within the experimental uncertainties for our set-up.