A convenient method to simulate and visually represent two-photon power spectra of arbitrarily and adaptively shaped broadband laser pulses

We report a convenient model of second harmonic (SH) spectra resulting from shaped laser fields created in pixilated pulse shaping devices based on a matrix of all possible pair-wise, field interactions. This model is useful in estimating simultaneous two-photon absorption (TPA) and second harmonic generation (SHG) cross sections. A simple geometric interpretation of the model allows for a graphical representation of the SH spectrum as a set of vectors in the complex plane. This interpretation can be used to illustrate the effects of applied phase (and amplitude) shaping on second-order observables and it provides a facile explanation of the spectral focusing mechanism observed in numerous control problems.


Introduction
Adaptive laser pulse shaping, an idea first introduced by Judson and Rabitz [1], offers rich opportunities in photophysical and photochemical research because it provides a robust methodology to manipulate complex chemical systems whose excited-state behavior is not addressable by state-of-the-art control theory. Several recent reviews [2]- [6] (and the references therein) offer a glimpse into the power of this technology. However, it remains difficult to extract control mechanisms and other salient physical information from the results of adaptive control experiments. This is an important focus of current research in the community.
The interpretation of adaptive optimization results can arise when the heuristic learning methodology does not distinguish between a 'good' pulse feature and those which contribute negligibly to the fitness [7]- [9]. Unless the converged control fields exhibit obvious temporal features [10]- [19] or there is a well-formed theory suggesting a suitable parameterization [20]- [23], the complexity of optimal pulse shapes can obscure the physical or chemical interpretation of a control mechanism. There has been some elegant work to minimize this discrepancy [9,24]. However, great care must be taken in any parameterization of an optimization experiment to avoid introducing experimental bias. Such bias may ultimately limit what is learned from a technique that otherwise permits and exploits pulse-shape complexity.
One class of control problems that have been extensively explored is concerned with the manipulation of second-order responses in atomic and molecular systems. Two-photon absorption (TPA) and second harmonic generation (SHG) have been investigated in a number of control scenarios [25]- [33]. The development of perturbation theory approaches to this problem builds on the work of the Noordham group [34] and the Silberberg group [35,36]. In some cases, perturbation theory of second-order responses as it pertains to pulse shaping applications can be simplified by analytic approximations such as Taylor series expansion of the applied phase [27]. However, the very nature of adaptive pulse-shaping makes the use of truncation techniques in theoretical models risky as even routine experiments can involve essentially arbitrary forms of the applied phase (and/or amplitude). In our view, descriptions are needed that both provide an intuitive understanding of how phase and amplitude functions become manifest in the shape of the second harmonic (SH) power spectrum and are capable of capturing the possible highdimensionality of actual experiments. These should emphasize patterns that herald control.
With this in mind we report a simple model of SH spectra resulting from shaped laser fields created in pixilated pulse shaping devices based on a matrix of all possible pair-wise field interactions. This is useful in estimating simultaneous TPA and SHG cross sections. A simple geometric interpretation of the model allows for a graphical representation of the SH spectrum as a set of vectors in the complex plane. This interpretation can be used to illustrate the effects of applied phase (and amplitude) shaping on second-order observables and it provides a facile explanation of the spectral focusing mechanism observed in numerous control problems.

Results and discussion
In the absence of any real intermediate, the observed yield, Y , of a second-order process in a molecule or material system resulting from simultaneous TPA can be expressed as Y ∝ g (2) TPA ( )S (2) Here, g (2) TPA ( ) is the two-photon excitation spectrum and S (2) ( ) is the SH power spectrum of the laser. In perturbation theory, g (2) TPA ( ) is a property of the system and in the case of idealized SHG it is assumed to be a constant function of frequency. The quantity S (2) ( ) is given by (2) In the above expression, is the frequency dependent phase and E(ω)is the Fourier transform of the electric field E(t). In this paper, lower case Greek letters (e.g. ω) are used where the frequency range corresponds to the fundamental laser pulse and upper case Greek letters (e.g. ) are used where the frequency range corresponds to second-order observables (i.e. 2ω). Equation (2) states that the intensity of the SH spectrum at a given frequency is determined by the square modulus of the integral of the product of all combinations of the electric field (i.e. E(ω m ) and Eω n )) whose frequencies add to (i.e. = ω m + ω n ). It is straightforward to use equations (1) and (2) to numerically simulate the second-order response of an idealized system given arbitrary amplitude and phase inputs [27,29]. However, from a conceptual standpoint this formulation is somewhat unsatisfying. The spectral response at any one frequency depends on the interaction of phase and/or amplitude shaping inputs over the range of frequencies ω and it is generally difficult to intuit how a particular SH spectrum results from phase and/or amplitude inputs. This is especially true when the optimal fields from adaptive control experiments have complex phase shapes.

