PCCP View Article Online

PAPER

View Journal | View Issue

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

Exploring the Hamiltonian inversion landscape Cite this: Phys. Chem. Chem. Phys., 2014, 16, 15615

Ashley Donovan and Herschel Rabitz* The identification of quantum system Hamiltonians through the use of experimental data remains an important research goal. Seeking a Hamiltonian that is consistent with experimental measurements constitutes an excursion over a Hamiltonian inversion landscape, which is the quality of reproducing the data as a function of the Hamiltonian parameters. Recent theoretical work showed that with sufficient experimental data there should be local convexity about the true Hamiltonian on the landscape. The present

Received 20th May 2014, Accepted 15th June 2014 DOI: 10.1039/c4cp02209b

paper builds on this result and performs simulations to test whether such convexity is observed. A gradientbased Hamiltonian search algorithm is incorporated into an inversion routine as a means to explore the local inversion landscape. The simulations consider idealized noise-free as well as noise-ridden experimental data. The results suggest that a sizable convex domain exists about the true Hamiltonian, even with a modest

www.rsc.org/pccp

amount of experimental data and in the presence of a reasonable level of noise.

1 Introduction Understanding and predicting the dynamics of a quantum system requires quantitative knowledge of the associated Hamiltonian. A large number of Hamiltonian inversion procedures have been developed utilizing experimental data, including the Rydberg– Klein–Rees method,1–5 functional sensitivity analysis,6,7 tomography procedures,8–11 Bayesian analysis,12 compressed sensing,13 and the use of measurement time traces.14 An optimal Hamiltonian identification procedure was also proposed based on quantum control experiments performed in a closed-loop fashion.15,16 In the latter approach, electromagnetic fields were tailored to shrink the distribution of Hamiltonians capable of reproducing the data within its error bounds. The searches for the optimal fields as well as the family of acceptable Hamiltonians were performed using genetic algorithms.17,18 Research is still needed to understand the fundamental aspects of quantum mechanical data inversion and the nature of efficient algorithms that can guide the experiments and perform the inversion in real time as the data is being generated. The search for Hamiltonians that successfully reproduce experimental data can be viewed as an excursion along an associated inversion landscape, which maps the Hamiltonian parameters to a data-dependent inversion cost function. As background, landscapes have been extensively studied in the context of quantum control by considering the observable as a function of the control variables. The topology of these quantum control landscapes is of prime interest, and under certain assumptions, these landscapes have been found to be devoid of suboptimal Department of Chemistry, Princeton University, Princeton, NJ 08544, USA. E-mail: [email protected]

This journal is © the Owner Societies 2014

extrema, or traps, thereby permitting unfettered discovery of optimal controls.19,20 Thus, quantum control landscapes are convex with an associated (level) set of optimal solutions,21 which is ideal for the purpose of finding at least one optimal control. Such favorable landscape topology has been numerically confirmed through the use of gradient-based search algorithms21–24 and in the laboratory.25,26 For the goal of data inversion, the highly nonlinear relationship between the Hamiltonian features and system observables makes it difficult to predict whether multiple Hamiltonians may exist that reproduce the experimental data within an acceptable tolerance level. However, theoretical work has shown that in the absence of noise and with access to a sufficient amount of experimental data, there should be local inversion landscape convexity about the true Hamiltonian.27,28 Importantly, when this circumstance occurs then the convexity is accompanied by a unique Hamiltonian, unlike the optimal control situation mentioned above where multiple solutions generally exist. For data inversion, this result implies that when starting with a trial Hamiltonian that lies within a suitable neighborhood of the true Hamiltonian, then it should be possible to use a local search algorithm to identify the true Hamiltonian. With this theoretical finding in mind, the present paper employs simulations to investigate the inversion landscape. Specifically, we seek to explore the boundary of the aforementioned local convexity through the use of a myopic Hamiltonian search algorithm that utilizes only a modest amount of quantum control experimental data. Fig. 1 provides a schematic overview of the Hamiltonian identification procedure utilized in this work. In the experimental component (A), a set of K control fields ek(t), k = 1,. . .,K are generated. Field noise, denoted Zek, may enter each experiment; here ek(t) + Zek is meant to accommodate any form of

Phys. Chem. Chem. Phys., 2014, 16, 15615--15622 | 15615

View Article Online

Paper

PCCP

experimental data sets are considered. These simulations aim to test whether minimization of the inversion cost function J leads to either level sets of multiple Hamiltonian solutions or convexity about the true Hamiltonian within a local neighborhood. Section 4 provides concluding remarks and discusses future directions for study.

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

