## Abstract

We show through numerical simulations and experimental data that a fast and simple iterative loop known as the Fienup algorithm can be used to process the measured Maker-fringe curve of a nonlinear sample to retrieve the sample’s nonlinearity profile. This algorithm is extremely accurate for any profile that exhibits one or two dominant peaks, which covers a wide range of practical profiles, including any nonlinear film of crystalline or organic material (rectangular profiles) and poled silica, for which an excellent experimental demonstration is provided. This algorithm can also be applied to improve the accuracy of the nonlinearity profile obtained by an inverse Fourier transform technique.

© 2004 Optical Society of America

## 1. Introduction

Until very recently, the second-order optical nonlinearity spatial profile of thin films was measured by the Maker fringe (MF) technique [1]. In this classical technique, a fundamental laser beam is focused onto the nonlinear sample and the intensity of the second-harmonic (SH) signal generated within the nonlinear region is measured as a function of the incidence angle of the fundamental beam. This process produces a signal, the MF curve, which is proportional to the square of the Fourier transform (FT) magnitude of the sample’s nonlinearity profile *d*(*z*). However, since the phase of this Fourier transform is not known the Fourier transform cannot be inverted uniquely to obtain *d*(*z*). To overcome this difficulty, we have recently demonstrated a number of techniques where two nonlinear samples, either identical or different, are pressed against each other and the MF curve of this sandwich is measured [2,3]. As a result of interference between the fundamental and SH signals within the two samples, the MF curve now includes information not only about the magnitude but also the phase of the Fourier transform. With straightforward processing of this MF curve it is now possible to recover the complete complex FT (amplitude and phase), which is then inverted to retrieve the profile *d*(*z*) of the two samples unambiguously. The concept of forming a sandwich structure has also been used to estimate the spatial width of the nonlinearity profile of poled silica samples [4]. However, this approach assumed a uniform nonlinearity profile and did not involve a FT formalism.

In this letter, we present a new, much simpler and more accurate technique to uniquely recover a nonlinearity profile from the classical MF measurement of a *single* sample instead of sandwich structures. This novel technique makes use of an iterative loop known as the Fienup algorithm [5]. We show that for a broad range of profiles, this algorithm enables the accurate recovery of the missing FT phase and thus of the nonlinearity profile. This approach constitutes a substantial improvement over previous techniques because of the simplicity of both the measurement and the computer code, and because of the greater speed of the data processing. The error in the recovered profile is also reduced. Furthermore, we show that this same iterative technique can be used to improve the accuracy of the profiles retrieved by any of the IFT techniques [2,3].

## 2. Iterative processing of d(z) using the Fienup algorithm

