Optical Second Harmonic Generation in Anisotropic Multilayers with Complete Multireflection Analysis of Linear and Nonlinear Waves using ♯ SHAARP. ml Package

Optical second harmonic generation (SHG) is a nonlinear optical effect widely used for nonlinear optical microscopy and laser frequency conversion. Closed-form analytical solution of the nonlinear optical responses is essential for evaluating the optical responses of new materials whose optical properties are unknown a priori . A recent open-source code, ♯ SHAARP. si , can provide such closed form solutions for crystals with arbitrary symmetries, orientations, and anisotropic properties at a single interface. However, optical components are often in the form of slabs, thin films on substrates, and multilayer heterostructures with multiple reflections of both the fundamental and up to ten different SHG waves at each interface, adding significant complexity. Many approximations have therefore been employed in the existing analytical approaches, such as slowly varying approximation, weak reflection of the nonlinear polarization, transparent medium, high crystallographic symmetry, Kleinman symmetry, easy crystal orientation along a high-symmetry direction, phase matching conditions and negligible interference among nonlinear waves, which may lead to large errors in the reported material properties. To avoid these approximations, we have developed an open-source package named Second Harmonic Analysis of Anisotropic Rotational Polarimetry in Multilayers ( ♯ SHAARP. ml ). The reliability and accuracy are established by experimentally benchmarking with both the SHG polarimetry and Maker fringes predicted from the package using standard materials. The ♯ SHAARP. ml can be accessed through GitHub (https://github.com/bzw133/SHAARP.ml).


Introduction
The development of coherent laser light over a broad frequency spectrum from nearinfrared and visible to terahertz (THz), ultraviolet, and X-rays regimes [1][2][3][4] has driven much of science and technology in the past decades, ranging from sensing, communications, biomedical instruments, imaging, and most recently nuclear fusion research. [5][6][7][8][9][10] Since the discovery of lasers in 1960 and the nonlinear optical effect in 1961 11,12 , nonlinear optics has been a primary source for generating a continuously tunable electromagnetic spectrum. In the last two decades, quantum communications and computing have relied on using nonlinear optics to generate entangled photons and to achieve ultrafast all-optical switching. [13][14][15] Optical second harmonic generation (SHG) refers to the nonlinear optical process where two photons of the same energy (ℏ ) combine to generate a new photon of higher energy (2ℏ ) in a nonlinear optical (NLO) medium. This phenomenon is described by the nonlinear polarization, 2 = (2) , generated in the NLO material at 2 frequency by the electric field of the incident light, at frequency . 16 Here, (2) is the second-order nonlinear optical susceptibility represented by a third-rank tensor (with 18 independent components). If the × × × N/A √ a R and T refer to reflection and transmission, respectively. b SI represents single interface. Numbers reflect the number of layers. c High symmetry means samples are oriented along a high-symmetry direction. d p-or s-refer to the electric fields of electromagnetic waves either parallel or perpendicular to the plane of incidence, respectively. e MR represents multiple reflections of waves, e&o represents homogeneous waves at their corresponding frequency, or 2 (e for extraordinary and o for ordinary waves), and 2 stands for nonlinear polarization that gives rise to SHG effects. f N/A refers to not applicable.

Results and Discussion
Theoretical background Figure 1a presents the ray diagram of linear and nonlinear waves through a multilayer system adopted in ♯SHAARP.ml. Without loss of generality, we assume the first layer (M1) to be SHG active. In a more general case, all layers can (but need not) be SHG active in experiments. When a monochromatic plane wave at frequency is incident upon the system, the electromagnetic properties of the plane wave inside the system are governed by the wave equation at frequency, × × + � � 1 1 � 1 2 � 1 3 � 2 1 � 2 2 � 2 3 � 3 1 � 3 2 � 3 3 where , � and are respectively the electric field inside the medium at frequency, anisotropic dielectric tensor components in the lab coordinate system (LCS), and magnetic permeability tensor at frequency. The will be assumed to be vacuum permeability for a nonmagnetic system, ~0 , where is the identity matrix. The subscripts and are dummy indices describing the direction of each tensor component of the anisotropic dielectric susceptibility tensor in the LCS, denoted as � LCS . Note that � LCS can be complex to account for absorption. Four coordinate systems are utilized, namely, principal coordinate system (PCS), crystal physics coordinate system (ZCS), crystallographic coordinate system (CCS), and lab coordinate system (LCS). In PCS, the complex dielectric susceptibility tensor is diagonalized.
ZCS is the orthogonal coordinate system in which the property tensors are defined, such as dielectric susceptibility tensor, SHG tensor, piezoelectricity tensor, etc. 45 The CCS describes the coordinate system formed by the basis vectors of the unit cell (which are not necessarily orthogonal), and LCS is an orthogonal coordinate system of the model system with the plane of incidence (PoI) coincides with the L1-L3 plane as shown in Figure 1a. Note that PCS, ZCS, and LCS are orthogonal coordinate systems, while the CCS can be non-orthogonal depending on the crystal symmetry. Equation (1) is a generalized eigenvalue problem that can be solved routinely. 46 The resulting eigenvalues and eigenvectors are related to the effective refractive indices and electric field directions for both ordinary and extraordinary waves. Due to reflectance at various interfaces, both forward and backward propagating waves exist in the heterostructure.
The resulting backward propagating wavevectors can be described as The superscripts e , o , F, and B , respectively, represent extraordinary, ordinary, forwardpropagating, and backward-propagating waves. M represents the th medium in the heterostructure. Similarly, the full electromagnetic properties of backward propagating waves can be obtained using Equations (1) and ( The optical dipolar second harmonic generation is defined by the generation of nonlinear polarization at 2 frequency when the NLO materials are pumped by the incident electric fields at frequency. The nonlinear polarization is defined as where M 2 , M , 0 , (2) , and are nonlinear polarization, fundamental electric field, vacuum dielectric permittivity, second-order nonlinear optical susceptibility, wave vector of the source wave, and position vector, respectively. Since arbitrary layers can be SHG active, M 2 will appear when the th layer is SHG active, as denoted by the subscript M . The generated nonlinear polarization is often known as the source wave that gives rise to the nonlinear optical effects. It is important to note that during the propagation of fundamental fields, the nonlinear polarization is generated throughout the entire optical path of M , according to equation ( The propagation of 2 is confined to the propagation of the fundamental wave at that generates it, and the corresponding 2 is hence called the bound wave or inhomogeneous wave.
On the other hand, the SHG wave generated by the bound wave can freely propagate governed by the direction specified by Snell's law at 2 , hence it is called the free wave or the homogeneous wave.
The anisotropic three-wave mixing phenomena is revealed in equation (10), where material anisotropy is taken into account. In each SHG active medium (M ), the forward and backward nonlinear wavevectors can thus be identified as S,2 = 2 eF, , 2 oF, , eF, + oF, , 2 eB, , 2 oB, , eB, + oB, , eF, + eB, , eF, + oB, , oF, + eB, , and oF, + oB, . The wavevectors for the ten nonlinear polarizations in the th layer are thus denoted as ( eFeF,2 , oFoF,2 , eFoF,2 , eBeB,2 , oBoB,2 , eBoB,2 , eFeB,2 , eFoB,2 , oFeB,2 , oFoB,2 )Mi for clarity, as shown in Figure 1a. For example, a nonlinear polarization eFoB,2 is formed when a forward propagating extraordinary wave ( eF, ) and a backward propagating ordinary wave ( oB, ) are combined. However, the wave mixing terms containing both forward and backward waves, such as eFeB,2 and eFeB,2 , are often dropped or ignored in existing literature due to a large phase mismatch. 29,30 Although these terms form standing waves propagating parallel to the layers, the standing waves at both the top and bottom surfaces of each layer can still contribute to the boundary conditions. For example, a nonlinear polarization ( eFeB,2 ) can be generated by a mixture of eF, and eB, at top or bottom surfaces leading to additional components in the boundary conditions. Therefore, we have implemented the mixing term in ♯SHAARP.ml, resulting in, at most, ten distinct nonlinear polarizations of different combinations of wavevectors for each SHG active layer. These ten waves are shown as ten different arrows in Fig. 1a.
The particular solutions of equation (11) can be obtained using the method described in previous work. 21 For example, the electric field of the nonlinear polarization induced by the mixture of two forward extraordinary waves can be written as eFeF,2 = eFeF,2ω ( eFeF,2 • −2 ) , where eFeF,2ω is a vector describing the direction and magnitude of the resulting bounded electric field due to the nonlinear polarization. Thus, all electric and magnetic fields generated by the ten distinct nonlinear polarizations can be uniquely identified by solving equation (11). On the other hand, the general solution of equation (11), which represents the homogeneous waves, can be calculated following the same procedure as solving equation (1) but at 2 frequency. Four nonlinear waves will be obtained to fully describe the multiple reflections of homogeneous waves, namely, ( eF,2 , oF,2 , eB,2 , oB,2 ) M , whose field strengths are determined using the boundary conditions to be described below.
The momentum conservation and energy conservation of the generated 2 waves lead to the following boundary condition: where is the phase difference for a forward wave propagating from top to bottom surface and

Outline of ♯SHAARP
The theoretical method described in the preceding section is implemented using Wolfram Mathematica with a user-friendly GUI and a detailed tutorial, which can be found in Ref. 50 .
Following the naming convention of our previous work, we named the newly developed software capable of modeling optical SHG of multilayer system as ♯SHAARP.ml. Figure 2 illustrates the calculation procedure of ♯SHAARP.ml. First, with a given point group symmetry, the dielectric tensor in the ZCS, and its orientation relative to the LCS coordinate system as inputs, one can conveniently obtain the mutual relations among the four coordinate systems within ♯SHAARP.ml, and thus define the geometry of the system. Then, by solving the wave equation  Case studies using ♯SHAARP.ml In the following, we present our experimental measurements of the SHG responses for a few typical nonlinear optical crystals and their heterostructures to demonstrate how they can be interpreted by numerical and semi-analytical analyses using ♯SHAARP.ml. In particular, we studied the Maker fringes of pure and Au-coated quartz single crystals and the SHG polarimetry of LiNbO3, KTP, and ZnO//Pt//Al2O3 heterostructure. We performed a predictive modeling of a bilayer consisting of two SHG active materials, namely, X-cut LiNbO3 on Z-cut quartz, which can be helpful in distinguishing the ferroelectric domain states of LiNbO3 from the SHG intensity map. These examples not only serve as benchmark tests of ♯SHAARP.ml against known NLO materials covering a wide range of types (uniaxial, biaxial, and absorbing) but also demonstrate the broad applicability of ♯SHAARP.ml to a variety of situations (e.g., Maker fringes, polarimetry, quantifying the effect of adopting different assumptions in the SHG modeling, analytical fitting to extract absolute values of SHG coefficients, and predictive simulations of SHG responses of NLO heterostructures).

Maker fringes of α-quartz single crystal
The study of α-quartz in nonlinear optics can be traced back to the discovery of second harmonic generation in 1961. 11 The first benchmark study for ♯SHAARP.ml is performed using the single crystalline α-quartz, which has been extensively investigated previously using the Maker fringes method. 22,24,28,30 The SHG coefficient 11 has been measured to be 0.3 pm/V. 51 In this case study, we demonstrate the capability of ♯SHAARP.ml in obtaining the semi-analytical expression for Maker fringe response and benchmark analysis with both existing models in the literature 24,28 and our experimental investigations. Figure 3 shows the comparison among numerical simulation results from ♯SHAARP.ml with various modeling conditions and existing results using analytical methods. 24,28 The Maker fringes condition is summarized in Figure 3a.
The fundamental wavelength ( ) is 1064 nm and the generated SHG signal from a 300 µm Xcut quartz is analyzed. Both the fundamental and SHG waves are p-polarized. Two widely applied Maker fringes models are utilized for comparison, namely the JK (Jerphagnon & Kurtz 24 ) method and HH (Herman & Hayden 28 ) method. The JK method was developed for an isotropic medium with an assumption that only forward propagating waves are involved. 24 The HH method extended this model to a birefringent uniaxial system with multiple reflections of homogeneous waves (free waves) at 2 frequency, but not for the inhomogeneous waves or linear waves. ♯SHAARP.ml involves multiple reflections for both linear and nonlinear waves (homogeneous and inhomogeneous) and thus can be reduced to JK or HH methods by making the corresponding assumptions. Schematics of the assumptions made for the three approaches can be found in Supplementary Note 1, Figure S1. Figures 3b and 3c illustrate the three Maker fringes patterns obtained from the HH method (denoted as analytic HH) and numerical analysis using ♯SHAARP.ml with both JK and HH modeling conditions, denoted as ♯SHAARP(JK) and ♯SHAARP(HH). 24,28 The blue dots, yellow and green lines correspond to analytic HH, ♯SHAARP(JK) and ♯SHAARP(HH), respectively. All three Maker fringe patterns are consistent with the literature. 28 In particular, analytic HH and ♯SHAARP(HH) show good agreement, demonstrating ♯SHAARP.ml can accurately reproduce the prior results. Figure 3c shows the magnified area of the dashed box region in Figure 3b. By enabling the multiple reflections of homogeneous waves at 2 frequency, ♯SHAARP(HH) produce additional fine fringes at from 20° to 30°, which are absent for ♯SHAARP(JK). This difference indicates that the interference between forward and backward homogenous 2 waves results in these fine fringes. In contrast to the fine fringes originating from the interference of the fundamental waves, the broader envelope in the SHG intensity with respect to (interval ranging across tens of degrees visible in Figures 3b, 4b and 4d) carry the essential information associated with the interference between the homogeneous and inhomogeneous waves. This interference originates from the phase difference between the source waves ( S,2 ) and the homogeneous waves ( e,2 and o,2 ) accumulated throughout the bilayer structure, and thus, the broader envelope is extremely sensitive to the changes in the crystal thickness and refractive indices at both and 2 frequencies. Therefore, SHG Maker fringes can be utilized as a sensitive probe of wafer uniformity. 53 For example, with a thickness variation of 1 , the Maker fringes change drastically, as demonstrated in Supplementary Note 4 (see Figure S4). It is worth noting that the crystal thicknesses determined in Figures 4b and 4d are slightly different, i.e., 123.6 and 121.2 , respectively, due to the change of probing positions and nonuniform thickness across the sample (10 μm variation across a 10 mm × 10 mm sample), as confirmed by the stylus profilometry. In addition, we note that the example presented in Figure 4c and d also illustrates the capability of ♯SHAARP.ml in handling multiple layers with strong reflections.
The phase difference between two propagating waves is critical to determining their interference, e.g., being constructive and destructive for in-phase and out-of-phase situations respectively. With ♯SHAARP.ml, we show that different ways to compute the relative phase  Figure S5). Such discrepancy may come from the fact that a small beam size comparative to the crystal thickness is used in the experiment, where a sizeable beam overlap and finite resolution of angles are essential for the interference to become observable in the experiments.
Therefore, for the quartz case, taking only the vertical phase along 3 direction will be sufficient in the SHG analysis throughout the current work.

LiNbO 3 and KTP Single Crystals
LiNbO3 and KTiOPO4 (potassium titanyl phosphate, KTP) have been widely studied for decades owing to their excellent nonlinear optical properties. [54][55][56] Their well-established nonlinear optical susceptibilities make the two crystals suitable for benchmarking analysis.
Utilizing the partial analytical expressions generated by ♯SHAARP.ml, the experimental polarimetry results can be analyzed to extract relative ratios of SHG coefficients, and the absolute SHG coefficients of the two single crystals can be obtained using α-quartz as the reference.
LiNbO3 crystallizes in a trigonal structure with the point group 3m and has a bandgap of around 3.8 eV. 55 in-plane isotropy in this orientation. As the crystal is tilted towards a larger incidence angle, the projection of 33 to the 1 increases, leading to an increase in the p-polarized SHG intensity, as seen in Figure 5c. By fitting two LiNbO3 crystals with different orientations and using quartz as the reference, the extracted ratios and absolute values of the SHG coefficients of LiNbO3 are summarized in Table 2, which agree well with previously reported values 24,29 .
KTP adopts an orthorhombic crystal structure with a point group of mm2. It is classified as a biaxial material with distinct optical responses along all three crystal physics axes. Thus, a careful analysis of full anisotropy and the presence of two optical axes are critical in optical modeling. In this study, we used two KTP slabs simultaneously, namely ~370 μm X-cut ( (100) orientation) slab and ~570 μm Y-cut ((010) orientation) slab, to analyze the full SHG tensor.
Both c axes are placed along the 2 direction, and their two optical axes lie in 1 -3 plane (see the experimental orientations in Supplementary Note 6, Figure S6). 45,58 Figures 5e and 5f are the SHG polar plots for p-and s-polarized SHG response, respectively. Four incident angles are utilized to identify five unknown SHG susceptibilities uniquely ( i = 0°, 10°, 20°,and 40°).
Using partial analytical expressions generated by ♯SHAARP.ml, the SHG polarimetry fittings show good agreement between the theory and experimental data, and the extracted ratios and absolute values of SHG coefficients of KTP are summarized in Table 2. As discussed in previous work, 21 the symmetry assumptions, such as that of isotropy, can lead to errors of up to 30% in the ratios between SHG coefficients, depending on the anisotropy of the materials. In this work, our discussion will focus on the influence of Kleinman's symmetry Kleinmann symmetry assumed. Comparing the four cases, we found most of the obtained SHG ratios vary within 20-30%, which are commonly comparable to the error bars. The NMR case is close to the ♯SHAARP.ml case, implying that the HH method may be a good approximation for studying KTP with photon energies below its bandgap. The KS case, however, can introduce relatively large deviations in the obtained coefficient ratios such as a 60% error for 31 / 32 in KTP.

ZnO//Pt//Al 2 O 3 thin films
ZnO has been widely studied for decades for electronics, photonics, and optoelectronics applications owing to its large piezoelectric coefficients, large exciton binding energies, wide optical bandgap, and good chemical and thermal stability. [64][65][66] Recently, ZnO with Mg substitution (Zn1-xMgxO) has been shown to possess ferroelectricity, paving its way toward waveguides and quasi-phase-matched (QPM) frequency conversion devices. 15,43 Though the nonlinear optical process in ZnO has been extensively explored in both bulk and thin films forms, its nonlinear optical susceptibilities have been reported with a large scatter in the values from less than one pm/V to hundreds of pm/V, indicating either sample variations or inconsistent modeling of the SHG data. [67][68][69][70] In this work, we select 159 nm ZnO//200nm Pt//0.5mm Al2O3 as an example to demonstrate the capabilities of ♯SHAARP.ml in probing thin films on substrates with a bottom electrode and the importance of multiple reflections in the analysis.
As described in earlier work, ZnO was grown using RF magnetron sputtering and formed a stack of ZnO//Pt//Al2O3, as shown in Figure 6a. 43  Further reference against a wedged X-cut LiNbO3 yields the absolute SHG coefficients of the entire SHG tensor. Figure 6e summarizes the absolute SHG coefficients obtained from ♯SHAARP.ml in comparison with the cases under various assumptions (the meanings of the notations are consistent with the previous section). The "♯SHAARP.ml" case yields the absolute 33 = 6.6 ± 2.2 pm/V, which is close to early reported values for films and single crystals (~7.15 pm/V). 68 This indicates the film under study has good qualities and low optical loss.
Comparing the results from ♯SHAARP.ml with those from KS and No FB , the obtained absolute SHG coefficients are reasonably close. On the other hand, the multiple reflections play a more significant role in the analysis. As can be seen from the NMR case, the obtained nonlinear susceptibilities are greatly exaggerated by one order of magnitude. This is because the total SHG signals were attributed to the single propagation of nonlinear polarization from the top to the bottom surface instead of multiple bounces. To compensate for the path difference between NMR and FMR, the nonlinear susceptibilities have to be increased, leading to SHG of nearly 10 times higher than the actual value.

SHG active bilayers, LiNbO 3 //Quartz
The generated SHG signals, in general, contain both amplitude and phase information of materials, such as the direction of a static (zero frequency) spontaneous polarization, Ps, of ferroelectric materials. (Note that this static ferroelectric polarization is distinct from any optical polarization at optical frequencies we have discussed earlier). Two ferroelectric domains with antiparallel spontaneous polarizations (separated by a 180 • domain wall) will generate nonlinear optical polarizations with a phase shift, yet of the same amplitude. Thus, the corresponding SHG intensities are identical for the two domains, leaving the ferroelectric domain state indistinguishable based on the intensity alone. 16,33,75 The SHG interference contrast imaging has been developed to resolve this issue. [75][76][77][78][79] In this subsection, we employ ♯SHAARP.ml simulation to illustrate the basic idea of SHG interference contrast imaging, intimately (without an air gap in this example) placing a periodically poled X-cut LiNbO3 crystal on top of a Z-cut quartz crystal as a model system.
The principle of SHG interference contrast imaging is schematically shown in Figure 7a, where the blue and red are fundamental waves and SHG waves, respectively. 16,76 An additional quartz is placed beneath the LiNbO3 crystal (abbreviated as LNO) to generate the interference of the nonlinear waves through reflection. The nonlinear waves generated by LiNbO3 (denoted as L 2 ) and quartz (denoted as Q 2 ) will interfere to resolve the phase information of L 2 . Figure   7b shows  Looking forward, we expect that ♯SHAARP.ml can broadly streamline research in nonlinear optics. The complete and accurate analytical framework with editable assumptions from ♯SHAARP.ml can provide nonlinear optical solutions in an on-demand modality. As more integrated nonlinear optical devices and new topological superlattices are being developed, the capability of modeling these heterostructure can thus be an effective way to design, characterize, and optimize nonlinear optical response from complex systems. Furthermore, ♯SHAARP.ml provides a unique programmable platform for future extensions to new functionalities, such as other three-wave mixing processes, magnetic-dipole or quadrupole induced nonlinear optical effects, Gaussian beams with finite beam size, and inhomogeneous material systems.

Sample Preparation
Both α-quartz and LiNbO3 single crystals were obtained from MTI Corporation. The (112 � 0) and (0001) oriented LiNbO3, namely X-cut and Z-cut, were utilized in the analysis.
Since the definition of X-cut LiNbO3 from MTI is distinct from the orientations used in other analyses, 80,81 we have used the Miller indices for clarity. The X-cut and Y-cut KTP crystals were obtained from CASTECH Inc (Conex Systems Technology, Inc.). The ZnO//Pt//Al2O3 was prepared using RF magnetron sputtering, and the detailed growth procedure can be found in the earlier work. 43

Second-harmonic generation:
The second harmonic generation measurements were performed using a Ti: Sapphire femtosecond laser system with the central wavelength at 800 nm (1 kHz, 100 fs). The 1550 nm (1 kHz, 100 fs) was generated through an optical parametric amplifier, pumped by the 800 nm amplified laser. The SHG polarimetry measurements were performed using a combination of a

Data Availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Acknowledgment
The ♯SHAARP open-source package is supported by the US Department of Energy, Office of  Figure S1. The differences in the nonlinear optical modelings among JK, HH, and ♯SHAARP.ml methods. The dashed regions represent independent boundary conditions. Linear, Inhomo, and Homo represent linear waves at , inhomogeneous waves at 2 , and homogeneous waves at 2 , respectively.

Supplementary note 2:
The influences of averaging incident angle and thickness on the Maker fringes patterns are investigated, as shown in Figure S2. Experimentally, the lens with a focusing distance of 10 cm (f=10 cm) is used, and the fundamental beam has a diameter of 5 mm (d≈ 5 mm). Thus, the convergence angle is estimated to be −1 ( ) ≈ 3°. The spot size focused on the sample at = 800 nm is around 50 μm, and a near 10 μm thickness variation is observed across 10 mm × 10 mm sample. Thus, the thickness variation within the spot size is estimated to be 50 nm. Thus, averaged Maker fringes patterns with a 3° binning widow for and a 50 nm binning window for sample thickness (h) are investigated. The bandwidth of 800 nm laser was measured to have ±5 nm variation. The variation in the wavelength ( ) has a similar effect as the thickness variation, and both effects contribute to the phase accumulated throughout the crystal, i.e. = 2π ℎ. Thus, the variation in is estimated to result in ~0.6% variation in phase, equivalent to ~0.7 μm variation in ℎ for a 120 μm thick crystal.  Supplementary note 3: Figure S3. The complex refractive index of the Au layer determined using spectroscopic ellipsometry. The resulting layer thickness is determined to be 13.9 nm.
Supplementary note 4: Figure S4. The Maker fringes pattern of Quartz with different thicknesses. The dot corresponds to the experimental results shown in Figure 4b. The h is the crystal thickness used in the simulation. The thickness of the crystal near the probing area is determined to be 123.6 , while ±1 in h leads to a large change in the Maker fringes pattern.

Supplementary Note 5:
The phase calculation methods play a critical role in the predicted Maker fringes pattern. Figure   S5 compares two different phase calculation methods and explores the influence on the obtained