Assessment of optimal control mechanism complexity by experimental landscape Hessian analysis: fragmentation of CH2BrI

Optimally shaped femtosecond laser pulses can often be effectively identified in adaptive feedback quantum control experiments, but elucidating the underlying control mechanism can be a difficult task requiring significant additional analysis. We introduce landscape Hessian analysis (LHA) as a practical experimental tool to aid in elucidating control mechanism insights. This technique is applied to the dissociative ionization of CH2BrI using shaped fs laser pulses for optimization of the absolute yields of ionic fragments as well as their ratios for the competing processes of breaking the C–Br and C–I bonds. The experimental results suggest that these nominally complex problems can be reduced to a low-dimensional control space with insights into the control mechanisms. While the optimal yield for some fragments is dominated by a non-resonant intensity-driven process, the optimal generation of other fragments maa difficult task requiring significant additionaly be explained by a non-resonant process coupled to few level resonant dynamics. Theoretical analysis and modeling is consistent with the experimental observations.


Introduction
The optimal control of quantum system dynamics using tailored laser fields is experiencing growing success over a diverse array of experiments [1][2][3][4][5]. In many of these studies, adaptive feedback control (AFC) [6] is used to find an optimal solution with a search algorithm (e.g., genetic algorithm (GA)), which is proving to be of wide utility for manipulating complex systems even without prior detailed knowledge of the system Hamiltonian. While search algorithms are effective at finding an optimal solution, extracting insights about the control mechanism often calls for additional effort. The challenge in understanding the mechanism is tied to the common situation of operating with a large number of control variables and the nonlinearity of the dynamics with respect to the field. In some cases, prior knowledge of the system has permitted expressing the control in terms of a small number of parameters, and principal component analysis [7,8] has found use in seeking a basis set of control variables with reduced dimensionality [9][10][11]. An overall goal is to gain knowledge about the effective complexity of the quantum system expressed in terms of the nature of the participating states under control, which is often not easily deduced from the nominal experimental data. Quantum control landscape analysis, under reasonable assumptions for closed quantum systems, shows that the effective complexity is determined by the number and nature of the system states that participate in the laser-matter interaction [12,13]. This result has important implications in assessing the intricate features of quantum systems under control, which is the subject of this study.
The features of the quantum control landscape (i.e., the yield or objective J(x) as a function of the vector of control field parameters x) contain rich information about system under manipulation. In particular, the topology at the maxima of a quantum control landscape has important physical significance. Under optimization of the objective, the best laser field E opt (t), corresponding to a point at the top of the landscape, will exploit the favorable dynamical features of the system. The complexity of the controlled dynamics can be assessed by analyzing the sensitivity of the yield to variations in the control at the top of the landscape. At such a landscape critical point x 0 , the gradient satisfies form the Hessian matrix, which contains the sensitivity information to leading order in the variations of the control in the neighborhood of x 0 . Experimental determination of the Hessian can be carried out by statistically sampling the vicinity of the landscape optimum through small random perturbations δx to the control parameters (e.g., the spectral phases) and observing the resultant changes δ + J x x ( ) 0 . Diagonalizing the Hessian matrix then reveals the effective freedom of movement near the landscape extremum reflected in the set of eigenvectors with negative eigenvalues. While the number and nature of the negative eigenvalues of the Hessian matrix identify the dynamical complexity of the quantum system, the eigenvectors establish how the control parameters work cooperatively to achieve optimal control. Quantum control landscape analysis shows that, upon satisfaction of particular physical assumptions [13,14], for a closed quantum system with N accessible energy levels, the Hessian rank (e.g., the number of nonzero Hessian eigenvalues) will lie between 1 and 2N-2 for a state-to-state transition process. This theoretical prediction was validated by many simulations [15] and in an experimental study which measured the Hessian at the fitness maximum for electronic transitions of atomic Rubidium [16]. The latter study used a shaped fs laser pulse with sufficient spectral coverage of the two resonant pathways of Rb. The results showed 2N-2 = 6 discrete negative eigenvalues, revealing an effective N = 4 level system for atomic Rb, which is consistent with the lowest-order dynamics expected at the laser powers employed. Another experimental study of the control landscape for proton NMR control in HDO showed a Hessian rank of 2 at the landscape top and bottom critical points, which also agrees with the theoretical rank prediction for a two-level system [17]. At the other extreme are systems involving only non-resonant transitions, such as two photon absorption (TPA), second harmonic generation (SHG), and multiphoton excitation. The thoeretical and experimental coherent control studies of multiphoton non-resonant transitions in the weak field regime were pioneered by Meshulach and Silberberg [18,19]. The frequency-domain analysis of TPA was later extended to the intermediate field regime, where a broad family of solutions were identified that could enhance the absorption [20][21][22]. While makeing the transition from the weak to strong field regime, the static resonance and perturbative strategies in the frequncy domain is no longer suitable, but the dynamic insights may still be captured with a time-domain picture addressing the ground state depletion and dynamic Stark shifting [23,24]. In the weak field limit, the experimental and theoretical Hessian analysis of TPA were carried out at the landscape top, which produceed a full rank Hessian matrix and a continuous Hessian spectrum resolved to the pulse shaper spectral resolution [25,26]. In the present study, for general −  photon non-resonant absorption, we derive an analytical expression for the Hessian. Thus, we have available the Hessian control landscape signatures for the limiting cases for either resonant few level systems or non-resonant multiphoton driven systems. We will show that knowledge of these limiting cases provides useful references to aid in interpreting the observed Hessian spectrum of complex systems including for controlled molecular fragmentation and ionization, where both resonant and non-resonant transitions may be involved.
The molecule chosen for analysis in this study is CH 2 BrI, whose bond selective photochemistry has been investigated both experimentally [27,28] and theoretically [29][30][31]. The control of dissociative ionization for the halomethane family of molecules has been treated in several works [32][33][34][35][36][37][38]. An early use of AFC demonstrated selective bond breaking and enhancement of the product ratio of by a factor of ∼1.8 in the fragmentation of CH 2 BrCl with shaped laser pulses [32]. Another study showed that the + + CH I CH Br 2 2 ratio from CH 2 BrI can be controlled with a few variables through a linear combination of a suitable basis [10]. Pump-probe experiments in CH 2 BrI identified oscillations and anti-correlations in the temporal dependence of the ion yields of CH 2 Br + and the parent ion CH 2 BrI + which were explained in terms of a simple control mechanism involving a one-photon ionic resonance [33]. These findings were corroborated by quantum chemical calculations that also revealed the role of the Br-C-I bending mode dynamics on the low lying electronic states of the ion [34,35]. In addition to the ionic resonance acting as one of the control mechanisms for dissociative ionization, resonance effects in the high-lying neutrual exicted states were also found to play an important role in the resonance-enhanced strong field ionization in CH 2 BrI [39,40]. Our recent study investigated the optimal control of a family of halomethane molecules [38] finding that the ratios X + /CH 2 X + (where X represents a halogen atom: F, Cl, Br, I) can be greatly enhanced with shaped laser pulses; this work also showed systematic behavior with respect to the electronic characteristics of the chemical reagents and the nature of the 'photonic reagents' (shaped pulses). Another set of experiments on controlled dissociative ionization of CH 2 XY determined the control landscape expressed in terms of 2nd, 3rd and 4th order polynomial phase variables. The flexibility in controlling the ratio X + /CH 2 X + over a large dynamical range was found to be correlated to distinct features in the landscapes of X + and CH 2 X + . In the present paper, we further explore the control landscapes for these fragment ion species (i.e., X + and CH 2 X + from the CH 2 BrI molecule), but now in the high dimensional pixelated phase basis and focus on extracting mechanistic insights from simple experimental landscape Hessian analyses.
The remainder of the paper is organized as follows. Section 2 describes the procedure to obtain noise-robust Hessian matrices using a statistical sampling method, which is employed in both the experiments and simulations given here. Section 3 presents experimental landscape Hessian analysis (LHA) results of CH 2 BrI for the optimal control of fragment yields and species ratios. In section 4 we evaluate the Hessian matrices for the idealized models: (1) a resonant few level system, (2) non-resonant  -photon absorption, and (3) the mixed situation in which nonresonant multiphoton transitions take place coupled with resonant transitions. These model simulations are performed as a foundation to understand the qualitative features of Hessian spectra in complex situations, such as for the current experiments with CH 2 BrI. The overall goal of the paper is to introduce LHA and present its implementation. The focus will be on the extracted Hessian information, and selected detailed features of the optimal control fields will be discussed in keeping with this goal. The collective experimental and simulation results provide a basis to suggest that rather simple models, possibly even of a quantitative nature may be generated to describe the nominally complex processes occurring during laser driven dissociative ionization of polyatomic molecules. Section 5 provides concluding remarks.

