text:
# Non-Gaussian Resistance Noise near Electrical Breakdown in Granular Materials

## 1 Introduction and Model

Non-Gaussian distributions of several quantities characterizing the behavior of complex systems in non-equilibrium steady-states, have been evidenced since long time in many experiments . The general properties of these anomalous distributions and their link with other properties of non-equilibrium complex systems is far from being fully understood. Rather, the study of these distributions is nowadays attracting an increasing interest in the literature . Here, we study the distribution of resistance fluctuations of a conducting thin film with granular structure near electrical breakdown. Electrical breakdown phenomena typically occur in conductors stressed by high current densities and they consist in an irreversible and dramatic increase of resistivity of the material and thus they are associated with a conductor-insulator transition . In our study we make use of the Stationary Network Under Biased Percolation (SNUBP) model . This model provides a good description of many features associated with the electrical instability of composites materials  and with the electromigration damage of metal lines , two important classes of breakdown phenomena. The film is modeled as a resistor network which reaches a steady state determined by the competition between two biased stochastic processes, breaking and recovery. Systems of different sizes and with different levels of internal disorder are considered, where the disorder is mainly related to the fraction of broken resistors within the network. The resistance and its fluctuations are calculated by Monte Carlo simulations which are performed under different stress conditions. Resistance fluctuations are found to deviate from Gaussianity near electrical breakdown. We analyze and discuss this non-Gaussianity in the spirit of the Bramwell-Holdsworth-Pinton (BHP) distribution, recently introduced in the context of the study of highly correlated systems near criticality . In particular, we show that the deviations from Gaussianity observed near the breakdown are related to finite size effects, and consequently vanish in the thermodynamic limit. Nevertheless, near the percolation critical point the non-Gaussianity persists in the large size limit and it is described by the universal BHP distribution.

According to the SNUBP model, a conducting film with granular structure is described as a two-dimensional resistor network. Precisely, we consider a square-lattice of $`N\times N`$ resistors. This lattice lies on an insulating substrate at a given temperature $`T_0`$, which acts as a thermal bath. Each resistor can be in two different states: (i) regular, corresponding to a resistance $`r_n=r_0[1+\alpha (T_nT_0)]`$ and (ii) broken, corresponding to the effectively “infinite” resistance, $`r_{OP}=10^9r_n`$ (in the following we call this state defect). Here, $`\alpha `$ is the temperature coefficient of the resistance and $`T_n`$ the local temperature. This latter is determined by taking into account Joule heating effects and thermal exchanges between neighbor resistors : $`T_n=T_0+A[r_ni_n^2+(3/4N_{neig})(r_li_l^2r_ni_n^2)]`$, where $`i_n`$ is the current flowing in the nth resistor, $`N_{neig}`$ the number of nearest neighbors over which the summation is performed. The parameter $`A`$ represents the thermal resistance of each resistor and determines the importance of Joule heating effects. By taking the above expression for $`T_n`$ we are neglecting time dependent effects in heat diffusion studied by Sornette et al. . A constant stress current $`I`$ is then applied through perfectly conducting bars at the left and right sides of the network. The two biased processes consist of stochastic transitions between the two possible states of each resistor and they occur, through thermal activation, with probabilities: $`W_{Dn}=exp[E_D/k_BT_n]`$ and $`W_{Rn}=exp[E_R/k_BT_n]`$, characterized by the two activation energies, $`E_D`$ and $`E_R`$  ($`k_B`$ is the Boltzmann constant). The network time evolution is obtained by a Monte Carlo simulation which updates the network resistance after breaking and recovery processes, according to the iterative procedure described in details in Ref. . The sequence of successive configurations provides a resistance $`R(t)`$ signal after an appropriate calibration of the time scale. Then, depending on the stress conditions ($`I`$ and $`T_0`$) and on the network parameters (sizes, activation energies and other parameters related to the material like $`r_0`$ and $`\alpha `$), the network reaches a steady state or undergoes an irreversible electrical failure. This last possibility is associated with the achievement of the percolation threshold, $`p_c`$, for the fraction of broken resistors. Therefore, for a given network at a given temperature, a threshold current value, $`I_B`$, exists above which electrical breakdown occurs . Below this threshold, the steady state of the network is characterized by fluctuations of the fraction of broken resistors, $`\delta p`$, and of the resistance, $`\delta R`$, around their respective average values $`<p>`$ and $`<R>`$. In particular, we underline that the ratio $`(E_DE_R)/k_BT_0`$ determines the average fraction of defects for a given value of $`I<I_B`$ and thus the level of disorder inside the network. In the following section we analyze the results of simulations carried out by considering networks of different sizes, with different levels of disorder and stressed by different currents. In all cases we take $`T_0=300`$ (K), $`E_D=0.170`$ (eV), $`r_0=1`$ ($`\mathrm{\Omega }`$), $`\alpha =10^3`$ (K<sup>-1</sup>), $`A=5\times 10^5`$ (K/W) (these values are chosen as physically plausible). $`N`$ ranges between $`50÷125`$, while $`E_R`$ between $`0.0259÷0.164`$ (eV).