The Fienup algorithm [5] is one of several iterative algorithms that have been proposed to recover a function from partial information about this function in the space (or time) and/or Fourier domain [6]. This algorithm uses the known (i.e., measured) FT magnitude of some unknown function *f*(*t*), together with the known properties of this function, such as being real or causal, and iteratively corrects an initial guess for *f*(*t*). In the context of the present work, we use it to recover the nonlinearity spatial profile *d*(*z*) of a thin (up to hundreds of microns) nonlinear sample from the measured MF curve of this sample. Note that without loss of generality we can assume that *d*(*z*) is a real and causal function, i.e., *d*(*z*)=0 for *z*<0, where *z*=0 defines the edge of the nonlinear wafer. The MF measurement of a single sample provides the magnitude of the spectrum |*D*_{M}
(*f*)| of the FT of *d*(*z*), where *f* is the spatial frequency, [3] but it does not provide the phase spectrum *ϕ*(*f*) of this FT. The role of the Fienup algorithm is to recover this missing FT phase. To start, one makes an initial guess, *ϕ*
_{0} (*f*), for the unknown FT phase. We have found empirically that this initial guess does not strongly impact the convergence of the algorithm, so for convenience it can simply be *ϕ*
_{0} (*f*)=0. As illustrated in Fig. 1, the first step in the algorithm is to calculate numerically the IFT of the complex quantity |*D*_{M}
(*f*)exp(*jϕ*
_{0}(*f*)). The second step is to take the real part of this IFT and retain only the *z*≥0 region, which gives a function *d*
_{1}(*z*) that constitutes a first estimate of the nonlinearity profile. The third and final step is to compute the FT of this profile, |*D*
_{1}(*f*)|exp(*jϕ*
_{1}(*f*)). The phase *ϕ*
_{1}(*f*) of this FT provides a new estimate for the missing FT phase *ϕ*(*f*). At this point the FT of *d*(*z*) has a known (measured) amplitude |*D*_{M}
(*f*)| and a best-estimate (calculated) phase *ϕ*
_{1}(*f*). The previous three-step cycle is then repeated using *ϕ*
_{1}(*f*) instead of *ϕ*
_{0}(*f*) as the new FT phase, which yields a second estimate *d*
_{2} (*z*) for the profile and a second estimate of the FT phase *ϕ*
_{2}(*f*). This process is iterated *n* times until convergence is achieved, i.e., until the average difference between the profiles *d*
_{n-1}(*z*) and *d*_{n}
(*z*) obtained during two consecutive cycles is less than a preset value, for example 0.1%. What the algorithm has done is reconstruct a more accurate spectrum *ϕ*_{n}
(*f*) than the initial guess for the originally unknown FT phase.

The Fienup algorithm has a number of interesting properties that have a strong bearing on its applicability to the present problem. First, although convergence has not been proved rigorously, we have found empirically that the algorithm always converges, as also pointed out by others [6]. Second, the Fienup algorithm converges to the correct solution for an extremely wide range of profile shapes. This is illustrated in Fig. 2, which shows a few standard profile shapes (buried Gaussian, rectangular profile, exponential, etc.) and the functions retrieved by the processing algorithm. With all of these profiles the Fienup algorithm works exceedingly well and recovers a profile that is extremely close to the original profile, in fact so close that in Fig. 2 the retrieved profiles cannot be distinguished from the original profiles. The average error in the recovered profile is as less than 0.004%, except for the rectangular profile, for which the error is 0.008%. These accuracies were achieved after only 100 iterations, which took only a few seconds on a 500-MHz computer. Note that these profiles were selected because they are physical nonlinearity profiles either measured or predicted by various theoretical models in nonlinear thin films of glasses, polymers, and crystals. Rectangular profiles, in particular, are a common occurrence in nonlinear crystalline films such as LiNbO3 and nonlinear organic materials.

To understand why these particular profiles (as well as countless others) are so well predicted by the Fienup algorithm, it is useful to take a closer look at the algorithm’s basic principle. The problem of retrieving a nonlinearity profile from the sole knowledge of its MF spectrum is analogous to the recovery of a real one-dimensional function from its FT magnitude alone. It is well known that, in general, the FT magnitude is not sufficient to recover the function. However, there are exceptions to this rule, namely families of functions for which the FT phase can be recovered from the FT magnitude alone, and visa versa [7,8]. One important such family is what is known as minimum-phase functions (MPFs) [8]. MPFs are characterized by having a *z*-transform with all its poles and zeros on or inside the unit circle. As a result of this property, the FT phase and the logarithm of the FT magnitude of an MPF are the Hilbert transform of one another. Consequently, the FT phase of an MPF can always be recovered from its FT amplitude, and an MPF can always be reconstructed from its FT magnitude alone. This reconstruction can of course be done by taking the Hilbert transform of the logarithm of the function’s FT magnitude to obtain the FT phase, then inverting the full (complex) FT. However, this direct approach is not the preferred solution because of difficulties in its implementation, in particular phase unwrapping.[6] A second approach is to use the Fienup algorithm: although it has not been proven mathematically, we and others [6] have found empirically that this algorithm always converges to a minimum-phase function. In other words, of the infinite family of FT phase functions that could be associated with the known FT amplitude, the algorithm converges to the one and only one that has the minimum phase. This duality between the Fienup algorithm and MPFs is the reason behind the outstanding agreement observed in Fig. 2: with the exception of the rectangular profile, all profiles in this figure happen to be minimum-phase functions.