Introducing LHA: methods and procedures
The fitness landscape δ in the neighborhood of an optimum x 0 is approximated by the local Hessian H, which is the matrix of second derivatives: where δx is a vector of small control parameter variations, and x 0 is the vector of parameters representing the optimal control field at a landscape maximum, which may be obtained from an AFC experiment. The gradient is neglected, as we are at or very close to the landscape maximum in the laboratory. Due to experimental noise, a determination of the Hessian matrix at x 0 through finite differences was found to be unreliable. Instead, we employ a stable statistical method [16,41] based on taking second moments of the experimental data. The landscape around x 0 is first sampled using a Gaussian distribution σ N (0, ) as a separable product in each of the control variables = … x i , 1, 2, i with σ being the standard deviation of the sampling distribution. The Hessian matrix elements H ij at the optimum can be determined from the sampled points by a yield-weighted correlation matrix of the form where 〈 〉 · denotes an average over the sampled data and ε ij consists of residual 4th or higher order derivative terms, which are assumed to be insignificant. In the following, we will use an example to illustrate the procedure to extract the Hessian and its eigenvalues and eigenvectors particularly considering experimental noise issues.
From equation (2.1), it is evident that a Hessian which is strictly diagonal indicates that the frequency components (pixels in the pulse shape) act independently of each other in their impact on the landscape fitness value. In contrast, the presence of off diagonal elements of the Hessian indicates specific cooperative relations between the frequency components for their impact on the fitness. In practice, even the cases that have a strong diagonal Hessian character can also have some modest off-diagonal features, while we will also see clear cases with strong off-diagonal Hessian character. Figure 1 presents the steps in the statistical determination of the Hessian illustrated for an experimental TPA signal, where x 0 is a length 80 vector of spectral phase variables equally spread over the range from 755 to 819 nm. A few percent of the amplified shaped fs laser pulse is focused into a GaP photodetector (Thorlabs Inc.), which gives rise to a TPA signal by absorbing two ∼790 nm photons. The fs laser system and pulse shaper are described in section 3. The transform-limited pulse corresponding to a reference flat phase (i.e., spectral phase ≡ x 0 = 0) is identified, which maximizes the TPA signal. The sampling is then performed by small simultaneous variations of all spectral phase components drawn from a Gaussian product distribution σ N (0, ) around zero. This procedure generates the L = 40 000 data points shown in figure 1(a) where the maximum TPA signal is normalized to 1. Since the true Hessian is retrieved when σ → 0, the choice of σ should be as small as possible. On the other hand, due to experimental noise, σ should be large enough to ensure that the variation of the objective yield δ is significantly larger than the noise level. In this example, we choose σ = 0.3 to satisfy the above two criteria, and the resultant variation of the objective yield lies in a window to ∼20% below the maximum, as shown in figure 1(a), which is significantly larger than the statistical noise of ∼2% in the experimental TPA signal. A subset of M data points are randomly picked from the full set L ( < M L) to determine the Hessian matrix elements using equation (2.2), followed by diagonalization of the matrix to obtain the eigenvalues and eigenvectors. This statistical bootstrap resampling procedure [42] is performed at least 20 times using M new random data points in order to determine the mean eigenvalues and their standard deviations at a particular M. This process is repeated by selecting different sampling sizes M (here, from 3000 to 30 000 with small steps). The statistical errors for the eigenvalues can be estimated from the data averaging in the bootstrap procedure with the error decreasing with increasing M. The statistical errors introduced from finite data sampling in this procedure should be distinguished from the experimental noise, and computer simulations indicate that the extracted eigenvalues are stable with respect to reasonable experimental noise (not shown here). An estimate of the true Hessian eigenvalues at the landscape top may be obtained by extrapolating the results for → ∞ M . Since the statistical method is based on Monte Carlo sampling, its rate of convergence is expected to scale as ∼ M 1 [43]. Additional error could be introduced in the extrapolation, which depends on the selected window and method used for the extrapolation. Figure 1 . Although, bootstrapping is effective for extracting the eigenvalues, it does not readily apply to obtaining the eigenvectors, especially in cases where neighboring eigenvalues are nearly degenerate, corresponding to non-unique eigenvectors. Therefore, later in this paper, we simply use all the available sampling points L to calculate a best estimate of the eigenvectors 〉 v {| } n ; the complete set of eigenvectors for the TPA example will be shown later in figure 6(a). A final 'reconstructed' Hessian matrix is formed by where λ { } n are the extrapolated eigenvalues obtained by the bootstrap procedure. The latter reconstruction is convenient for qualitative graphical purposes, as the resultant Hessian matrices look similar, but less noisy, than that calculated directly from all the sampling points L. Figure 1(d) shows the reconstructed Hessian matrix, which is dominated by diagonal elements with a Gaussian-like distribution around the center wavelength. The wavelength scaling is mapped to 80 pixels as explained in section 3.
This statistical method for Hessian determination at the landscape maxima is employed for all the experiments in section 3, as well as for the simulations in section 4 to assure consistency.