## 2 Results

We report in Fig. 1 the resistance of a $`75\times 75`$ network as a function of time and at increasing currents. The evolutions $`R(t)`$ in Fig. 1 are obtained by taking $`E_R=0.103`$ (eV) for the activation energy of the recovery process. This value of $`E_R`$ provides a network with an intermediate level of disorder. This figure displays two important features of the electrical response of a conducting film: i) the linear regime at low currents is followed by a nonlinear regime where the average resistance increases significantly with the current; ii) the amplitude of the resistance fluctuations increases strongly close to the breakdown. Precisely, the lowest curve in Fig. 1 corresponds to the linear regime, i.e. it is obtained for $`I=0.001(A)<I_0`$ where $`I_0`$ is the current value associated with the onset of the nonlinearity of the I-V characteristic . The second curve corresponds to the nonlinear regime, $`I=0.45(A)>I_0`$, the third to the threshold for the breakdown, $`I=I_B=0.70(A)`$, and finally the highest curve to a network undergoing an electrical breakdown at a current $`I=0.75(A)>I_B`$. This last curve has been shifted upwards by $`0.1(\mathrm{\Omega })`$ for graphical reasons. A detailed analysis of the behavior of the average resistance and of the relative variance of resistance fluctuations as a function of the current can be found in Ref. .

We have investigated the Gaussianity of the steady state signals shown in Fig. 1. The results of the analysis are reported in Fig. 2 for $`I=0.001`$ (A) (linear regime) and in Fig. 3 for $`I=I_B=0.70`$ (A) (threshold current for breakdown). Precisely, in these figures we show on a lin-log plot the product $`\sigma \mathrm{\Phi }`$ as a function of $`(<R>R)/\sigma `$, where $`\mathrm{\Phi }`$ is the probability density function (PDF) of the distribution of $`\delta R`$ and $`\sigma `$ is the root mean square deviation from the average. This normalized representation, by making the distribution independent of its first and second moments, is particularly convenient to explore the functional form of the distribution. These results, together with all the others presented in this paper have been obtained by considering time series containing about $`1.2\times 10^6`$ resistance values. Moreover, for comparison, we have also reported the Gaussian distribution (which in this normalized representation has zero mean and unit variance) and the Bramwell-Holdsworth-Pinton (BHP) distribution . These authors have proposed a universal non-Gaussian PDF for the fluctuations of global quantities of systems at criticality. Denoting $`m`$ a fluctuating quantity (for example the magnetization of a ferromagnet like in Ref. ), $`<m>`$ and $`\sigma _m`$ its mean value and root mean square deviation respectively, $`P(m)`$ its PDF, $`y(m<m>)/\sigma _m`$ the normalized variable, $`\mathrm{\Pi }(y)\sigma _mP(y)`$ the normalized PDF and $`xb(ys)`$, the BHP distribution takes the following expression :

$$\mathrm{\Pi }(y)=K[e^{xe^x}]^a$$
(1)
where $`a=\pi /2`$, $`b=0.936\pm 0.002`$, $`s=0.374\pm 0.001`$ and $`K=2.15\pm 0.01`$. This expression can be considered as a generalization of the Gumbel distribution, which is often associated with the occurrence of rare events. In Figs. 2 and 3 together with the PDF obtained for networks with linear sizes $`N=75`$ (small circles) are also displayed the PDF corresponding to $`N=50`$ (stars), $`N=100`$ (squares) and $`N=125`$ (down triangles). Figure 2 shows that in the linear regime the distribution of $`\delta R`$ is Gaussian for all system sizes. Near breakdown, at $`I=I_B`$, (Fig. 3) the PDF shows non-Gaussian tails. However, this non-Gaussianity progressively vanishes when the system size is increased. Therefore, when networks with intermediate level of disorder are considered, as in the case of Figs. 2 and 3, deviations from Gaussianity are related to the finite size of the system. The PDF in this case is well fitted by the Eq. (1) once the parameters $`a`$, $`b`$, $`s`$ and $`K`$ are considered as fitting parameters, as described in Ref. . Such phenomenological fits can provide helpful tools in the study of failure precursors in the case of finite size systems.