In summary, the Fienup algorithm always converges to *the* minimum-phase function whose FT magnitude is |*D*_{M}
(*f*)|. Since this solution is unique, if we knew *a priori* that the profile to be reconstructed was an MPF, then we would be certain that the solution provided by the Fienup algorithm is the correct profile. Conversely, if the original profile is *not* a minimum-phase function, then strictly speaking the algorithm will not converge to the correct profile. Since in general we do not know whether the original profile is an MPF, it follows that we also do not know whether the recovered profile is the correct one.

In spite of this apparent limitation, we have observed through computer simulations that the Fienup algorithm converges to a solution surprisingly close to the correct profile for an extremely wide range of profile shapes. The reason for this broad success is that a large number of nonlinearity profiles happen to be either an MPF or close to (in some sense) an MPF. To understand intuitively which physical functions are likely to be a minimum-phase function, let’s denote an MPF by *d*_{min}
(*n*), where *n* is an integer that numbers sampled values of the function variable, e.g., distance *z* into the sample in the present case of a nonlinearity profile. Because all physical MPFs must be causal (although not every causal function is an MPF), *d*_{min}
(*n*) must be equal to zero in at least half of the space variation, for e.g., for *n*<0, as is the case of nonlinearity profiles. Second, the energy of a minimum-phase function, which is defined as $\sum _{n=0}^{m-1}$|*d*_{min}
(*n*)|^{2} for *m* samples of the function *d*_{min}
(*n*), must satisfy the inequality [8];

for all possible values of *m*>0. In Eq. 1, *d*(*n*) represents any of the functions that have the same FT magnitude as *d*_{min}
(*n*). This property suggests that most of the energy of *d*_{min}
(*n*) is concentrated around *n*=0 [8]. Stated differently, any profile with either a single peak or a dominant peak will be either a minimum-phase function or close to one and thus work extremely well with the Fienup algorithm. There could be other function besides functions with a dominant peak that are MPFs, but so many nonlinearity profiles, including the most common ones, satisfy this criterion, i.e., exhibit a dominant peak, that this sub-class of profiles is well worth investigating by itself. This sub-class of MPFs is at the origin of the broad success of the Fienup algorithm in this application. This criterion is obviously satisfied by the functions illustrated in Fig. 2, which again are all MPFs except for the rectangular profile. Although a rectangular profile is not truly an MPF, because it has a single peak it is expected to be close to an MPF, which is indeed the case (almost all the poles and zeros of its *z*-transform are on or inside the unit circle, and the remaining few are just outside). This explains the success of the Fienup algorithm with rectangular profiles (see Fig. 2): even though a rectangular profile is not a true MPF, it is close enough to one that the recovery is still amazingly good, with an average error of only 0.008%. The two-peak profile in Fig. 2 has a dominant peak (the secondary peak is only ~1/3 as high as the main one) and it is, as expected, a minimum-phase function. As a result, for this function the average error in the recovery is also extremely low (under 10^{-5} %). The nonlinearity profile of poled silica, which typically exhibits a sharp dominant peak just below the sample’s anodic surface [3], also satisfies the MPF criterion, which suggests that the algorithm should work well for poled silica, as will be demonstrated below.

To explore the robustness of this algorithm, we evaluated its accuracy for a number of functions other than minimum-phase functions. In a first series of simulations, uniform random noise (~14% peak-to-peak) was added to the profiles of Fig. 2, which were thus no longer MPFs. Each noisy profile was then subjected to the iterative loop of Fig. 1, assuming zero initial FT phase. In all cases, the profiles recovered after 100 iterations were in excellent agreement with the original noisy profile (average error under ~1.4%), as illustrated in Fig. 3 for an exemplary profile (buried Gaussian). The Fienup algorithm clearly works well even in the presence of noise.