2 Physical model and inversion algorithm Fig. 1 Overview of the Hamiltonian identification procedure. In the experimental component (A), a set of K fields (denoted ek(t) for k = 1,. . .,K) are chosen. A general form for field noise is denoted by Zek (e.g., additive or multiplicative noise can be incorporated through the field’s amplitudes and/or phases). The fields separately interact with the system through the Hamiltonian, H, and produce a set of measurements, denoted Ok. While multiple measurements could be taken for each field (say, over multiple time points t A [0,T]), in this work only Ok(T) is considered. Uncertainly in the measurement is accounted for with additive noise ZOk, though other noise models may also be considered. In the inversion component (B), a ˜ is chosen. The fields from A (i.e., the best trial set of Hamiltonians H estimates from the laboratory) are used to obtain computed measurements ˜ k(T). The goal of the inversion component is to search through the O ˆ k(T) within ˜ that reproduces the data O space of Hamiltonians to identify H experimental uncertainty; this is done through minimization of a cost 2 P ^k ðTÞ  O ~k ðTÞ over H ˜. function J ¼ ð1=KÞ O

2.1

The Hamiltonians H(t) considered in this work are of finite dimension N and assume the dipole approximation, H(t) = H0  me(t).



K   1X ~k ðTÞ 2 ; ^k ðTÞ  O O K k

(1)

whose minimization to an acceptable degree implies identification of a Hamiltonian from the data in (A). In the present work, the search for the Hamiltonian is performed using a gradientbased algorithm that is a variation of the D-MORPH technique.29 Importantly, as D-MORPH is a myopic search algorithm, it can provide insight into local inversion landscape features in a way that a genetic or stochastic algorithm could not. In particular, the breadth of the neighborhood around the true Hamiltonian that permits local convexity can be systematically explored. The paper is organized as follows. Section 2 lays out the physical model and develops the search algorithm used to identify elements of a system’s Hamiltonian. Section 3 presents numerical simulations that illustrate the procedure described by Fig. 1. There are many issues to consider in reaching a fully implementable inversion algorithm. The key issue explored here is an assessment of a practical level of convexity with a modest amount of data. Idealized noise-free as well as noise-ridden

15616 | Phys. Chem. Chem. Phys., 2014, 16, 15615--15622

(2)

In eqn (2), H0 is a diagonal field-free Hamiltonian matrix, m is the transition dipole matrix, and e(t) is the time-dependent applied field. The matrix structure of the Hamiltonian could arise from various discretization procedures; here the perspective considers N eigenstates of H0 as forming a basis. The observable used for the simulated experiments performed in Fig. 1A is the state-to-state transition probability Pi-f, Pi-f = |h f|U(T,0)|ii|2

k

noise (e.g., additive or multiplicative). The fields ek(t), k = 1,. . .,K separately interact with the system through the Hamiltonian H and produce a corresponding set of observables, denoted Ok(T), for measurement at terminal time T. Uncertainly in the kth measurement is accounted for by the noise ZOk. In the inversion ˜ is chosen. The input fields component (B), a trial Hamiltonian H ˜ ˜ , corresponding from (A) are used to compute Ok(T) based on H to the same type of measurement from (A). The inversion cost function J is defined,

Physical model

(3)

at time t = T where the unitary propagator U(t,0) satisfies the ¨dinger equation time-dependent Schro @ i h Uðt; 0Þ ¼ HðtÞUðt; 0Þ @t

(4)

for t A [0,T]. In the present inversion algorithm, the experiments in Fig. 1A are not performed to optimize a particular experimental outcome for the inversion objective. Rather, the experiments are treated as a means to collect a set of measurements that hopefully reflect the expected sensitive interaction between a field e(t) and a quantum system. In this regard, Pi-f provides a very demanding test for inversion purposes, as it just relies on data from the modulus of a single matrix element of U and at one time T. A more elaborate algorithm would seek particular measurements that enhance the sensitivity of the data and controls to the Hamiltonian. But, given the responsive nature of quantum dynamics to the Hamiltonian features, we will show that working with simple physically motivated controls and Pi-f already can function well for inversion purposes. The general goal of the Hamiltonian inversion algorithm in this paper is to identify the elements of the matrices H0 and/or m that reproduce the experimental Pi-f measurements. The search for these matrix elements is performed using the gradient-based D-MORPH algorithm presented below. 2.2

D-MORPH inversion algorithm

Background for the inversion algorithm is provided by a previous study considering optimization of Pi-f with the D-MORPH algorithm22 using Hamiltonian structural controls, expressed in terms of the elements of H0 and m. In the latter work, an arbitrary fixed field e(t) was introduced and the associated Hamiltonian

This journal is © the Owner Societies 2014

View Article Online

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

PCCP

Paper