We have also considered the role of disorder in the breakdown process and its effect on the distribution of resistance fluctuations. At given values of $`E_D`$ and $`T_0`$, the parameter controlling the disorder inside the network, i.e. the average fraction of broken resistors $`<p>`$, is the activation energy $`E_R`$. Figure 4 shows the average fraction of broken resistors versus the external current for different values of $`E_R`$. The data are obtained for a $`75\times 75`$ network. The different sets of symbols refer to steady states of the network. The region above the solid curve corresponds to electrical breakdown (in this region the stable states of the network are non-conducting states). Thus, Fig. 4 represents the phase diagram of the system in the $`<p>I`$ plane. It can be shown  that all the curves in Fig. 4 collapse onto a single one after normalization of $`I`$ to $`I_0`$ and of $`<p>`$ to $`<p>_0`$ (the defect fraction in the limit of vanishing current). Moreover, the relative variation of the defect fraction, $`[<p><p>_0]/<p>_0`$, scales as the ratio $`I/I_0`$ with a quadratic exponent . Furthermore, we note that the electrical breakdown corresponding to the values of $`E_R`$ considered in Fig. 4, is associated with a first order phase transition . These results agree with the behavior observed in electrical breakdown experiments, performed in the Joule regime of composites . Nevertheless, it can be shown  that when $`E_R`$ reaches its maximum value, $`E_{R,MAX}`$, (determined by the stability condition obtained for a system of given size in the vanishing current limit ), the conductor-insulator transition becomes of the second order. This change in the nature of the transition, when going from small to high disordered systems, also agrees with the predictions of Andersen et al. . Therefore, we have reported in Fig. 5 the normalized PDF calculated at the breakdown threshold for different values of $`E_R`$. Again, a $`75\times 75`$ network is considered. We can see that the non-Gaussianity of the distribution at $`I=I_B`$ increases significantly when increasing the $`E_R`$ value. However, as already shown in Fig. 3, this non-Gaussianity vanishes in the large size limit. Nevertheless, when the value of $`E_R`$ is nearly equal to $`E_{R,MAX}`$, the system approaches the critical point and the PDF tends to achieve the BHP form, independently of system size. This is shown in Fig. 6 which displays the normalized PDF calculated at the breakdown threshold, $`I=I_B`$, for $`E_R=0.164`$ (eV). The three sets of data are obtained for $`N=75`$, $`I=0.0065`$ (A) (crosses), $`N=100`$ and $`I=0.009`$ (A) (down triangles) and $`N=125`$ and $`I=0.011`$ (A) (full circles) and correspond to networks near their critical point (the data for $`N=50`$ are not shown in Fig. 6 because $`50\times 50`$ networks are always unstable for this value of $`E_R`$). We can see that in this case the PDF becomes independent of the sizes of the system and it is well fitted by the BHP distribution.

## 3 Conclusions

We have studied the distribution of the resistance fluctuations of conducting thin films with granular structure near electrical breakdown. The study has been performed by describing the film as a resistor network and by using the SNUBP model . We have considered systems of different sizes and with different levels of internal disorder. A non-Gaussianity of the fluctuation distribution is found near electrical breakdown. This non-Gaussianity increases with the degree of disorder of the network. However, we have shown that this non-Gaussianity is related to the finite size of the system and that it vanishes in the large size limit. Nevertheless, near the critical point the non-Gaussianity persists in the large size limit and is well fitted by the universal Bramwell-Holdsworth-Pinton distribution .


format: text

text:
# Extracting quantum dynamics from genetic learning algorithms through principal control analysis

## Abstract

Genetic learning algorithms are widely used to control ultrafast optical pulse shapes for photo-induced quantum control of atoms and molecules. An unresolved issue is how to use the solutions found by these algorithms to learn about the system’s quantum dynamics. We propose a simple method based on covariance analysis of the control space, which can reveal the degrees of freedom in the effective control Hamiltonian. We have applied this technique to stimulated Raman scattering in liquid methanol. A simple model of two-mode stimulated Raman scattering is consistent with the results.