Matrix-based approximation
Our laboratory employs a home built pulse shaper in the 4f-geometry of a zero dispersion compressor utilizing a pixilated liquid crystal spatial light modulator (SLM) to input phase and/or amplitude functions. In an effort to draw connections between experimental control variables and the mechanisms exploited by adaptive algorithms or open-loop phase and/or amplitude shaping tools, we have found it advantageous to develop a mathematical model where the amplitude and phase functions are indexed according to the SLM pixels. This is most easily achieved by approximating the frequency separation between pixels (|ω n+1 − ω n | = ω) as constant, where n corresponds to pixel number (in a typical laboratory experiment 100-200 pixels are used). Although the 4f-geometry of a pulse shaper does not disperse the pulse linearly in frequency at the Fourier plane, this approximation is reasonable when the bandwidth is only tens of nanometers in the near IR. In a typical experiment in our lab, the measured ω deviates from the average value of ω by a maximum of ∼3% over the full-width at half-maximum (FWHM) of the laser pulse. In making this assumption, the Jacobian is treated as constant, which has the effect of slightly exaggerating the contributions of the redder frequencies to the SH spectrum. However, for the simulations herein, the effects are negligible given the bandwidth of our laser.
The utility of this approximation is that it leads to a discrete parameterization of the electric field (in frequency) indexed by the pixilation of the SLM, where ω n , A n and φ n are the central frequency, applied amplitude and applied phase at pixel n, respectively. Note the input range for A n is [0, 1] and the input range for φ n is [0, 2π ]. As . E 4 (ω 4 ,φ 4 ) E 5 (ω 5 ,φ 5 ) . E 5 (ω 5 ,φ 5 ) Figure 1. Representation of the sums s 3 and s 5 as dictated by the operator SDE 3 [M] (see equation (5)) and SDE 5 [M] for a hypothetical five-pixel SLM. All elements in the i = 3 set (connected by blue arrows) consist of two field components that sum to the same total frequency 3 . This is easily seen by imagining that the matrix element index numbers are equal to the frequency of the individual field components: 3 + 1 = 2 + 2 = 1 + 3 = 4 = 3 . Similarly, all elements in the i = 5 set (connected by red arrows) consist of two field components that sum to the same total frequency 5 . Note, here each A n = 1 and amplitude is not explicitly written.
stated, the intensity of the SH spectrum at a given frequency is determined by the square modulus of the integral of the product of all combinations of the electric field (i.e. E(ω m ) and E(ω n )) whose frequencies add to (i.e. = ω m + ω n ). Given our discrete parameterization, all possible nonzero product combinations of the electric field can be easily written by specifying a matrix M (see figure 1), and each matrix element can be readily calculated using equation (3). Because pixel index is proportional to frequency (a good approximation with our bandwidth, vide supra), all elements in sets that are 'perpendicular' to the matrix diagonal involve two electric field components whose individual frequency values sum to the same total frequency i . The index i refers to the set of 'perpendicular' matrix elements m mn . For example, when i = 1 the set consists only of element m 11 , when i = 2 the set consists of elements m 12 and m 21 , when i = 3 the set consists of elements m 13 , m 22 and m 31 , etc. Note there are 2n − 1 of these sets, where n is the number of SLM pixels. Because each of the elements in a set i consists of components of the electric field whose individual frequencies sum to the same total frequency i , equation (2) can be approximated with the matrix-based formulation: Here, is an operator we define to sum each element of the ith set. The operator name 'SDE' is chosen to refer to 'Sum of Degenerate Elements'. Essentially, SDE i [M] is a Riemann sum approximation to the integral in equation (2). To help clarify this, a schematic representation of the sums s 3 and s 5 as prescribed by the operator SDE 3 [M] and SDE 5 [M] for a hypothetical five-pixel experiment are shown in figure 1. The blue and red arrows represent the direction of the summation. Using this methodology the yield of a second-order process can be readily approximated, where g (2) TPA ( i ) is the value of the two-photon excitation spectrum at i . This matrix-based approximation could easily be extended to third-(or higher) order processes by calculating all possible field interactions with a three-dimensional (or higher) array and summing over the appropriate elements.