In a second series of simulations, we investigated profiles with several peaks of comparable magnitude. Because none of the peaks is dominant, such profiles do not satisfy Eq. (1) and are not minimum-phase functions. We found that with two peaks of comparable magnitude the algorithm’s prediction is only marginally degraded and still quite acceptable. For example, when the peaks in the two-peak profile of Fig. 2 are given comparable height, the recovered profile is still essentially indistinguishable: the average error in the recovery increases to only 0.1%. With three peaks of comparable height (see Fig. 4) even though the Fienup algorithm is not nearly as accurate as in the previous examples, the recovered profile still provides a decent estimate of the original profile. These simulations demonstrate that the Fienup algorithm can be successfully applied to a wide range of profiles, and that it is certainly not limited to minimum-phase functions.

Another important observation about this recovery technique is that for all the functions we have tested, if the nonlinearity strength at the origin *d*(*z*=0) is increased to a value much larger than the rest of the profile, for example *d*(0)=5·max{*d*(*z*)}, the convergence of the algorithm improves substantially in terms of both accuracy and speed. This improvement is certainly expected because Eq. (1) is now satisfied, even if it were not satisfied prior to this change, which means that the function becomes an MPF. As an example, for the three-peaked nonlinearity profile of Fig. 4, if *d*(0) is increased to *d*(0)=10·max{*d*(*z*)}, the new recovered profile (red dashed curve in Fig. 4) is significantly closer to the original profile. This powerful result opens the possibility of recovering absolutely *any* nonlinearity profile by depositing on it a stronger and very thin nonlinear material, such as LiNbO_{3}. The film should be preferably thin, not because its thickness affects convergence (it does not) but because it is easier both to deposit and to measure its MF curve.

It is important to point out two minor limitations of this processing technique. First, it cannot recover the exact location of the profile within the sample, i.e., how deeply *d*(*z*) is buried below the surface of sample. Second, it cannot determine unequivocally the sign of the nonlinearity profile. Consequently, if *d*_{n}
(*z*) is the solution provided by the Fienup algorithm for a given nonlinear sample, then all ±*d*_{n}
(*z*-*z*
_{0}) functions are also solutions. However, in many cases these limitations are fairly inconsequential because being able to determine the shape of a nonlinearity profile is much more important than determining its sign or exact location. Furthermore, if need be both the sign and location of *d*(*z*) can be determined by other means, for example by using the two-sample IFT technique [3].

## 3. Application to experimental poled silica samples

To verify the applicability of this processing technique to practical nonlinear materials, we tested it against a wafer of poled silica. The sample, a 25×25×0.15 mm wafer of fused silica (Infrasil), was thermally poled in air at 270 °C and 4.8 kV for 15 min [9,10]. In a first step, the MF spectrum of this sample was measured and the Fienup algorithm was applied to it assuming zero initial FT phase to recover a first nonlinearity profile. In a second step, the same sample was characterized by the two-sample IFT technique [3], which provided a second, absolute nonlinearity profile. The hope was, of course, that the two profiles were the same. The setup used for these MF measurements can be found in earlier work [3,11].

Figure 5 shows the measured MF spectrum (open circles), which again is proportional to the FT amplitude |*D*_{M}
(*f*) |of *d*(*z*). The maximum achieved internal propagation angle (*θ*
_{max}) in this measurement is ~87°. However a large *θ*
_{max} of this sort is not a necessity, e.g., even *θ*
_{max}~78° is shown to supply adequate information for a fairly accurate recovery of the nonlinearity profile.[12] The FT phase and the profile recovered from this spectrum using the Fienup algorithm are plotted in Fig. 6 (green curve) and Fig. 7 (dotted curve), respectively. For comparison, on the same two figures we also plotted the FT phase (blue curve) and the profile (solid curve), respectively, recovered by the two-sample IFT technique. The profiles retrieved by the Fienup algorithm and by the IFT technique are in excellent agreement: both exhibit a sharp nonlinearity coefficient peak of magnitude *d*
_{33}≈-1.05 pm/V just below the surface of the sample, a sign reversal at a depth of about 12 µm, and a wider positive nonlinear region extending to a depth of about 45 µm. These observations are in keeping with other profiles obtained by other IFT techniques in similar poled samples (for further discussions about the physics of these observed nonlinearity profiles, please see Refs. 3 and 13–15). The FT phase spectra recovered by these two very different techniques (Fig. 6) are also in very good agreement with each other. These excellent agreements lend support to both our processing technique and the two-sample technique.