preprint: A FOCUS Preprint

The central challenge of coherent control of quantum dynamics is to find the optimal path to guide a quantum system from its initial state to some target final stateTannor and Rice (1988); Brumer and Shapiro (1992). Several theoretical methods have been developed to aid this search Tannor et al. (1986); Peirce et al. (1988), and there has been considerable experimental success as wellShnitman et al. (1996); Gordon (1995). However, in all but the simplest systems, the search is hampered by incomplete knowledge of the system Hamiltonian. Strongly coupled systems such as large molecules in condensed phase are so complicated that it is nearly impossible to calculate optimal pulse shapes in advance.

Feedback learning algorithms overcome this limitation by using the physical system itself to explore its own quantum dynamics through an experimental search Judson and Rabitz (1992). A typical search experiment compares the ability of several thousand different shaped laser pulses to transform the system $`|\psi `$ from its initial state at time $`t=0`$ to some desired target state $`|\chi `$ at a later target time $`t=T`$. Examples of transformations that have been studied include molecular photodissociation, atomic photoexcitation and photoionization. The pulse shapes are selected through a fitness-directed search protocol, such as a genetic algorithm Holland (1992). The fitness is a measured quantity proportional to the objective functional $`J[H;x_i]`$, which is the square of the projection of $`|\chi `$ onto $`|\psi `$ at the end of the experiment:

$$J[H;x_i]=|\chi (T)|\psi (T)|^2.$$
(1)
$`J`$ depends on the Hamiltonian $`H`$ for the system evolution, which depends in turn on the laser electric field $`E(t)`$ determined by the settings of the $`n`$ pulse shape control parameters $`x_i`$, $`i=1\mathrm{}n`$. $`J`$ reaches its extreme value for the optimal pulse. This pulse can be calculated using optimal control theory if $`H`$ is known Peirce et al. (1988); otherwise, it must be discovered through the learning search.

Several recent papers have suggested modifications or extensions of learning feedback that can measure properties of the system Hamiltonian Geremia et al. (2001); Mitra and Rabitz (2003); Daniel et al. (2003); Geremia and Rabitz (2002). Here we propose a different approach based on analysis of the trial experiments. We will show that the ensemble of trial pulse shapes can reveal essential features of the dynamics.

Genetic algorithms and similar evolutionary search strategies have many different variations Davis (1991). Our implementation starts with approximately 50 randomly generated optical pulse shapes produced by spectrally filtering an ultrafast laser pulse Tull et al. (1990). Each pulse shape is described by a column matrix of control parameters $`x_i`$ called a genome consisting of about 25 numbers (genes), each encoding the amplitude and/or phase of a different segment of the optical spectrum. The control target is measured for each pulse shape. Then the algorithm creates a new generation of pulse shapes by combining attributes of the fittest members of the previous generation Pearson et al. (2001). After several generations, the pulse shapes usually cluster near high fitness regions of the search space. When the algorithm finds a pulse shape or several shapes that cannot be improved over many generations, the search stops, and the highest fitness pulse shape is declared the optimal solution to the search. We test 1000 to 10,000 pulse shapes in a typical experiment. We maintain a record of every pulse shape, its fitness, and its parentage (genealogy).

The learning algorithm achieves control without prior knowledge of the system Hamiltonian, and has far more degrees of freedom $`n`$ than the minimum required for control. The number of possible solutions is exponential in $`n`$. In a typical search, the phase of each color is adjusted by the spectral phase filter to a precision of about ten degrees, so there are $`2^5`$ possible values of each gene. This means that the number of possible solutions for a genome of length 25 is $`2^{5\times 25}4\times 10^{37}`$. Genetic algorithms can search this large state space with great efficiencyHolland (1992). Unfortunately, simply finding a good solution has not often provided significant insight into the system dynamics or Hamiltonian. The optimal pulse shape found by the learning algorithm, while sufficient to achieve control, is often complicated and may contain unnecessary features.

The conditions for reaching an extremum in $`J[H;x_i]`$ may only depend on two or three essential features of the control field $`E(t)`$. These features are not obvious in the successful genome because they may depend on all 25 genes. The Hamiltonian could be written in a much simpler form if these essential degrees of freedom ($`u_j`$) were found. Here we show how to establish the $`u_j`$ through covariance analysis of the pulse shapes evaluated during the learning searchJolliffe (2002). Covariance analysis is commonly used to reduce the dimensionality of and to find patterns in high dimensional data sets.