structure landscape was found to be free of suboptimal traps under particular assumptions. Many Hamiltonian structures could maximize the state-to-state transition probability. In the context of the present paper, this previous work can be viewed as an inversion with just one observation, K = 1; the landscape is convex, but with a multitude of Hamiltonians producing the same single ‘data point’ Pi-f E 1.0. This paper explores whether even a modest increase of K and an associated set of controls ek(t) and observations P ki-f can produce a unique Hamiltonian. In an inversion procedure, both H0 and m may be sought for identification. The present study assumes prior sufficiently accurate knowledge of H0 (e.g., from an appropriate form of spectroscopy) and searches for the transition dipole matrix m. The transition dipole matrices are assumed to be real symmetric with diagonal elements set to zero. For clarity of presentation, the D-MORPH algorithm written below addresses dipole matrix elements. Parallel equations that address H0 elements can be readily derived,22 and a D-MORPH search may be performed to concurrently identify H0 and m. The D-MORPH search seeks to minimize the inversion cost function J¼

K  2 1X k k ~ P^i!f  P~i!f ðHÞ K k

(5)

˜ for a data set P ki-f, k = 1,. . .,K. over elements of m~ residing in H k k ˆ i-f = P i-f + Z represents an experimental measurement Here, P where the true value of P ki-f has additive noise Z (details are provided in Section 3). Full minimization of J would imply that ˜ ki-f (H ˜ ) values computed using the identified Hamiltonian all P ˜ equal the respective P ˆ ki-f values from the experiments within H the level of noise. A parameter s Z 0 is introduced to follow variations in the m~ matrix elements along an inversion landscape trajectory defined as m~(s). Minimizing the inversion cost function J (cf., eqn (5)) along a trajectory requires satisfying dJ X @J d~ mmn ¼ 0 (6) ds @~ m mn ds mn for all s. The inequality can be enforced by setting d~ mmn @J : ¼ @~ mmn ds

(7)

We have that  @ P~ k ð~ @J 2 X ^k i!f mÞ k Pi!f  P~i!f ¼ ð~ mÞ @~ mmn K k @~ mmn

! (8)

˜ ki-f /q~ mmn is:30 where qP k @ P~i!f

@~ mmn

ðT X 2 @~ m ¼  Im hqjki hkjU y ðt; 0Þ ek ðtÞUðt; 0Þjiidt h @~ m 0 mn kai

!

(9) †

with hq| = hi|U (T,0)| f i h f |U(T,0) and q~ m/q~ mmn is a matrix with 1 in the (m,n) and (n,m) positions and zero elsewhere. Eqn (7) and (8) assure a monotonic decrease in J until either J - 0 is approached, or the search stops at a suboptimal value of J; the latter outcome indicates that the landscape is not convex under

This journal is © the Owner Societies 2014

the posed conditions. Given a sufficient amount of data, a local convex neighborhood should exist wherein an initial ˜ converges to the true Hamiltonian H.27 HowHamiltonian H ˜ initially lying outside ever, the trajectory of a Hamiltonian H this domain may not converge to the true Hamiltonian H. As mentioned above, under specified control conditions for K = 1, global convexity should hold and at least one Hamiltonian should exist satisfying the data point Pi-f.22 But, the case of inversion as the goal rests on an analysis calling for a sufficiently large set of data K c 1 to assure local convexity of the inversion landscape about the true Hamiltonian as the unique solution when J B 0. The following numerical illustrations aim to explore these and other issues, including the impact of experimental noise on inversion quality.

3 Numerical illustrations The simulations presented below are meant to be an initial foray into understanding the Hamiltonian inversion landscape. The prior work27 considered the problem of identifying a nondegenerate quantum system’s real symmetric transition dipole matrix, where H0 is known and the diagonal elements of the dipole matrix are also known to be zero. The analytical results showed local inversion landscape convexity about the true Hamiltonian when (1) no noise was present, (2) a sufficiently large number of experimental measurements were utilized, (3) the trial Hamiltonians lie within some neighborhood around the true Hamiltonian, and (4) the time duration of the applied electromagnetic field, T, was arbitrarily large. These issues are selectively addressed below. The simulations also aim to illustrate the ability of a simple gradient search algorithm (D-MORPH) to correctly identify a quantum system’s Hamiltonian when a modest number of experiments with reasonable levels of noise-ridden data are incorporated into an inversion procedure. Any practical inversion procedure must rely only on the estimated variance of the identified Hamiltonian to assess the quality of the inversion, which is reported in the simulations here. As a check of the capability of the algorithm, the quality of the identification is also assessed in the simulations by comparing the inversion results to the true Hamiltonian (in this work, m matrix elements). All cases where J B 0 resulted in a small variance around the true Hamiltonian. Additionally, the present simulations utilize finite-dimensional systems and further assume knowledge of the system’s size (for example, only four-level Hamiltonians are considered when aiming to identify a four-dimensional system). Future studies will need to assess whether searching through a space of unknown dimension will significantly impact the inversion procedure. In this work, the state-to-state transition probability measurements are simulated using a set of applied fields. The fields have the form  ! L 8p T 2 X eðtÞ ¼ exp  2 t  al sinðol t þ fl Þ T 2 l