Visual interpretation
For our purposes, the matrix-based approximation of the SH spectrum is useful because it admits a simple geometric interpretation. A visual representation in the complex plane can provide an intuitive feel for how phase and/or amplitude shaping of the laser pulse affect the SH spectrum. This has allowed us to better understand various amplitude and phase switching methodologies that will be described elsewhere. In order to produce the visual representation in the discussion that follows, certain simulation parameters are fixed. Herein, a Gaussian function G n (ω n ) is used to define the electric field amplitude (|E n (ω n )| = G n (ω n )) of the fundamental at each pixel. The Gaussian is parameterized such that the input pulse is centered at 800 nm with an FWHM of ∼30 nm. We set n = 128 and define a spectral window {ω 1 , ω 2 , . . . , ω n } → {λ 1 , λ 2 , . . . , λ n } between 760 and 840 nm. This results in i = 255 sums of degenerate elements from the matrix M, which defines the spectral window { 1 , 2 , . . . , i } → { 1 , 2 , . . . , i } between 380-420 nm of the second-order response. In the discussions that follow, only the effects of phase shaping are considered. Therefore, A n = 1 for all pixels and A n is not included in subsequent equations.
Each matrix element m mn within M defines a vector in the complex plane, whose magnitude, |m mn |, is determined by the product |m mn | = |E m E n | and whose argument is determined by the angle θ mn = (φ m + φ n ). Because each m mn of M represents a vector in the complex plane, the operation SDE i [M] can be thought of as adding vectors 'head-to-tail' to make a new vector whose magnitude, when squared, is proportional to the expectation value of the SH spectrum of the laser pulse at i . Note, although amplitude shaping is not treated here, it would affect the expectation value of the SH spectrum by modulating the magnitude of |m mn |. Conversely, phase shaping affects the expectation value of the SH spectrum by modulating the direction of the component vectors m nm relative to each other. If one considers the case in which i = 128 given the simulation parameters above, then s 128 = SDE 128 [M] is determined by the sum of 128 elements (i.e. the set m mn corresponding to 128 ). The amplitude of this vector squared defines the spectral intensity at the center of the SH spectrum ( i ∼ 400 nm). It is noted that this visual interpretation is similar to the one used by Mitra and Rabitz in their work on Hamiltonian encoding. They have shown that vectors in the complex plane are useful for illustrating when different quantum pathways linking an initial state with a final state are in phase leading to constructive interference or out of phase leading to destructive interference [37]- [40].

Transform limited pulses.
In the case of a bandwidth limited pulse with zero applied spectral phase, the argument of the complex exponential in equation (7)  As a simple comparison, a bandwidth limited pulse where the applied phase at each pixel is π/2 is considered. In this case, θ mn = π everywhere, and each vector element is rotated 180 • in the complex plane and points in the negative direction along the real line. Figure 2(b) shows this case. Because all of the vectors are still in phase, the magnitude of s 128 that would result is unchanged and the expectation value of the SH is also unchanged. While this is a simple example, it illustrates that the addition of a constant phase causes all the vector elements to rotate through the same angle, which results in a macroscopic rotation of the vector sum as a whole. It is noted that in the previous two examples, all vector sums s i (including s 128 ) point in the same direction; e.g. in the positive real direction when the applied spectral phase at each pixel is zero and in the negative real direction when the applied spectral phase at each pixel is π/2.
It is well known that the application of a linearly increasing or decreasing phase (as a function of frequency) does not alter the temporal envelope of a laser pulse, bandwidth limited or otherwise. Therefore, application of linear phase does not affect the SH spectrum. With this in mind, it is interesting to consider how the visual representation of the vector sums s i = SDE i [M] illustrates this phenomenon. Here three cases are presented based on linear phase functions of the form φ(n) = a · n + b, where n is the pixel number (increasing with increasing wavelength), a is an arbitrary scalar and b is an arbitrary constant.
In Because of this phase angle, each component vector of s 3 is rotated through the same angle in the complex plane. As a result, the vector sum s 3 is rotated relative to s 5 , but s 3 is still of maximum possible amplitude. Analogously, the phase angle of the matrix element in s 1 (m 11 ) is four times the negative value of the scalar a and the resultant vector sum is rotated through an angle twice the size of the angle in s 3 . In this way, each successive vector sum s i is rotated through a constant angle relative to s i−1 although the magnitude of the sum remains constant. Consequently, the rotation of resultant vector sums s i relative to one another (resulting from the application of a linear phase) has no bearing on the SH spectrum. The first case ( figure 3(a)) should be compared with the second and third cases (figures 3(b) and (c)). In the second and third cases, b = 0 so that there is a phase offset at the center of the spectrum. In figure 3(b), a = 0.02, which is the same as in the example discussed above. In figure 3(c) a 'steeper' linear phase is used (a = 0.03). In both cases figures 3(b) and (c), s 128 does not point in the positive real direction because there is a residual phase offset at the center of the spectrum. Additionally, the relative rotation of one vector away from another is greater in figures 3(c) than (b) because of the steeper linear phase. In both cases, the vector sums, s 128−i and s 128+i , fan outwards symmetrically from s 128 . In each case, the relative angle between consecutive vectors is constant and the magnitude of the angle depends on the value of the scalar a; but again, this has no bearing on the observed SH spectrum.