We propose to apply these techniques, not to sets of data, but to the control space for the experiment. Linear combinations of genes with high fitness should appear correlated in the fitness-driven genetic algorithm. These correlated linear combinations correspond to the principal components of the control field that direct the quantum dynamics under investigation. Principal control analysis is the application of covariance techniques to a fitness directed search.

Principal control analysis is implemented on our system by calculating the covariance matrix of the set of all pulse shapes in the search, defined by:

$$C_{ij}=\delta _i\delta _j\delta _i\delta _j,$$
(2)
where the expressions $`\delta _i=x_{i+1}x_i,i=1\mathrm{}n1`$ are the nearest neighbor phase differences. By using the phase differences in the analysis we remove the ambiguity associated with the unimportant global phase. The covariance matrix is not the only measure of correlation. We could also weight the terms of the covariance matrix using the fitness, or normalize each term to the individual gene variance. In this paper, we will use the simple covariance. This is appropriate because all the $`\delta _i`$’s are of the same type.

Once the covariance matrix is determined, we calculate its eigenvectors and eigenvalues. Each eigenvalue $`\lambda _j`$ measures the variance of the projections onto the corresponding eigenvector. This has a special meaning for a learning control search: it shows how far the control setting moved during the learning process. A small subset of eigenvalues usually contains most of the weight of the trace of the covariance matrix. The correlations of the projections with the pulse shape fitnesses allow us to determine which control directions were most important for the physical process under consideration. The controls expressed in the basis of the eigenvectors are uncorrelated: Each of these controls changes the fitness without correlation with the others over the search set.

We propose that these eigenvectors with the largest fitness correlation are the essential control directions ($`u_j`$). We expect the eigenvectors with the larger eigenvalues to be most strongly correlated with the fitness of a pulse shape solution. Conversely, a low correlation indicates those eigenvectors that have not contributed substantially to increasing the fitness during the search. These eigenvectors correspond to extraneous dimensions, which could be eliminated (i.e., their projection set to zero) without losing substantial control.

By projecting the GA solutions onto the $`k<n`$ eigenvectors that correlate best with the fitness (the principal controls), we reduce the dimension of the control space. The solutions with highest fitness, when expressed in the reduced basis of the principal controls, represent the essential features of the search solutions. The objective functional also takes on a simpler form in this basis:

$$J=J[H;u_{1,\mathrm{},k}]$$
(3)
where each $`u_j`$ can assume a range of values on the order of $`\pm \sqrt{\lambda _j}`$. The search for the specific target state is now a matter of optimizing each control $`u_j`$ over this range.

In summary, we propose to apply covariance techniques to the control space of learning feedback experiments. Control values derived from the genomes are analyzed by a covariance matrix, defined in Eq. 2. The matrix eigenvectors are independent control directions. The correlation of the fitness with the eigenvectors suggests which control directions are the most important. The corresponding eigenvalues indicate the necessary excursion along each axis. Therefore, searches conducted in the eigenvector basis are more efficient. Finally, the best solutions are found by projecting the optimal learning control solutions onto the important control directions.

We now apply this analysis to a well-studied control problem: the selective excitation of vibrational modes in liquid methanol. The experiment has been described previously Pearson and Bucksbaum (2004); Pearson et al. (2002, 2001). An intense shaped 800 nm ultrafast laser pulse (the pump laser) is focused into a cell containing methanol. Above a threshold fluence, the pump induces stimulated Raman scattering into either the symmetric or antisymmetric Raman-active C-H stretch mode. Either mode can be selectively excited by adjusting the shape of the pump pulse through phase shaping and/or amplitude shaping of its spectrum.

Figure 1 depicts the results of a phase-only feedback control experiment. Panel a (b) shows the Raman spectra produced by the optimal pulse for the symmetric (antisymmetric) stretch mode. This learning search included $`2720`$ different pulse shapes. Moderate fitness increases were observed for either target mode after $`25`$ generations. The fitness increase was greater for the symmetric mode. This is typical of our searches based on phase-only controlPearson et al. (2002).

A Wigner representation of the pulse shape solution found by the learning algorithm is plotted in panel c (d). The Wigner function is a spectrally resolved field auto-correlation:

$$W(\omega ,t)=𝑑\omega ^{}E(\omega \omega ^{})E^{}(\omega +\omega ^{})e^{2i\omega ^{}t}$$
(4)
Wigner representations are complete time-frequency spectrograms of the optical control field (up to a global phase), but the important features leading to control of the methanol are obscure. The inability to interpret the result easily is typical of many GA search solutions Levis et al. (2001); Brixner et al. (2001).

The principal control analysis of this problem begins with a single covariance matrix (Eq. 2) for the entire population of pulse shapes evaluated in the two feedback experiments. The physical system under control was the same in the two problems; only the target was different. We therefore expect the independent searches to be nearly spanned by a small number of eigenvector controls. The covariance matrix is not simple to interpret because the principal components in this problem are widely distributed among all of the control settings in the search space. However, the essential features of the control problem begin to emerge if the covariance matrix is diagonalized. This can be seen in Fig. 2a, where the eigenvalues of the covariance matrix are plotted in descending order. The algorithm has made large excursions only along the few control directions that have large eigenvalues. The phase functions associated with the three largest eigenvalues are shown in Fig. 2b. These control directions are also the most strongly correlated with fitness. Figure 2c (d) plots the correlation of eigenvector with fitness for the symmetric (antisymmetric) stretch. The three most significant eigenvectors are the same as the control directions that correlate most strongly with the fitness. Therefore we propose that the dimension of the search space can be reduced without inhibiting control.

Each pulse shape can now be re-expressed in the eigenvector basis. To arrive at the essential features of the best pulse, we calculate the projections $`\eta _k`$ of the optimal pulse shape onto the principal control directions $`u_k`$. When these $`u_k`$ are expressed in the original basis, their components are individual discrete frequencies that make up the field. Therefore, $`_{j=1}^k\eta _ju_j(\omega )`$ produces a pulse that contains traits necessary to achieve the target, with minimal extraneous features. Figure 1e (f) shows the Wigner plot of this essential pulse for the symmetric (antisymmetric) stretch mode. The essential pulse for the symmetric stretch preserves $`66\%`$ of the original pulse shape vector ($`_{j=1}^k\eta _j^2(\omega )=66\%`$). For the antisymmetric stretch, this number is $`88\%`$.

This procedure is independent of the specific nature of the physical system or the dynamical Hamiltonian; however, the principal control directions can be used to construct a simplified interaction Hamiltonian, since the laser electric field now only depends on a few parameters:

$$H(t)=H(E[u_{j=1\mathrm{}k};t]).$$
(5)
This can provide important constraints. For example, a recent paper demonstrated a control mechanism for SRS in methanol based on periodicity in $`I(t)`$ Pearson and Bucksbaum (2004). The Fourier transform of $`I(t)`$ then reveals the most important Raman coupling frequencies in the problem Meshulach and Silberberg (1998). Figure 3 (top) shows the magnitude of the $`FT[I(t)]`$ for the optimal pulse shape found by the learning algorithm for the symmetric stretch mode (dashed) and for the same pulse shape projected onto the principal control directions $`u_k`$ (solid). Figure 3 (bottom) shows similar plots for the antisymmetric mode.

For both modes, projecting the optimal pulses onto the principal control directions enhances the frequency components around $`3`$ THz. The coupling frequency found through principal control analysis agrees with the model based on the mode separation in methanol. Additionally, comparing the phases of the Fourier transforms for each mode yields a phase difference of $`\pi /2`$ in the region around $`3`$ THz, consistent with the model described in reference Pearson and Bucksbaum (2004).

In conclusion, we have shown how covariance analysis of a genetic search algorithm can uncover essential features of the dynamical Hamiltonian. Although our example involved only phase-shaping of the optical field, this technique should be applicable to any system where fitness-directed learning algorithms have been used to reveal the path from an initial quantum state to a target. Principal control analysis can also be incorporated into the experimental search protocol. By discovering the principal controls, it should be possible to search the space more efficiently, and to test ideas about the system dynamics as the search is proceeding. The method should be most useful in cases where the dynamics can be described by only a few principal degrees of freedom, which are linear combinations of the control parameters of the search space.

We thank Anna Amirdjanova, Jayson Cohen, Vladimir Dergachev, Chitra Rangan, Ben Recht, and Tom Weinacht for valuable discussions. This work was supported by the National Science Foundation under grant 9987916.