(10)

Phys. Chem. Chem. Phys., 2014, 16, 15615--15622 | 15617

View Article Online

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

Paper

PCCP

with L amplitudes, frequencies, and phases (al, ol, fl, respectively). All units are treated as dimensionless. In the simulations, T = 10 and the frequencies ol are chosen to be resonant with the system’s transitions between the known energy levels of H0. The field parameters are the sets of amplitudes and phases, which are often utilized in quantum control experiments.31–33 In each D-MORPH search, an initial trial transition dipole matrix m~ must be chosen. The details on how the field parameters and the trial dipole matrices were chosen is discussed in each case below. Each of the simulations below addresses particular aspects of the inversion procedure including the degree of local convexity about the true solution, the impact of various noise levels and the utility of increasing the amount of data. A clear pattern of the algorithm’s effectiveness will be evident from the collective results. 3.1

Idealized noise-free experiments

The goal of the following noise-free simulations is to test whether a gradient-based search algorithm can identify a Hamiltonian that reproduces the measurements and explore whether the myopic search yields a unique solution. The first set of simulations aim to identify a quantum system of dimension N = 4, where the diagonal elements of H0 are listed in ref. 34 and the true transition dipole matrix is randomly chosen as 0 1 0 1:4601 1:2416 0:9739 B C B 1:4601 0 1:7332 0:5390 C B C C: m¼B (11) B C B 1:2416 1:7332 C 0 1:4425 @ A 0:9739 0:5390 1:4425

0

During the dipole search, the diagonal elements were assumed to be zero, yielding 6 unknown Hamiltonian parameters. The number of P ki-f  P k1-4 measurements, k = 1,. . .,K, was taken to be three times the number of unknown parameters (here, K = 3  6 = 18). This entailed using 18 noise-free fields of the form in eqn (10), with L = 6 arising from the 6 unique energy transitions from H0. The 18 fields were characterized by distinct sets of phases and amplitudes randomly chosen from the domains [0,2p] and [0,1], respectively; the fields produced a range of P1-4 values. Upon performing multiple inversion tests starting with different trial dipoles m~, a consistent dipole should lie at the intersection of these distinct landscape excursions. With just 18 measurements we will assess the practical degree to which the inversion gives a good quality unique solution. For each inversion, a random trial transition dipole was used as a starting point on the inversion landscape. The inversion landscape was traversed by morphing over the transition dipole matrix elements in an effort to reach J B 0 following the D-MORPH procedure described in Section 2. The first set of 40 inversions was performed with trial dipole matrix elements randomly (uniformly) chosen to be within a 10–15% window about the true values in eqn (11). Knowledge of the dipole elements to at least this level of precision could be obtained, for example, via observed spectral intensities combined with electronic structure calculations. In each of

15618 | Phys. Chem. Chem. Phys., 2014, 16, 15615--15622

the 40 inversions, the same set of 18 fields ek(t) and P k1-4, k = 1,. . .,18 values were used. The D-MORPH search over dipole elements was driven towards a target J value of 108, and all 40 optimizations achieved this value. To characterize the quality of convergence to the true dipole, the distance metric Mm was determined, rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   (12) m  mÞ2 : Mm ¼ Tr ð~ The 40 inversion simulations produced a small range of Mm values, varying from 0.0011 to 0.0015. Excellent convergence quality was observed, and no simulations from within the starting window resulted in J - 0 producing an incorrectly identified dipole matrix within this level of precision. The second set of inversion simulations involved using trial dipoles lying within a 25–50% window about the true solution. Virtually identical results were obtained over 40 simulations using the same 18 data points and fields as in the previous set of inversions. All optimizations reached the target J B 108 and Mm ranged from 0.0011 to 0.0016. Simulations were also run where the initial dipole domain was expanded to comprise a 50–75% window about the true solution. Of these 40 inversion simulations performed with K = 18, a total of 13 yielded J values greater than 0.009 and none of these latter cases identified the true dipole. The remaining 27 runs produced J B 108 and showed excellent convergence to the true dipole. These overall results support the notion of local landscape convexity within a neighborhood about the true solution,27 especially notable here with just 18 data points and fields. While excellent results were observed using a system of dimension N = 4, it is important to consider the effect of system size on inversion quality. A system of dimension N = 8 was considered, with 28 dipole elements for identification. The known values of H0 are listed in ref. 35. A total of K = 3  28 = 84 simulated experimental P k1-8, k = 1,. . .,84 measurements were utilized for a single inversion. The trial m~ matrix elements were randomly chosen to be 10–15% away from the true values. The D-MORPH optimization procedure was able to achieve J B 108, and Fig. 2 plots the true values of m as squares and