Chirped pulses.
The addition of linear chirp or quadratic phase is the simplest applied phase function that results in an observable change in the SH spectrum at i . Now the individual vectors m mn are rotated relative to each other in the complex plane, shortening the resultant vector s i = SDE i [M]. Consider, for example, a sequence of five phase functions of the form φ(n) = a · n 2 + b, where n is the pixel number (increasing with increasing wavelength), a is a scalar used to increase the chirp and b is an arbitrary constant. Figure 4 shows the arrangement of vectors m mn head-to-tail (red) and the resultant vector s 128 (blue) as linear chirp is increased from left to right. In each of these plots b = 0. It can be seen in each of the plots that quadratic phase or linear chirp causes the individual vectors to rotate with respect to one and the other (i.e. to get out of phase) leading to tighter spiraling of the individual vectors plotted head to tail. Because the rotation leads to dephasing of the component vectors, the magnitude of the sum necessarily decreases. As a result, the expectation value of the SH spectrum at 128 decreases as well. Note, the addition of any nonlinear applied phase of the form φ(n) = a · n m + b, where m is an integer greater than 1, will also necessarily change the observed SH spectrum (and s 128 ).
It is common to think of phase functions representing orders of chirp as even or odd polynomials relative to the window of pixels. This can be achieved easily using φ(n) = a · (n − c) m + b such that the quantity |(n − c)| has the same value at the first and last pixel. In this case, some care must be taken in considering s 128 because individual vectors will only rotate with respect to one and the other when m is an even integer thereby changing the SH intensity at 128 . When m is an odd integer (greater than 1) s 128 and therefore the SH intensity at 128 will not change because for each matrix element in the sum, as was observed with the matrix elements of s 128 in the case of an applied linear phase, the argument of the complex exponential is constant [21,35,36]. Of course for either even or odd (>1) m, the total SH intensity will decrease as a is increased.
From a conceptual standpoint, it proves useful to connect shaping parameters (such as a and b above) to the various qualities of component vectors m mn and the vector sums s i . In the case above of an applied quadratic phase φ(n) = a · n 2 + b, (b = 0), there are two distinct features that can be conceptualized as 'motions' (here as function of a). Firstly, there is a 'rotation' of the vector sum s i as a whole, which results from the fact that the numerical first derivative of the phase function (i.e. the difference between the applied phase at consecutive pixels) changes as a is varied. As seen in figure 4 (comparing panels from left to right) the resultant blue vector s 128 rotates clockwise as a in increased. This is analogous to the rotation of the vector sum from the application of a linearly increasing phase in figure 3. Note that the rotational motion does not impact the expectation value of the SH spectrum. Note also that one can always observe the total phase θ mn = (φ m + φ n ) accumulated (modulo 2π) for any element m mn of the sum s i = SDE i [M] by identifying the direction of the individual vector m mn (one of the red head-to-tail components) in the complex plane. Secondly, there is a 'decay' of the vector sum magnitude, which results from the fact that the numerical second derivative (i.e. the difference between consecutive values of the numerical first derivative) also changes as a is varied. This rotates the component vectors with respect to each other. It is the decay that accounts for observed changes in the expectation value of the SH spectrum. Note in the simple case of a linear chirp (quadratic phase), the normalized rate of decay (as a function of a) is the same for all vectors s i because the numerical second derivative of each applied phase function is constant. Consequently, application of a linear chirp (as a function of a) causes the SH spectrum to decrease homogeneously. However, by applying phase functions that contain higher order terms (as a function of some scalar) it is possible to vary the normalized rate of decay of the various vectors s i , which can be used to shift the maximum value of the SH spectrum to the red or the blue (see section 2.3 below). Additionally, as will be shown in a forthcoming paper, a description of the resulting vector sum properties, in the case where the applied amplitude and phase have been parameterized to produce phase-locked pulse pairs, provides useful insight into various applications of Fourier transform spectroscopy.