Experimental implementation of LHA: an example of controlled framgmentation of CH 2 BrI
A detailed description of the experimental setup is given in [38]. Briefly, our experiments make use of a Ti:Sapphire amplified laser system (KMLab, Dragon, 3 kHz) producing 1mJ, sub 30 fs pulses centered ∼790 nm. The laser pulses were introduced into a pulse shaper containing a programmable dual-mask liquid crystal spatial light modulator (SLM) with 640 pixels and a resolution of ∼0.20 nm/pixel (CRI, SLM-640) operated with only phase modulation (i.e., at 100% amplitude transmission). Half of the pixels (320 pixels) are used and every group of 4 pixels are bundled to form a single parameter resulting in 80 independent control parameters to have complete coverage of the effective laser bandwidth. Pixel 1 and 80 correspond to ∼755 nm and ∼819 nm, respectively, and the maximum spectral intensity is located near pixel 40 (∼787 nm). This scale is shown in figure 1(d) and the subsequent Hessian plots in the paper. The laser pulse is focused into a vacuum chamber with a spherical lens ( f = 20 cm), resulting in a peak intensity of ∼ × − 1 10 W cm 14 2 . Samples of CH 2 BrI were introduced into the vacuum chamber through an effusive leak valve to maintain a constant pressure of 1 × 10 −6 Torr throughout the experiment. A time-of-flight spectrometer (TOF-MS) is connected to the interaction region, which measures the mass-to-charge ratio of ions produced in the laser focus. The resulting TOF spectrum is recorded with a digital oscilloscope (LeCroy 104MXi) and averaged over at least 1000 laser shots.