Fig. 2 True dipole matrix element values (squares) and identified values (circles) for a system of dimension N = 8 where no noise was present. The dipole matrix contained 28 elements to be identified. The trial dipole elements were randomly chosen to be 10–15% away from the true values. A total of K = 3  28 = 84 simulated experimental measurements of P1-8 were used for this single inversion. Despite the large number of unknown parameters, excellent convergence is observed.

This journal is © the Owner Societies 2014

View Article Online

PCCP

Paper

the identified values as circles. Excellent dipole identification is observed, despite the large number of dipole matrix elements.

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

3.2

Experiments with noise

The next set of simulations include noise in both e(t) and Pi-f. ˜l and f~l were randomly chosen from the The values of a domains [0,1] and [0,2p], respectively. These values may be thought of as the best estimate of their true values in the laboratory, and they are used in the computational inversion with each trial dipole. To account for an experimentally reasonable degree of field uncertainty, random noise Zal and Zfl of 1–2% is incorporated in the fields, via al and fl, for generating the simulated experimental data: ˜l al = (1 + Zal)a

(13)

fl = (1 + Zfl)f~l.

(14)

and

This procedure is utilized for each of the fields k = 1,. . .,K. Noise is also added to each Pi-f = P1-4 simulated measurement through ˆ1-4 = (1 + ZP )P1-4 P 1-4

(15)

where ZP1-4 assumes a random value of 1.5–2% for the kth ˆ1-4 is forced to stay within the interval experiment, though P [0,1]. The noise model implemented in these simulations likely will not fully represent a realistic view of dealing with all aspects of noise encountered in an experimental setting; the goal here is to utilize the present model to explore the impact of noise on practical landscape convexity as well as the ability of the myopic search algorithm to identify the true Hamiltonian under these circumstances. Fig. 3A shows the final value for J versus the distance Mm of the identified transition dipole matrix from m in eqn (11) for 40 inversion simulations. Each inversion used the K = 18 simulated experiments from Section 3.1. The elements of the 40 trial dipole matrices were randomly chosen to be 10–15% away from the true system. The D-MORPH optimization procedure was employed until J(s + 1)  J(s) B 1010, with the target J B 108; here the parameter s labels the steps in the D-MORPH procedure. The final J values over the 40 inversions ranged from 2.0  104 to 9.2  104, and Mm was in the range 0.069–0.24. Fig. 3B shows the distribution of identified dipole matrix elements as circles with the true values displayed as squares, where the index (1–6) on the abscissa respectively corresponds to m(1,2), m(1,3), m(1,4), m(2,3), m(2,4), m(3,4). The mean values for each identified dipole element, denoted m(i,j), are m(1,2) = 1.46, m(1,3) = 1.23, m(1,4) = 0.97, m(2,3) = 1.71, m(2,4) = 0.53, m(3,4) = 1.43 with standard deviations of 0.040, 0.048, 0.026, 0.065, 0.040, 0.041, respectively. There is surprisingly good convergence given the noise incorporated into both the experimental fields and measurements. Importantly, as with the noise-free simulations, no systems approached J B 0 with a significantly erroneous dipole, indicating that the proposed local convexity27

This journal is © the Owner Societies 2014

Fig. 3 Results for 40 inversion simulations for a system with a Hamiltonian of dimension N = 4, where the target m is listed in eqn (11). The trial dipole matrices (i.e., the initial guess for solving eqn (5)) are randomly chosen to be 10–15% away from the true solution. In these simulations, noise was incorporated into the applied field amplitudes and phases (cf., eqn (13) and (14), respectively) as well as the measured P1-4 value (cf., eqn (15)). The noise was 1–2% for the amplitudes and phases and 1.5–2% for P1-4. In each simulation, K = 18 simulated experimental measurements were utilized for inversion. (A) Shows the inversion cost function J versus the distance between the identified and true dipole matrices, Mm (cf., eqn (12)). (B) Shows the distribution of identified dipole matrix elements (true values are shown as squares that can be readily identified with the particular elements in eqn (11)). While the presence of noise prevents complete convergence to either J B 0 or Mm B 0, the dipole elements are identified very well. The results suggest local landscape convexity about the true system.

holds even when noise is incorporated into experimental fields and measurements. Fig. 4 shows results from 40 inversions with the same degree of noise as above, but with the initial m~ trial matrix domain