So far we have always assumed a zero initial FT phase *ϕ*
_{0}(*f*) when using our processing algorithm. Now that we have access to the FT phase recovered by the two-sample technique (blue curve in Fig. 6), we can run our processing algorithm using this spectrum instead of zero as a better initial guess for this FT phase. This operation is equivalent to using our algorithm to post-process the FT phase recovered by the IFT technique, with the hope to obtain an even more accurate nonlinearity profile. The profile obtained by this post-processing technique after 100 iterations is shown as the dot-dashed curve in Fig. 7. Comparison to the profile obtained with the two-sample IFT technique (Fig. 7, solid curve) shows that post-processing did not modify the overall profile shape, but that it significantly smoothed out the artificial oscillations introduced by the IFT technique. The post-processed profile of the two-sample technique (dot-dashed curve) and the profile obtained with the Fienup algorithm assuming zero initial phase (dotted curve) are very close to each other: the average difference between them is only 0.14%, which clearly demonstrates the validity of both approaches. The similarity between the profiles before and after post-processing confirms that the IFT technique came very close to recovering the actual profile. It also demonstrates the usefulness of our processing algorithm in another application, namely post-processing a profile obtained by an IFT technique. To better illustrate the beauty of this new post-processing technique, we show in Fig. 5 the MF spectrum of the post-processed nonlinearity profile, calculated numerically (solid curve). The agreement between this theoretical curve and the measured MF data (open circles) is very good.

We have tested this post-processing technique with our other two IFT techniques[2,3] and found that in all cases post-processing significantly attenuates unphysical oscillations in the profile. Post-processing is consequently a powerful tool to improve the accuracy of nonlinearity profiles recovered by an IFT technique. This processing step is also relatively fast: on a 500-MHz computer 100 iterations typically take only ~10 s, compared to 5–10 min. for the data processing in an IFT technique. Note also that the Fienup loop converges much faster if the thickness *W* of the nonlinear region is known. Knowing this thickness means that at the outset *d*(*z*) can be set to 0 not only over the *z*<0 space but also for *z*>*W*, which restricts the range of *z* values over which *d*(*z*) is unknown. Thus the iterative loop does not need to recover as many discrete values of *d*(*z*) and it converges faster.

## 4. Conclusions

We have shown through numerical simulations and experimental data that a fast and simple iterative processing loop known as the Fienup algorithm can be used to process the measured Maker fringe curve of a nonlinear sample to accurately retrieve both the shape and the magnitude of the sample’s nonlinearity profile. This algorithm works extremely well with profiles that are minimum-phase functions. For common minimum-phase functions such as a buried Gaussian or an exponential, the average error is less than 0.004% after 100 iterations. We also demonstrated that the algorithm converges to a solution very close to the exact profile not only for minimum-phase functions but also for noisy minimum-phase functions and more generally for any profile that exhibits a dominant peak, which covers a very wide range of practical nonlinearity profiles. This is true in particular of a rectangular profile, which is not a true minimum-phase function but does exhibit a dominant peak. The Fienup algorithm is therefore predicted to work well (average error of 0.008%) for deposited films of nonlinear crystals (e.g., LiNbO_{3}) or nonlinear organic materials. This is also true of profiles with two comparable peaks (average error of 0.1%) and even three-peak functions (error of a few %). Furthermore, the same algorithm can be applied to the nonlinear profile obtained by an IFT technique to improve the profile accuracy, in particular to remove unphysical oscillations introduced by the IFT technique data processing algorithm. The validity of these conclusions was verified experimentally with a nonlinear sample of poled silica. The profile recovered from the measured MF curve of this sample using the Fienup algorithm was found to be in excellent agreement with the profile of the same sample measured using an IFT technique.