Ratio optimization
In previous work, several groups have demonstrated that adaptive femtosecond pulse shaping optimizations can improve the measured yield of emission from a molecule (following TPA) relative to SHG (the ratio TPA/SHG) by exploiting a two-photon spectral focusing mechanism [25,27,29]. In these experiments, the adaptive algorithm exerts control by finding pulse shapes that are best able to focus (or shift) the maximum intensity of the SH spectrum into the spectral region where the TPA cross section of the molecule is greatest. In this way, the algorithm converges upon pulse shapes that improve the ratio (TPA/SHG) over a homogeneously distributed SH spectrum because the conversion efficiency of SHG is ostensibly constant over the SH spectral window. Using the matrix-based approximation, it is straightforward to simulate a ratio optimization experiment (TPA/SHG) and to plot SH spectra during the optimization. Furthermore, the visual representation of the vector sums s i and the component vectors m mn demonstrate how the various components of the shaped electric field combine to focus the SH spectrum into the appropriate spectral range to effect control.
For the following simulations, we calculate SHG and TPA based on an electric field parameterized by E n (ω n , φ n ) = |E n (ω n )|e iφ n , where |E n (ω n )| = G n (ω n ) (i.e. the same Gaussian envelope used above). To simulate SHG, we assume that the conversion efficiency is uniform over the SH spectral range. To simulate TPA, we use a linear function of frequency to approximate a simple molecular TPA spectrum; namely, g (2) TPA ( i ) = β i + γ , where β and γ are arbitrary constants chosen such that g (2) TPA ( 1 ) = 1.0 and g (2) TPA ( 255 ) = 2.0 (in arbitrary units). A hybrid evolution-strategy/genetic-algorithm [29] is used to optimize the simulated feedback (TPA/SHG) over 300 generations of 60 individuals. Because the magnitude of the arbitrarily chosen TPA spectrum is greatest at shorter wavelengths (higher frequencies), the control mechanism can be understood phenomenologically in terms of the algorithm shifting and focusing the SH spectrum toward the blue edge of the SH range of the fundamental. As an example, the SH spectra (derived from the matrix-based approximation equation (5)) of the best pulses from generation 0, 25, 50 and 299 from a typical optimization are shown in figure 5 (black lines). For comparison, the SH spectrum of a bandwidth limited pulse is shown in figure 5 (gray dashes; off scale). In the arbitrary units, the SH spectrum of a bandwidth limited pulse has a maximum intensity of 1.0 at ∼400 nm. Clearly the SH spectrum focuses to the blue edge of the spectrum relative to the red to maximally accommodate the TPA spectrum while minimizing SHG. The visual interpretation of the matrix approximation provides an insightful look into the underlying physics. The individual vectors contributing to s 89 , s 90 , s 91 , s 92 , s 93 and s 94 are shown for generations 0, 25, 50 and 299 in figure 6. The component vectors for these sums correspond to the SH intensity in the spectral region ∼393.9 to ∼394.6 nm and were chosen because of their proximity to the maximum intensity of the SH spectrum of the optimized pulse (see figure 5). For comparison, the vectors contributing to s 128 for generations 0, 25, 50 and 299 are shown in figure 7. The component vectors for this sum correspond to the SH intensity at ∼400 nm and were chosen because they are illustrative of the spectral region into which the SH spectrum is not focused.
These results demonstrate visually how the manipulation of the applied phase causes different SH intensities at the various SH frequencies.  figure 7. In this latter case, as the experiment progresses, it is evident that the algorithm is able to worsen the alignment of the component vectors during the optimization. In other words, successive vectors tend to point in different directions because of the higher order terms contained in the applied phase on the red side of the spectrum. This misalignment of the component vectors translates to a smaller vector sum and the expectation value of the SH spectrum decreases at these frequencies over the course of the optimization. In other cases (not shown) the algorithm does not obviously worsen the alignment of the component vectors in the spectral region where the SH intensity is not focused, it simply does not meaningfully improve the SH intensity above what is achieved by a randomly chosen pulse. In totality, the alignment of the component vectors on the blue side of the spectrum coupled with a decrease (or no corresponding improvement) in alignment on the red side of the spectrum leads to the spectral focusing mechanism.
As noted earlier, our visual interpretation resembles one used by Mitra and Rabitz in their work on Hamiltonian encoding [37]- [40]. In this example shown above involving spectral focusing, the similarity is clear in the sense that both use alignment of vectors in the complex plane to show how multiple pathways between starting and end points can constructively or destructively interfere by way of achieving control. We can confirm that the head-to-tail vector sums are illuminating the expected control mechanism by consideration of the phase functions explored and exploited by the adaptive algorithm. Compare the unwrapped phase functions and the numerical second derivatives (i.e. the difference between the difference of the applied phase at each pixel) of the best pulses from generations 0 and 299, which are shown in figure 8. (GEN 0 left; GEN 299 right; phase functions top panel; second derivatives bottom panel). It is evident that the algorithm achieves spectral focusing by manipulating the applied phase over the course of the optimization such that it is nearly linear (i.e. the second derivative approaches zero in generation 299) on the blue side of the spectrum, while maintaining higher order terms (i.e. the second derivative is still nonzero and rapidly oscillating in generation 299) on the red side of the spectrum. This result follows naturally from the discussion in the previous sections, where it was observed that the application of linear phase did not attenuate the SH spectrum, whereas the application of nonlinear phase did attenuate the SH spectrum. These results are consistent with the experimental observations reported by the Joffre group [27]. We note that the second derivative of the phase function is substantially nonzero at the very blue edge of the SH range; however, because the intensity of the electric field is relatively small over this range, this portion of the spectrum does not contribute significantly to the observed fitness. Therefore, the algorithm does not meaningfully manipulate these pixels over the course of an experiment.