Fig. 4 Results for 40 inversion simulations for a Hamiltonian of dimension N = 4, where the true matrix m is given in eqn (11). The trial dipole matrices (i.e., the initial guess for solving eqn (5)) are randomly chosen to be 50–75% away from the true solution. In these simulations, noise was incorporated into the applied field amplitudes and phases (cf., eqn (13) and (14), respectively) as well as the measured P1-4 (cf., eqn (15)) to simulate experimental conditions. The noise was 1–2% for the amplitudes and phases and 1.5–2% for P1-4. In each inversion, K = 18 simulated experimental measurements were utilized. (A) Shows the inversion cost function J versus the distance between the identified and true dipole matrices, Mm (cf., eqn (12)) and (B) shows the distribution of identified dipole matrix elements (true values shown as squares). There are two distinct regions visible in (A), and (C) focuses on the 24 out of 40 runs that yielded J o 0.0011. These runs show much better convergence to eqn (11), as shown in (D). The results suggest local landscape convexity about the true system, as simulations producing J - 0 all gave good quality dipole inversion. Importantly, those runs in (C) and (D) that show good convergence are comparable to the results in Fig. 3 where the trial dipoles matrices evidently were all within the locally convex domain.

Phys. Chem. Chem. Phys., 2014, 16, 15615--15622 | 15619

View Article Online

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

Paper

greatly expanded, from 10–15% to now 50–75% away from eqn (11). These simulations again incorporated the input amplitudes and phases from the K = 18 experiments used previously, with noise added to both the fields and P1-4 values to simulate experimental conditions. Fig. 4A shows the inversion cost function J versus Mm and Fig. 4B presents the distribution of identified matrix element values (following the dipole element index from Fig. 3). In Fig. 4A, two distinct regions can be identified: 24 of the 40 runs comprise a set of simulations where J was r0.0011, while the remaining 16 yielded larger values of J up to 0.027. Using all 40 runs, the mean extracted dipole element values are m(1,2) = 1.34, m(1,3) = 1.19, m(1,4) = 1.22, m(2,3) = 1.73, m(2,4) = 0.58, m(3,4) = 1.60. The associated standard deviations are 0.54, 0.65, 0.55, 0.72, 0.63, and 0.74, which are all an order of magnitude larger than those from Fig. 3B. Fig. 4C enlarges the results from the 24 aforementioned simulations achieving J r 0.0011 from Fig. 4A, and Fig. 4D shows the corresponding matrix element distribution. The mean values of the extracted dipole elements from Fig. 4D are m(1,2) = 1.46, m(1,3) = 1.26, m(1,4) = 0.97, m(2,3) = 1.74, m(2,4) = 0.54, m(3,4) = 1.45 with standard deviations shrinking to 0.035, 0.054, 0.032, 0.050, 0.040, and 0.053. The results from Fig. 4C and D suggest that good convergence can be obtained even when trial matrices are chosen from a wider domain, provided that those achieving a sufficiently small J value are utilized. In this regard, the results in Fig. 4C and D are much like those in Fig. 3A and B, respectively, as well as for the standard deviations of the identified dipole elements. Fig. 5 shows results from 40 inversions with trial dipole matrices taken from a domain of 50–75% away from eqn (11).

Fig. 5 Results for 40 inversion simulations for a Hamiltonian of dimension N = 4, where the true matrix m is given in eqn (11). The trial dipole matrices (i.e., the initial guess for solving eqn (5)) are randomly chosen to be 50–75% away from the true solution. In these tests, noise was incorporated into the applied field amplitudes and phases (cf., eqn (13) and (14), respectively) as well as the measured P1-4 (cf., eqn (15)) to simulate experimental conditions. Compared to Fig. 4, the noise was increased to 2–5% for amplitudes, phases, and P1-4 values. In each inversion, K = 18 experimental measurements were utilized. (A) Shows the inversion cost function J versus the distance between the identified and true dipole matrices, Mm (cf., eqn (12)) and (B) shows the distribution of identified dipole matrix elements (true values shown as squares). There are two distinct regions visible in (A), and (C) focuses on the 27 out of the 40 runs that yielded J o 0.005. These runs show much better convergence to the true dipole matrix, as shown in (D).

15620 | Phys. Chem. Chem. Phys., 2014, 16, 15615--15622

PCCP