## Acknowledgments

This work was supported by Litton Systems, Inc., a wholly owned subsidiary of Northrop Grumman.

## References and links

**1. **P. D. Maker, R. W. Terhune, M. Nisenhoff, and C. M. Savage, “Effects of dispersion and focusing on production of optical harmonics,” Phys. Rev. Lett. **8**, 21–22 (1962). [CrossRef]

**2. **A. Ozcan, M. J. F. Digonnet, and G. S. Kino, “Inverse Fourier transform technique to determine second-order optical nonlinearity spatial profiles,” Appl. Phys. Lett. **82**, 1362–1364 (2003). [CrossRef]

**3. **A. Ozcan, M. J. F. Digonnet, and G. S. Kino, “Improved technique to determine second-order optical nonlinearity profiles using two different samples,” Appl. Phys. Lett. **84**, 681–683 (2004). [CrossRef]

**4. **C. Corbari, O. Deparis, B. G. Klappauf, and P. G. Kazansky, “Practical technique for measurement of second-order nonlinearity for poled glass,” Electron. Lett. **39**, 197–198 (2003). [CrossRef]

**5. **J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. **3**, 27–29 (1978). [CrossRef] [PubMed]

**6. **T. F. Quatieri Jr. and A. V. Oppenheim, “Iterative techniques for minimum phase signal reconstruction from phase or magnitude,” IEEE Trans. Acoust., Speech, Signal Processing **29**, 1187–1193 (1981). [CrossRef]

**7. **M. Hayes, J. S. Lim, and A. V. Oppenheim, “Signal reconstruction from phase or magnitude,” IEEE Trans. Acoust., Speech, Signal Processing **28**, 672–680 (1980). [CrossRef]

**8. **V. Oppenheim and R. W. Schafer, *Digital Signal Processing*, (Prentice Hall, 2002), *Chap. 7.*

**9. **R. A. Myers, N. Mukerjkee, and S. R. J. Brueck, “Large second-order nonlinearity in poled fused silica,” Opt. Lett. **16**, 1732–1734 (1991). [CrossRef] [PubMed]

**10. **Y. Quiquempois, P. Niay, M. Douay, and B. Poumellec, “Advances in poling and permanently induced phenomena in silica-based glasses,” Current Opinion in Solid State & Materials Science **7**, 89–95 (2003) [CrossRef]

**11. **A. Ozcan, M. J. F. Digonnet, and G. S. Kino, “Cylinder-assisted Maker-fringe technique,” Electron. Lett. **39**, 1834–1836 (2003). [CrossRef]

**12. **A. Ozcan, M. J. F. Digonnet, G. S. Kino, and Edward L. Ginzton Laboratory, Stanford University, Stanford, CA 94305, are preparing a manuscript to be called “Detailed analysis of inverse Fourier transform techniques to uniquely infer second-order nonlinearity profile of thin films.”

**13. **M. Mukherjee, R. A. Myers, and S. R. J. Brueck, “Dynamics of second-harmonic generation in fused silica,” J. Opt. Soc. Am. B **11**, 665–669 (1994). [CrossRef]

**14. **T. G. Alley, S. R. J. Brueck, and R. A. Myers, “Space charge dynamics in thermally poled fused silica,” J. Non-Cryst. Solids. **242**, 165–176 (1998). [CrossRef]

**15. **D. Faccio, V. Pruneri, and P. G. Kazansky, “Dynamics of the second-order nonlinearity in thermally poled silica glass,” Appl. Phys. Lett. **79**, 2687–2689 (2001) [CrossRef]