Conclusions
In conclusion we have developed new tools to represent and facilitate understanding of how broadband laser pulses that have been phase-(and/or amplitude-) shaped in pixilated SLM devices can affect nth-order perturbative phenomena. Second-order non-resonant processes are the focus of the current work and the methodology is based on a matrix of pair-wise products between field components in the fundamental laser spectrum. We have defined an operator that sums elements 'perpendicular' to the matrix diagonal where each of these elements (of a given sum i) corresponds to the components of the electric field whose individual frequencies sum to the same total frequency i in the SH spectrum of the shaped laser field. In this way, we readily generate a vector in the complex plane whose magnitude is proportional to the twophoton intensity at that frequency i . More importantly, we can visualize component vectors representing each pair-wise pathway to i in the SH spectrum. The patterns that emerge in the head-to-tail arrangement of these vectors illuminate the accumulation of phase in achieving a particular two-photon power at i . In a simple simulated adaptive control experiment, we have shown that these patterns of head-to-tail vector components illuminate how control is being achieved by the heuristic algorithm; namely, second-harmonic power spectrum focusing. In future work this methodology will be exploited in experiments on complex systems that exhibit complex adaptively discovered optimal phase functions. We hope to identify frequencies where multiphoton response is important and where particular phase-dependence is being exploited. Another direction of interest is to generalize the methodology to consider multiphoton absorption cases where one-photon resonances are energetically accessible in the fundamental spectrum. It would be significant if patterns emerge in the head-to-tail arrangement of vector components that herald phase jumps in the vicinity of a resonance [41]- [45]. This would help draw mechanistic connections between closed-loop adaptive experiments and open-loop phase control. Such connections would be valuable interpreting adaptive control in systems being driven with weak to strong fields [45,46] where both perturbative and dynamic Stark effects can be present.