The noise levels of Zal, Zfl, and ZP1-4 were increased to 2–5%. Similar to Fig. 4A, two regions can be seen in Fig. 5A, which plots J versus Mm: one region with J 4 0.005 and one with J o 0.005. A total of 27 out of the 40 simulations comprise this latter set. The distribution of identified matrix elements from all 40 simulations is shown in Fig. 5B. The mean extracted dipole elements over all 40 runs are m(1,2) = 1.23, m(1,3) = 1.12, m(1,4) = 1.06, m(2,3) = 1.82, m(2,4) = 0.72, m(3,4) = 1.58 with standard deviations of 0.56, 0.41, 0.35, 0.76, 0.42 and 0.46, respectively. Fig. 5C shows J versus Mm for the 27 out of 40 runs with J o 0.005, and Fig. 5D shows the corresponding distribution of identified matrix elements. The mean extracted dipoles from this subset are m(1,2) = 1.46, m(1,3) = 1.28, m(1,4) = 0.97, m(2,3) = 1.68, m(2,4) = 0.55, m(3,4) = 1.49 with standard deviations reduced to 0.10, 0.14, 0.58, 0.14, 0.097, and 0.14, respectively. The convergence improves with this subset of simulations but is inferior to that observed in Fig. 4D where less noise is present. To explore the impact of increasing the amount of data used for inversion, Fig. 6 shows results from 40 simulations where (a) trial dipole matrices were taken from a domain 50–75% away from eqn (11), (b) the incorporated noise was on the order of 2–5%, and (c) K = 30 experimental measurements were used compared to the 18 utilized in the previous sets of simulations. Conditions (a) and (b) are the same as the illustrations in Fig. 5. Fig. 6A displays J versus Mm, and again two regions are observed: one region with J 4 0.0021 and one with J o 0.0021. The latter set comprised 26 out of the 40 runs. Fig. 6B shows the distribution of identified matrix elements for all 40 runs; the mean extracted dipole values are m(1,2) = 1.42, m(1,3) = 1.32, m(1,4) = 0.99, m(2,3) = 2.00, m(2,4) = 0.51, m(3,4) = 1.57 with standard deviations are 0.44, 0.72, 0.25, 0.67, 0.17, and 0.54, respectively.

Fig. 6 Results for 40 inversion simulations for a Hamiltonian of dimension N = 4, where the true m is given in eqn (11). The trial dipole matrices (i.e., the initial guess for solving eqn (5)) are randomly chosen to be 50–75% away from the true solution. In these simulations, 2–5% noise was incorporated into the applied field amplitudes, phases, and Pi-f values (cf., eqn (13)–(15) respectively). In each inversion, K = 30 experimental measurements were utilized for inversion. (A) Shows the inversion cost function J versus the distance between the identified and true dipole matrices, Mm (cf., eqn (12)) and (B) shows the distribution of identified dipole matrix elements (true values shown as squares). There are two distinct regions visible in (A), and (C) focuses on the 26 out of the 40 runs that yielded J o 0.0021. These runs show much better convergence to the true dipole matrix in (D), as well as improvement over the results from Fig. 5 which utilizes less data.

This journal is © the Owner Societies 2014

View Article Online

PCCP

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

Fig. 6C zooms in on Fig. 6A to consider the 26 cases where J o 0.0021, and Fig. 6D displays the associated matrix element values. The mean extracted dipole values are m(1,2) = 1.47, m(1,3) = 1.23, m(1,4) = 0.96, m(2,3) = 1.76, m(2,4) = 0.54, m(3,4) = 1.42 and the distribution is much smaller compared to 6B, with standard deviations decreasing to 0.11, 0.044, 0.077, 0.11, 0.049, and 0.068, respectively. Overall the convergence to m (cf., eqn (11)) is improved using more experiments compared to Fig. 5, though importantly K = 30 still constitutes a very modest amount of data to consider for inversion.

4 Conclusions This work explored the Hamiltonian inversion landscape using simulated experiments, where the state-to-state transition probability at a single time was measured after the system interacted with a randomly chosen, resonant applied field. Two classes of inversion simulations were considered. First, noise-free inversions were performed where both the fields and the measured transition probabilities were assumed to be precisely known. The results of these simulations indicated that convergence to a landscape minimum corresponds to identifying the true Hamiltonian (i.e., a level set of multiple solutions was not observed), reflecting the proposed local convexity.27 The second class of simulations incorporated noise into both the field and measurements. The results show that good identification can be achieved with a relatively small amount of noisy data. It was observed that as the neighborhood is expanded from which the initial (trial) dipoles are chosen, some trials residing outside the evident domain of convexity are unable to converge to the true Hamiltonian. However, the subset of runs that did reach small J values all corresponded to good quality extracted dipole matrix elements. This feature is very important in practice as the convex domain will not be known a priori, and focusing on those runs of good convergence provides a means to deal with this situation. The convergence towards the true Hamiltonian is likely aided by the high sensitivity of the quantum dynamics data. It is interesting to note that quantum control landscapes also display convexity, but with non-unique control solutions existing due to an infinite null space about the optimal objective value.21 Such a multitude of solutions to choose from is ideal for the purpose of achieving optimal quantum control. The convexity on the inversion landscape, however, lends itself to a single solution (i.e., the true Hamiltonian) lying at the intersection of dynamic events reflected in the set of data. A practical algorithmic bottleneck to the performance of experimental inversion in sync with laboratory timescales is the computational component. In this regard, the results here suggest that ‘open-loop’ inversion can be feasible with a local search algorithm, including in the case where field and data noise is present. This work did not make any attempt to optimize the experiments,15,16 as there is substantial computational cost in implementing such a closed-loop inversion optimization procedure; however, future work may consider such a setup using a local Hamiltonian search algorithm. The systems considered in

This journal is © the Owner Societies 2014

Paper