LHA of individual fragment species from CH 2 BrI
The first step towards experimental implementation of the LHA technique is to identify the top of the landscape for a particular product species by finding an optimal pulse. All of the optimizations were performed using a GA with the control field parameters as the 80 spectral phases in the pixel basis. We focus on the four major ionic fragments from dissociative ionization of CH 2 BrI: Br + (m/q = 79, 81), CH 2 Br + (m/q = 93, 95), I + (m/q = 127), and CH 2 I + (m/ q = 141). The absolute yield J(x) for each ionic fragment is separately optimized to locate its field x 0 . The resulting TOF spectra for each individual optimization are shown in figures 2(a)-(d), and the vertical scales for all the spectra are kept identical. The optimal pulses for Br + and I + are similar and significantly longer in duration than the transform-limited pulse, thereby resulting in TOF spectra with the optimal pulses having enhanced yields of X + and reduced yields of CH 2 X + . The optimal pulses for CH 2 Br + and CH 2 I + are much shorter, but can still produce a respective yield better than the transform-limited pulse. While the TOF spectra optimizing CH 2 Br + and CH 2 I + look similar to each other, distinct fragmentation patterns are evident when compared to those optimizing Br + or I + .
Following the method in section 2, statistical sampling is then performed through small random phase variations with a Gaussian distribution σ = N (0, 0.3) around the optimal phase solution for each fragment. The data δ is taken from the integrated signal of the selected fragmentation ion peak. The optimal yield J x ( ) 0 is re-evaluated after each measurement of δ , and the fitness used to calculate the Hessian is taken as the normalized ratio δ + J x x J x ( ) ( ) 0 0 in order to reduce the effects of slow long term signal variation. In each case, the total random sampling size is ∼ L 20 000. Using the methods described in section 2, we first examine the reconstructed Hessian matrices in figures 3(a)-(d). All of the cases in figure 3 show matrices dominated by diagonal elements with negative values as well as some significant off-diagonal features. The Hessian matrices for Br + and I + (figures 3(a) and (b)) show more prominent off-diagonal features and couplings than the Hessian matrices corresponding to CH 2 Br + and CH 2 I + (figures 3(c) and (d)). In addition, the off-diagonal features for Br + and I + are discretely located at specific wavelengths (e.g., ∼780 nm, ∼793 nm for Br + as well as other selected values), indicative of cooperative effects [16] also reflected in the dominant Hessian eigenvectors (see figure 5).The assignment of the frequencies and couplings for maximizing any of the species requires additional analysis drawing on the electronic states and potential energy surfaces from quantum chemical calculations, which is beyond the scope of this paper. However, the LHA technique directly provides valuable dynamical information for building models drawing in the relevant physics. The Hessian of CH 2 Br + shows significant dense off-diagonal coupling from ∼760 to ∼785 nm, while CH 2 I + only exhibits modest off-diagonal features, and the continuous diagonal distributions are consistent with non-resonant intensity-driven processes [25,26]; also see the appendix. Interestingly, a blue shift of the highest amplitude diagonal element with respect to the laser center wavelength was observed for CH 2 Br + and CH 2 I + , which is distinct from the Hessian of TPA in figure 1(d) and other intensity-driven processes (see section 4.2) whose maximum diagonal element should coincide with the laser central wavelength having the peak intensity. The shift toward favoring energy blue photons evidently changes the statistical distribution of the ionic state populations and allows for access to more dissociative states giving rise to the respective fragment yields. In contrast, the optimal yields of Br + and I + show no pronounced blue shift, with their appearance energies Figure 2. TOF spectra from CH 2 BrI obtained by applying the optimal pulses that respectively maximize (a) Br + , (b) I + , (c) CH 2 Br + , (d) CH 2 I + , identified by arrows with 'opt'. For clarity, only the portion of the spectra that cover these four fragment ions is shown. The distinct fragmentation patterns and absolute yields for each fragment indicate that the landscape optima differ from each other, although there is similarity within the pairs (a), (b) and (c), (d).
being much higher than those for CH 2 Br + and CH 2 I + . As a result, achieving the optimal control of Br + and I + fragments is likely governed by a different mechanism involving indirect ionization without the need for high energy blue photons.
The key feature of LHA derives from the assessment of the Hessian spectra (i.e., eigenvalues and the associated eigenvectors). In each case, the eigenvalues are displayed in sequential order by their magnitudes. Interpreting such spectra nominally would be a complex task, but the experimental data here contain distinct features whose origin can be understood based on simple laser driven physical dynamics with supporting evidence from the model analysis in section 4. Figure 4 shows the Hessian spectra for the various fragments using the procedure described in section 2. With each case being at its respective landscape top, all of the Hessian spectra in figure 4 contain negative eigenvalues; in general, such spectra could be negative semi-definite. The most noticeable feature is that, the Hessian spectra of Br + and I + are quite distinct from those of CH 2 Br + and CH 2 I + , although each spectrum can be divided into a region of rapid (R) and slow (S) variation. The Hessian spectra for Br + and I + have a clearly pronounced R region indicating the presence of a few outstanding (discrete) eigenvalues along with a continuous S domain. In contrast, the spectrum for CH 2 Br + is nearly continuous with an evident change from its rapid variation R to a much slower trend S occurring at eigenvalue ∼25. The case of CH 2 I + is similar to CH 2 Br + , but CH 2 I + has two distinctive discrete eigenvalues in the R region, which are less prominent than those for Br + and I + , followed by a slowly varying The spectral errors are estimated to be similar to those shown in figure 1(c). All of the spectra can be divided into a region of rapid (R) eigenvalue variation (even to the degree with some outstanding discrete lines in (a), (b) and (d)) and into a region of slow (S) variation. The break between the R and S regions, shown as straight blue lines, lies at eigenvalue number ∼5, ∼6, ∼25 and ∼4 in figures (a)-(d), respectively. These features of the LHA spectra may be directly interpreted in terms of the control mechanisms, especially involving strong non-resonant processes, along with a few resonantly coupled states in cases (a), (b) and possibly (d).
continuum S domain. The magnitude of a particular eigenvalue weights the significance of the coordinated control variations identified by the associated eigenvector.
The paper is focused on presenting a physical interpretation of LHA, especially including distinctions with purely intensity driven processes, as with TPA in figure 1. The distinctions concern situations where: (i) outstanding discrete eigenvalues are present along with a continuous portion of the Hessian spectrum and (ii) evident off-diagonal coupling appears in the Hessian. In the case of TPA in section 2, we showed that the Hessian is nearly diagonal with only weak off-diagonal coupling. In addition, the spectrum for such a purely non-resonant process is Gaussian-like, which may be analytically established for the general case of  -photon non-resonant excitation (see figure 10 and appendix A). These results provide the foundation to attribute the continuous portion of the spectrum to the effect of a dominant nonresonant intensity driven process present in all of the fragments resulting from maximization of the dissociative ionization yields for CH 2 BrI in this study. The label 'dominant' refers to the impact of a random control variation collectively having a large overall response from the many small non-resonant eigenvalues compared to that from the few strong discrete eigenvalues. For the cases of CH 2 Br + and CH 2 I + , this intensity driven process is especially dominant. However, their Hessian spectra are distorted from a Gaussian-like shape (see figures 1(c) and 10) as well with CH 2 I + exhibiting two distinct eigenvalues, collectively indicating that the optimal control of these two products also involves dynamical mechanisms other than a simple intensity driven process. The relative importance of the non-resonant intensity driven contributions can be read directly from figure 4 as , which follows the ascending order of their appearance energies [44]. The discrete outstanding eigenvalues evident in the Hessian spectra of Br + , I + and to a lesser degree in CH 2 I + , suggest in these cases that at least a few discrete resonantly coupled states are involved (see section 4 and [13,14]) along with nonresonant dynamics leading to the landscape maximum. The TPA Hessian has strong diagonal dominance, which is also present in the Hessians of the CH 2 BrI fragments in figure 3. The detailed features of the CH 2 BrI Hessians is only partially explained by intensity effects, as will become apparent from an analysis of the respective Hessian eigenvectors below. The eigenvectors will also reflect the evident off-diagonal couplings in all of the CH 2 BrI Hessians, especially for the cases with Hessian spectra containing discrete components.
The observed distinctions in the Hessian spectra for CH 2 Br + and CH 2 I + over that of Br + , I + indicates that different mechanisms and fragmentation pathways are involved in achieving optimum control over their respective yields. Examination of the SHG-FROG measurements revealed longer optimal pulses for Br + or I + than those for CH 2 Br + and CH 2 I + (not shown here). Our previous study [45] on the controlled fragmentation of halomethane molecules suggested that the products of CH 2 X + and X + have distinct landscape structures and that their relative yields can be readily manipulated [38] in a low dimensional space of coordinated field variations over the spectral bandwidth. Observations of the TOF trace with a chirp scan in the latter study also showed that CH 2 X + yields favor a short and intense pulse (i.e., near zero chirp) while X + yields are maximized at a larger chirp. The latter increase in the pulse duration reduces the CH 2 X + yield significantly due to the lower peak intensity and opening up of the additional dissociation channel into the products CH 2 + and X + . The X + yield may be enhanced via the latter channel, which favors a longer pulse over the ever present intensity effect. A one-photon ionic resonance in the dissociative ionization of CH 2 BrI was identified [33] in a pump-probe experiment, which led to modulation of the CH 2 Br + yield. Our Hessian spectra do not show a clean signature of a few level resonance entering CH 2 Br + , but instead show up in Br + and I + . These results may appear to be contradictory, but the region of the landscape for CH 2 Br + studied with pump-probe is not at its global maximum where the present study is focused. A resonance playing a role in the control mechanism can be prevalent, as long as two electronic excited states are coupled with ∼800 nm radiation in the present circumstance. Thus, in the Hessian spectrum, the lines standing out of the continuum are expected to reflect the signature of a resonant few level vibronic system coupled in some fashion to the non-resonant dynamics. Further understanding of the physical nature and couplings between these two processes is examined in the modeling studies in section 4. Importantly, even the simple models considered there are already capable of showing the general origin of the Hessian spectral signatures we see for the fragments arising from CH 2 BrI. Additional future refinement and analysis could build on this work seeking more specific physical models quantitatively tuned to capture the LHA features of each fragment.
We also examined the eigenvectors of the Hessian matrices for the four fragments. In each case, the first eigenvector corresponding to the most negative eigenvalue (blue curve), along with the optimal spectral phase (red curve) is plotted in figures 5(a)-(d). We found that for (a) Br + and (b) I + the first eigenvector contains features covering a wide range of spectral components within the laser bandwidth which reflect a complex cooperative interaction amongst the spectral phase components, while for (c) CH 2 Br + and (d) CH 2 I + , the first eigenvector is more focused within a narrower window of the blue shifted spectral region. These observations in the eigenvectors are consistent with the structures seen in the respective Hessian matrices, where cases (a) and (b) show strong off-diagonal couplings and cases (c) and (d) show significant blue shifted features. In addition, we also compare the first eigenvector with the optimal spectral phase and calculate the correlation coefficients between these two vectors in each case. The correlation coefficient is defined as , where x and y are the elements of the two vectors (i.e., respective eigenvector and spectral phase in pixels), σ is the associated standard deviation and n is the number of points used to calculate the correlation coefficients. The pixel index that corresponds to the spectral wavelength is shown as the upper x-axis lables in figure 5. We find that there is strong correlations between the first eigenvector and the optimal spectral phase for case (a) = (|corr| 0.71), where the first eigenvalue is outstanding in magnitude. Less correlation is presented in case (b) = (|corr| 0.31) and (d) = (|corr| 0.30), while almost no correlation exists in case (c) = (|corr| 0.01).
As pointed out earlier and in the caption to figure 4, all of the Hessian spectra show a region of initial rapid variation R followed by a break to a domain of slow variations. This behavior is directly linked to the degree of off-diagonal coupling in the Hessian matrices shown in figure 3. As a comparative illustration, figure 6 shows the complete set of eigenvectors for the case of (a) TPA as a reference and (b) the fragment CH 2 Br + . A clear, spectrally centered 'V' shape structure is observed for the pure non-resonant case of TPA. This overall feature can be interpreted as the nth eigenvector having spectral character focused at those wavelengths intersecting the V. As the eigenvalues in figure 1(c) diminish in magnitude with increasing index n, it is evident from the eigenvector V structure that the physics naturally reflects utilization of the laser intensity profile centered at λ ∼ 787 nm. A similar, but blue shifted, V structure is observed for CH 2 Br + within the first ∼25 eigenvectors, which coincides with the break between the rapid R and slow S behavior in figure 4(c  Then, we have that where 1 s is a diagonal 80 × 80 matrix with 1 on the diagonal for the spectral domain of 755-764 nm and 785-819 nm, while the remaining diagonal elements are zero over 765-784 nm, and λ¯is the mean eigenvalue in the slow domain S. As a result, we can identify CH 2 Br + as having special intensity driven dynamics associated with its R domain centered at ∼775 nm. In addition there is a second intensity driven piece linked with the slowly varying spectral domain S. The full eigenvector matrices for the other three fragments (Br + , I + and CH 2 I + ) arising from CH 2 BrI are not shown, but they all contain discernible features, especially with the first few distinguishable eigenvectors prescribing a strong signature of the Hessian off-diagonal coupling structure in figures 3(a) and (b). Thus, although all of the spectra in figure 4 have R and S domains, the case of CH 2 Br + stands out as both of these domains appear to be associated with physically distinct intensity driven processes. The combined discrete and continuous spectra in the respective R and S domain of Br + , I + and CH 2 I + may be interpreted as signatures of laser driven dynamics of a few resonantly coupled levels tied to a strong non-resonant process.
The Hessian analysis exibits special characteristics in the optimal control of individual fragment species of (a) Br + , (b) I + , (c) CH 2 Br + , and (d) CH 2 I + from dissociative ionization of CH 2 BrI. The combination of results and analyses of the Hessian matrix, Hessian spectrum and Hessian eigenvectors, show degrees of both common and distinct features that highlight the underlying control mechanisms for the four different cases. The common features lie in the dominant diagonal elements in the Hessian matrices, the existance of slowly varying Hessian spectrum component (S), and a 'V' shape structure in the Hessian eigenvectors, all of which point to an intensity driven non-resonant process that prevails in all the cases. It should be noted that, the diagonal features in the Hessian matrix could also arise from resonant transitions at frequencies corresponding as well to non-resonant transitions, as often encountered in complex molecules. The disctinctive features are reflected in the off-diagonal elements in the Hessian matrices, the outstanding eigenvalues forming a rapid (R) varying region in some of the Hessian spectra as well as the correlations between the optimal pulse and the first eigenvector. Cases (a) and (b) show more of these additional distinctive features than cases (c) and (d), indicating that the effect of resonance processes are also involved and mixed with the non-resonance process is more apparent in the former two cases. Importantly, we see a transition of the leading nonresonance effect in the control mechanism to the presence of a resonance contribution with the increasing appearance energy in the control target (i.e., < < < + + + + CH Br CH I I Br ) 2 2 . In the next section, we will consider two additional cases, in which the species ratios rather than the individual yields are chosen as the control objectives.

LHA of species ratios from CH 2 BrI
Additional LHA studies were performed for objectives considering optimal fragment ratios instead of the absolute yield of a single fragment ion. The fragmentation products from competing processes are frequently treated for ratio optimization (e.g., selective bond breaking of C-Br and C-I in the present case), and in these cases the optimal fitness is often less dependent on the laser intensity and the resulting optimal fields are usually more complex [32,38]. Here, we consider two objectives that optimize the fragment ratios of either CH 2 Br + / CH 2 I + or Br + /I + , which respectively correspond to breaking the weak bond versus the strong bond and vice versa in the second case. Importantly, such ratio cost functions form legitimate landscapes, but they are the ratio of observables rather than a single observable. Thus, the Hessian of a ratio can be influenced by the local Hessian of each fragment giving rise to a more complex scenario to assess. Considering the new fitness J to be a ratio of two observed fragments A and B (i.e., = J x ( ) . Note that the optimal control x 0 for the ratio may differ from either control that optimizes the individual fragments of A or B. Taking the second derivative of is the normalized ratio at the maximum following the procedure discussed earlier. The latter result shows that, the Hessian of the ratio (J = A/B) at its landscape top is proportional to the difference of the Hessians for A and B. In particular, 0 , with the proportionality factor being B x 1 ( ) 0 . It is important to note that the Hessians H x ( ) A 0 and H x ( ) B 0 here are not evaluated at the landscape tops for the individual species A and B, but at the optimal field x 0 of the ratio. However, we will show in one of the cases below that the landscape maximum and local structure of A, B and A/B appear close to each other. In this particular circumstance, the Hessian spectrum of A/B can be approximated by the difference of the Hessian spectra for A and B evaluated at there individual maxima, which approximately coincide.
The AFC experiments for the ratio objectives identified optimal pulse solutions that gave rise to the ratios of 2.53 ± 0.02 for CH 2 Br + /CH 2 I + and 0.35 ± 0.01 for Br + /I + , which outperform the transform limited pulse by 5% and 32%, respectively. While the optimal pulse for the CH 2 Br + /CH 2 I + ratio is near transform limited, that for Br + /I + results in much lower laser peak intensities and smaller absolute yields for both fragments. The LHA for these two cases produce Hessian spectra which are respectively shown in figures 7(a) and (b) as red dots with drop lines. In addition, we also show the reconstructed Hessian matrices for the ratios as insets. The Hessian spectra for (a) CH 2 Br + /CH 2 I + and (b) Br + /I + exhibit general features similar to the spectra of their respective numerators (i.e., the Hessian spectrum of CH 2 Br + shown in figure 4(c) and for Br + shown in figure 4(a)), except for the tail portion of the spectra for the ratios being much closer to zero. The Hessian matrices for the ratios also show similar features when compared to the Hessian matrices of their respective numerator in figure 3 for each case. While the Hessian matrix for CH 2 Br + /CH 2 I + in case (a) shows a slight blue shifted feature along with clear off-diagonal couplings in the vicinity of the diagonal elements from ∼770 to ∼795 nm , the Hessian for Br + /I + in case (b) shows a strong feature at ∼780 nm along with significant coupling facets through the off-diagonal elements. Based on the discussion in the paragraph above, in each case, we took the difference of the Hessian spectra for the numerator and denominator and compared that with the Hessian spectrum of the ratio. We found very good agreement in the case of CH 2 Br + /CH 2 I + but not for Br + /I + . The difference Hessian spectrum of CH 2 Br + and CH 2 I + is plotted in figure 7(a) as blue stars (*), scaled by a constant in keeping with the anticipated proportionality. In support of this finding, the optimal fields identified for CH 2 Br + , CH 2 I + and CH 2 Br + /CH 2 I + had analogous features, resulting in similar TOF spectra in signal strength and pattern. Furthermore, the three Hessian matrices H A , H B ,  figure 4(c) is also plotted in (a) as blue stars(*), which is scaled and matches well with the red dots. This striking similarity suggests that a common control mechanism is operating for maximizing CH 2 Br + and CH 2 I + as well as CH 2 Br + /CH 2 I + . The insets in both cases show the respective Hessian matrices for the ratio objectives exhibiting strong diagonal features with local off-diagonal couplings.
H A B nearly mutually commute with all the commutator matrix elements being ∼2 orders of magnitude smaller than the respective matrix elements in the pair of commuted Hessians. This result partially arises from the strongly diagonal aspect of all the Hessian matrices, but as noted the off-diagonal elements evident in figures 3(c), (d) and 7(a) were particularly reduced in magnitude upon commutation. Consistent with the near commutation of the Hessians, their associated eigenvectors were also similar (not shown). Therefore, this collective information leads to the conclusion that the landscape maximal domains for CH 2 Br + , CH 2 I + and CH 2 Br + / CH 2 I + are similar in location and structure. In contrast, the difference spectra from + H Br and + H I does not match that from + + H Br I (not shown), which represents the likely common situation of the control for the optimal ratio being different from the optimal control of either Br + or I + along with accompanying differences in the local landscape structures. Further insights may be anticipated from evaluating the Hessians + H Br and + H I at the optimal field of Br + /I + instead of their respective optima in figures 3(a) and (b), which will be investigated in the future. Although the control mechanisms in the ratio cases may differ from those for the single fragment yields, the same CH 2 BrI system is under control. Therefore, the features in figures 7(a) and (b) are consistent with the general findings for the individual spectra in figure 4.

Models and simulations
In order to gain insight into the processes underlying the experimental results in section 3, here we apply LHA to some simple idealized systems. The main goal is to explore the physical origin of the key Hessian features as described in section 3. Importantly, detailed modeling of photo-induced molecular fragmentation requires knowledge of multiple electronic potential energy surfaces involving higher lying electronic excited states and their couplings as well as the performance of nuclear dynamics [46][47][48][49]. Such simulations can be arduous, and we aim to demonstrate that even a simple physical model containing resonant and non-resonant processes can qualitatively explain the observed Hessian spectra. This result is important as it suggests that future analyses should be able to create specific physical models with varying levels of detail that reliably describe the dynamics, possibly by drawing in a suitable level of vibronic information about the system and fitting the LHA spectra to appropriate model parameters. As limiting reference cases, section 4.1 treats a pure resonantly coupled few level quantum system, and section 4.2 deals with a purely non-resonant multiphoton process. Finally, section 4.3 examines the mixed case containing both resonant and non-resonant transitions in the former two cases. The qualitative similarities between the experimental LHA eigenvalue spectra in section 3 and the latter case suggests that, at least in some situations, the optimal control of laser-induced dissociative ionization may be modeled in terms of non-resonant multiphoton absorption processes coupled to resonant transitions.

Resonant few level systems
The application of LHA to resonantly coupled few level quantum systems [12,13] is illustrated here through numerical simulations for a case similar in form to atomic Rb [16,50]. The model in figure 8(a) is a four-level system with free Hamiltonian H 0 and dipole moment μ. We denote by 〉 k | the eigenvector of H 0 associated to the eigenvalue ϵ k , = … k 1, , 4. The Hamiltonian is given by and the state ψ 〉 t | ( ) evolution in the interaction representation is described by The only nonzero dipole matrix elements are μ μ = = * 1.59 21 12 , μ μ = = * 2.24 31 13 , μ μ = = * 0.48 42 24 , and μ μ = = * 0.54 43 34 , where e is the electric charge. The observable signal J is taken as the population in state 〉 |4 at time T after the pulse is over In the simulations, equation (4.1) was solved in the rotating wave approximation and ψ 〉 t | ( ) was propagated using sequential steps of size Δt by An optimum field was first identified using a GA to maximize J. The spectral amplitude profile shown in figure 8(b) with a fluence of 0.8 − mJ cm 2 was used in all the numerical simulations in section 4. The optimal pulse temporal profile in figure 9(a) has structure spanning −1000 to 1000 fs, clearly distinct from the narrow transform-limited pulse (∼30 fs FWHM) associated with the spectral bandwidth in figure 8(b). The Hessian in figure 9(b) and its eigenvalues in figure 9(c) were determined by the procedure in section 2. The Hessianʼs clear diagonal features at ∼780 nm and ∼794 nm are associated with the 〉 → 〉 |1 |2 and 〉 → 〉 |1 |3 transitions, respectively. Also a weak diagonal peak appears at ∼776 nm corresponding to the transition 〉 → 〉 |2 |4 . There is no distinctive feature associated with the 〉 → 〉 |3 |4 transition, and this situation likely reflects the particular identified optimal solution amongst many others [12,13,16,50]. The nondiagonal Hessian elements are associated with couplings between the spectral components of the field [16] with clear features at wavelengths corresponding to interference between the two lowest order pathways (i.e., the feature at 780 nm, 795 nm). There is also an expected coupling feature at 776 nm, 780 nm reflecting the transitions of pathway 1 in figure 9. Three outstanding non-zero eigenvalues are evident in figure 9(b) along with three additional small eigenvalues, for a total of six, which is consistent with the prediction [14] of there being up to 2N-2 eigenvalues for the N = 4 levels in figure 8(a).

Non-resonant  -photon absorption
Consider the electronic transfer from a ground state 〉 |0 to an excited state 〉 |1 through a  -photon transition (see inset of  Hessian treatment for this process is given in appendix where it is shown that the Hessian eigenvalues are expected to form a continuous set as shown in figure 10. The latter spectra were generated from equation (A.7) with  = 6, ω  10 = 9.3 eV, using the statistical sampling procedure in section 2. The spectrum of the simulated six photon process in figure 10 is similar to the TPA experimental spectrum of figure 1(a). Continua eigenvalue spectra were observed to form a portion of all the experimental data in section 3, however the simultaneous existence of outstanding discrete Hessian eigenvalues in some of the experimental cases suggests that the mechanisms can not be fully explained in terms of non-resonant multiphoton absorption processes. The continuum portion of the experimental spectra in section 3 also show quantitative distinctions from the form in figure 10. Such differences likely lie in the detailed nature of the non-resonant processes, and this circumstance provides a basis to build future models with measured levels of structure in keeping with the Hessian features and the resultant spectra. The following section explores the Hessian matrix arising from a model that contains both resonant and non-resonant transitions.

Non-resonant  -photon absorption coupled to resonant excitations
Atomic and molecular transitions often take place over a broad range of energies much larger than the limited laser bandwidth available in the laboratory to induce the associated transitions. In this circumstance, laser experiments frequently exploit non-resonant multiphoton excitations along with any available single photon resonant transitions [33]. As shown in section 4.1 a system with a few resonantly coupled levels exhibits a discrete Hessian spectrum, while section 4.2 shows that a pure non-resonant multiphoton absorption process produces a continuum Hessian spectrum. In this section, we apply LHA to a system containing a combination of resonant and non-resonant transitions. The goal is to explore the impact of such models on the Hessian structure and distinctive spectral features, and in doing so we aim to provide insights into the experimental results of section 3 and thereby support Figure 10. Eigenvalues corresponding to a  -photon non-resonant absorption process driven by a transform-limited pulse corresponding to the spectrum shown in figure 8(b). The field was parametrized into 80 spectral intervals and a corresponding number of Hessian eigenvalues. The inset schematically depicts a  = 6 photon transition. the use of LHA as a simple experimental tool reflecting the underlying resonant/non-resonant nature of laser controlled dynamics. In performing this analysis there are a broad variety of models in the category of resonant and non-resonant processes. For example, there may be several contributing non-resonant transitions of different orders, the non-resonantly coupled levels may appear at various locations with respect to the resonantly coupled levels, and the detailed coupling within and between the resonant and non-resonant transitions can all affect the outcome in a quantitative fashion. In this regard, the model below is one amongst many that were considered, and we found that the qualitative character of the Hessian spectrum (i.e., its distinctive discrete-continuum feature) was weakly dependent on the model variation. However, the quantitative character of LHA is sensitive to the model details, forming a basis for future model identification from the experimental data.
The energy level diagram and transitions chosen for the model are shown in figure 11(a). To determine the state of the system ψ 〉 t | ( ) , equation (4.1) was solved numerically with    [51]. A simpler model with three levels (e.g., 〉 |0 , 〉 |1 , and a final target state 〉 |2 resonantly coupled to state 〉 |1 was also explored with LHA and yielded qualitatively similar results as the ones shown below. In relation to the experiments of section 3, states 〉 |0 and 〉 |1 may be interpreted as corresponding to CH 2 BrI and CH 2 BrI + , respectively. State 〉 |1 is then also coupled to states 〉 |2 and 〉 |3 through one-photon resonances and then to state 〉 |4 representing a particular fragmentation product. Beyond these comments, the model here is explored for its essential Hessian features, and as such it was not pursued for refinement to fit the experimental data. As shown in the appendix section A.1, the optimum phase for the non-resonant  -photon transition 〉 → 〉 |0 |1 corresponds to the pulse with maximum intensity (i.e., the transformlimited pulse). On the other hand, the optimum production of population in state 〉 |4 requires the use of a shaped pulse that deviates from being transform-limited as shown in figure 11(b). The optimal pulse has a dominant intense near-bandwidth-limited central peak and low intensity pre-and post-pulse features. The Hessian and its eigenvalues and eigenvectors were determined (see figures 11(c) and (d), respectively) using the procedure of section 2. In contrast to the Hessian for the resonantly addressable few level system (see figure 9(b)), the Hessian for the non-resonant/resonant model studied here shows stronger diagonal features whereas the off-diagonal features resemble the resonant case. The set of extracted eigenvalues reflects the combination of non-resonant and resonant situations present in the model with three outstanding eigenvalues along with a continuum of small-amplitude eigenvalues. There are also clearly defined rapid R and slow S varying domains (not labeled in figure 11(d)), as with the experimental Hessian in figure 4. The eigenvectors corresponding to the three outstanding eigenvalues show features centered at the resonant transition frequencies. The prominent off-diagonal Hessian structure was very well reconstructed (see section 2) using only the first three eigenvectors along with their corresponding eigenvalues. Within the scope of the simple model here, the LHA is consistent with all of the experimental cases in section 3; the three outstanding eigenvalues will disappear for the situation of the model reducing to that of being purely non-resonant as in section 4.2.
The simulated population dynamics of the model in figure 11 is shown in figure 12. The fieldʼs temporal profile E t | ( )| is also shown again for comparison in figure 12(a) (doted line). The high-intensity peak around t ∼ 0 in E t | ( )| induces the 〉 → 〉 |0 |1 transition through multiphoton absorption, and is also responsible for a small part of of the population transfer from 〉 |1 to 〉 |2 and 〉 |3 , while the relatively low intensity post-pulse structure in E t | ( )| for > t 0 completes the transfer from 〉 |1 to 〉 |2 and 〉 |3 and then transfers almost all of the population from 〉 |2 and 〉 |3 to 〉 |4 . The overall low yield in state 〉 |4 is likely a result of the limited fluence and the use of only phase modulation during the optimization. The low-intensity pre-pulse was found to have no effect in the population dynamics.
The detailed dynamics can be highly complex for laser driven dissociative ionization of polyatomic molecules. However, the outstanding features of the LHA spectra in section 3 and the simulations here suggest that these systems can be effectively described in terms of simplified physical pictures reflected in their key Hessian spectral features. As shown above, the distinctive LHA eigenvalue structures apparent in the experimental data can be qualitatively reproduced by a simple model that combines non-resonant and resonant dynamical transitions. Thus, considering the qualitative similarities between the experimental and model LHA features, we propose that the controlled dynamics arise from the existence of resonantly addressable energy levels coupled to a dominant non-resonant process for Br + , I + , and CH 2 I + , while optimization of CH 2 Br + appears to not exploit a resonant pathway. No attempt was made to introduce a suitable model for the optimization ratio of two products, but the results suggest that an expansion of the present model with competing final product states would likely accommodate that circumstance.

Conclusion
We introduce a general experimental LHA tool, LHA, with detailed descriptions of its implementation, which may be routinely employed in conjunction with AFC experiments to aid in elucidating the underlying control mechanisms driven by an optimal pulse. Performing LHA only calls for taking additional statistically sampled data about the optimal field followed by post facto analysis. The extracted Hessian matrix along with its eigenvalues and eigenvectors can provide valuable hints regarding the nature of the processes associated with the optimal control along with insight about the dynamical complexity of the system under control. It may be especially valuable in assessing the control of complex systems when the potential identification of a simple underlying model can be of high value. Additional applications are expected to provide experience in mechanistically interpreting the results of LHA laying out a clear foundation for that development.
We applied LHA to control mechanism assessment and dynamic dimensionality identification of a complex quantum system involving the dissociative ionization of the polyatomic molecule CH BrI I and + + Br I were selected. As shown in sections 3 and 4, the findings of the experimental studies and the modeling work demonstrate that the complexity of the underlying system dynamics can be reduced to a relatively simple physical picture with the differences in the mechanisms for each product reflected in their Hessian structure. In particular, our experimental results combined with theoretical analysis and the use of suitable models suggest that these nominally complex control problems can be interpreted in a simple mechanistic fashion involving a combination of nonresonant and resonant transitions. Beyond these key features additional subtle quantitative distinctions in the data provide a basis for future model refinement building on the present parsimonious treatment in section 4. Current work is underway to combine LHA with judicious vibronic information about a particular system to fit the Hessian data to practical models. Future works may benefit from systematically performing LHA on a series of optimal solutions at different laser energies to reveal information about the nature of the resonance transitions combined with model extraction from the data. The latter developments go beyond the scope of the present work and additional study is needed to establish LHA as a routine practical tool. Moreover, LHA is not limited to single observable objectives or the simple ratios discussed in this paper, as it can also be used for other complex objective fitness functions suitable for the particular needs. This new tool has its roots in the understanding of control landscape structure, which opens up applications of a wide ranging nature in quantum control where mechanistic insights are always valuable to ascertain.
where ω , and ϵ i is the energy of state 〉 i | . Let us also assume the existence of at 1, non-resonant with the field and itʼs first −  1 harmonics, i.e., ω ≫ This situation can be encountered, for instance, when all the 〉 k | levels lie well above 〉 |1 , i.e., ϵ ϵ ≫ k 1 , ∀ k, although the derivation below is not restricted to this case. With these assumptions, the  th order Dyson expansion term of the evolution operator is given by   Notice that the integral over −  t 1 can not be solved as the previous nested integrals, since ω ω ≈  o 10 and, thus,  E t ( ) oscillates as fast as ω t exp (i ) 10 and cannot be considered a constant. Equation (A.5) expresses the fact that a  -photon non-resonant transition from 〉 |0 to 〉 |1 is governed by the Fourier transform of  E t ( ) at frequency ω 10 [19]. A measurable signal resulting from multiphoton absorption would be proportional to Let us consider the landscape for  -Photon absorption in which an electron from the ground state 〉 |0 is excited to state 〉 |1 through the interaction with a laser field E(t). The fitness function of this landscape is given by where Ω E ( ) k is the Fourier transform of E t ( ) k at frequency Ω, and for simplicity we take α  = 1. In the second line of equation (A.7) it was assumed that → +∞ T and ⩽ = E t ( 0) 0 (i.e., the field is turned on at t = 0). If the carrier frequency of the laser field E(t) is ω o , then we have that ω ω ≈  o 10 for  -photon absorption. The variation of Ω  E ( ) with respect to the fieldʼs spectral phase ϕ ω ( ) can be calculated from   Moreover, since for large  we have that ω ω ω − ≈ 10 10 , then the Hessian eigenvalues for the  -photon process form a continuous spectrum whose dependence on ω is approximately proportional toω E Re ( ( )). Figure 10 shows eigenvalues corresponding to  -photon absorption and for the field ω E ( ) employed in the simulations of sections 4.1 and 4.3.