this work were of low finite dimension. Simulations also should be performed which consider substantially larger systems to understand the relationship between the number of unknown Hamiltonian parameters and the number of experiments needed to assure convergence. The proposed algorithm ultimately needs to be refined in a laboratory setting where real noise is involved and a host of experimental issues can come into play.

Acknowledgements A.D. gratefully acknowledges support from the Program in Plasma Science and Technology at Princeton University. The authors also acknowledge support from the NSF (CHE-1058644), ARO-MURI (W911NF-11-1-2068) and ARO (W911NF-13-1-0237).

References 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

R. Rydberg, Z. Phys., 1932, 73, 376. O. Klein, Z. Phys., 1932, 76, 226. A. Rees, Proc. Phys. Soc., London, 1947, 59, 998. R. Gerber, R. Roth and M. Ratner, Mol. Phys., 1981, 44, 1335. R. Roth, M. Ratner and R. Gerber, Phys. Rev. Lett., 1984, 52, 1288. T.-S. Ho and H. Rabitz, J. Chem. Phys., 1989, 90, 1519. T.-S. Ho and H. Rabitz, J. Phys. Chem., 1993, 97, 13449. C. D. Franco, M. Paternostro and M. Kim, Phys. Rev. Lett., 2009, 102, 187203. D. Burgarth and K. Maruyama, New J. Phys., 2009, 11, 103019. D. Burgarth, K. Maruyama and F. Nori, Phys. Rev. A: At., Mol., Opt. Phys., 2009, 79, 020305. K. Maruyama, D. Burgarth, A. Ishizaki, T. Takui and K. B. Whaley, Quantum Inf. Comput., 2012, 12, 0763. S. Schirmer and D. Oi, Laser Phys., 2010, 20, 1203. A. Shabani, M. Mohseni, S. Lloyd, R. Kosut and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 012107. J. Zhang and M. Sarovar, 2014, arXiv:1401.5780. J. Geremia and H. Rabitz, J. Chem. Phys., 2003, 118, 5369. J. Geremia and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2004, 70, 023804. D. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, 1989. S. Forrest, Proceedings of the Fifth International Conference on Genetic Algorithms, Morgan Kaufmann, 1993. H. Rabitz, M. Hsieh and C. Rosenthal, Science, 2004, 303, 1998. T.-S. Ho and H. Rabitz, J. Photochem. Photobiol., A, 2006, 180, 226. V. Beltrani, J. Dominy, T.-S. Ho and H. Rabitz, J. Chem. Phys., 2011, 134, 194106. A. Donovan, V. Beltrani and H. Rabitz, Phys. Chem. Chem. Phys., 2011, 13, 7348. K. Moore and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2011, 84, 012109. K. Moore and H. Rabitz, J. Chem. Phys., 2012, 137, 134113. J. Roslund and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2009, 79, 053417.

Phys. Chem. Chem. Phys., 2014, 16, 15615--15622 | 15621

View Article Online

Paper

31 D. Meshulach and Y. Silberberg, Nature, 1998, 396, 239. 32 T. Brixner, N. Damrauer, P. Niklaus and G. Gerber, Nature, 2001, 414, 57. 33 J. Herek, W. Wohleben, R. Cogdell, D. Zeidler and M. Motzkus, Nature, 2002, 417, 533. 34 H0(1,1) = 1.5074, H0(2,2) = 9.9261, H0(3,3) = 12.9429, H0(4,4) = 14.0673. 35 H0(1,1) = 2.1486, H0(2,2) = 5.1589, H0(3,3) = 6.0121, H0(4,4) = 6.0313, H0(5,5) = 6.0542, H0(6,6) = 8.7463, H0(7,7) = 9.9811, H0(8,8) = 12.9747.

Published on 23 June 2014. Downloaded by Temple University on 23/10/2014 05:00:34.

26 Q. Sun, I. Pelczer, G. Riviello, R.-B. Wu and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 89, 033413. 27 Z. Leghtas, G. Turinici, H. Rabitz and P. Rouchon, IEEE Trans. Autom. Control, 2012, 57, 2679. 28 C. L. Bris, M. Mirrahimi, H. Rabitz and G. Turinici, ESAIM: COCV, 2007, 13, 378. 29 A. Rothman, T.-S. Ho and H. Rabitz, Phys. Rev. A: At., Mol., Opt. Phys., 2005, 72, 023416. 30 V. Beltrani, J. Dominy, T.-S. Ho and H. Rabitz, J. Chem. Phys., 2007, 126, 094105.

PCCP

15622 | Phys. Chem. Chem. Phys., 2014, 16, 15615--15622

This journal is © the Owner Societies 2014

Exploring the Hamiltonian inversion landscape.

The identification of quantum system Hamiltonians through the use of experimental data remains an important research goal. Seeking a Hamiltonian that ...
1MB Sizes 2 Downloads 3 Views