♯SHAARP.si is an open-source package for deriving and simulating reflected optical second harmonic generation (SHG) from a single interface (si), typically the surface of a single crystal. Optical SHG describes the process where two photons of frequency $\omega$ interact with a nonlinear medium (a crystal) to create a photon at $2\omega$ , so called the SHG process.  (A followup package for multiple interfaces is currently being developed).

This package builds in the most general approach to both analytically and numerically solving the surface SHG response of a single crystal surface with arbitrary crystal symmetry, arbitrary orientation and a complex dielectric function (complex refractive indices).

As a very brief primer, the SHG interaction is given by $P_i^{2\omega} = d_{ijk}E_j^{\omega}E_k^{\omega}$ , where E are the fundamental electric fields of photons at $\omega$ frequency, P is the nonlinear polarization at $2\omega$ frequency created inside the crystal, and $d_{ijk}$ is a third rank polar tensor describing the nonlinear optical property of the crystal. The subscripts i, j, k are dummy subscripts denoting the polarization directions of the respective quantities; each of these indices can be 1, 2, or 3, that represent the orthogonal crystal physics axes denoted in ♯SHAARP as $(Z_1,Z_2,Z_3)$. The SHG tensor, $d_{ijk}$, thus has 3x3x3=27 terms; however crystal symmetry can significantly reduce the number of non-zero terms in the tensor. As described in the manual in detail, besides the crystal physics axes, there are four other sets of axes we will use in this package, namely, the crystallographic axes $(a,b,c)$, the lab axes, $(L_1,L_2,L_3)$, principal axes, $(Z_1^{princial},Z_2^{principal},Z_3^{principal})$ and the polarization axes (s, p); relationships between them is important to understand in order to usefully employ this package in experiments.

One critical application of the code is to provide analytical expressions to provide insight of intrinsic properties through the experimental observation. By such fitting, one can obtain structural information, polarization direction, as well as nonlinear optical susceptibilities. A second application is to simulate SHG response under various polarization states of incident photons, orientations and geometric considerations.

We request that you cite the following technical reference in any work for which you used ♯SHAARP:

Cite this work:

• Zu, R., Wang, B., He, J. et al. Analytical and numerical modeling of optical second harmonic generation in anisotropic crystals using ♯SHAARP package. npj Comput Mater 8, 246 (2022). https://doi.org/10.1038/s41524-022-00930-4

• GitHub
• Documentation
• Arxiv
• The ♯SHAARP package for multi-layers analysis (♯SHAARP.ml):
  • Code repository: https://github.com/bzw133/SHAARP.ml
  • Documentation of ♯SHAARP.ml: https://shaarpml.readthedocs.io/en/latest/
  • Tutorial videos ♯SHAARP.ml: https://youtu.be/YiKRm6i7kNs

♯SHAARP.si is written as a notebook using Wolfram Language and need to run with Mathematica®

• Download SHAARP_V#.zip in the file list to start using ♯SHAARP.

• To install Mathematica, please refer to the Installing Mathematica
• If you already have a Wolfram account, you can log in through the portal and download the software.

• Wolfram Player is a free software offered by Mathematica® to interact with Mathematica® Notebook.
• The download page can be accessed here
• Note: Wolfram Player provides access to most of the capabilities of ♯SHAARP.si except for Full Analytical expressions.

1. Unzip the file and open the SHAARP_si.nb.
   • Note: It is recommended to keep all the files in the same directory to access the full features of ♯SHAARP.si
   • Make sure that the Dynamic Evaluation has been enbled (it is enabled by default).
2. From the menu Evaluate -> Evaluate Notebook
   • Note: This process will clear out all the definitions from other notebooks and enable the "Notation" package for the analytical solutions.
3. After ~30s waiting time for initialization, you will see the main panel:
4. The main panel contains four parts
   • Input panels specify all the input parameters
   • Output panels give the output diagrams and equations
   • Progress bar show the progress after clicking the Update button
   • Update button execute the program by clicking it

• On the left-hand panel (LHP), scroll down to “Case Studies”.
• Click on LiNbO₃ (1120) MTI X-cut. All the parameters will be autofilled for this case.
• In the “Functionality”, select “SHG Simulation”. Leave the other LHP settings in their default settings.
• At the top R.H. corner, click “Update". (Please be patient, it may take a few seconds; Watch the progress bar on the top to track when the calculation is done.)
• The second row of the right-hand panel (RHP) should depict the SHG polar plots, $I^{2\omega}(\varphi,\psi)$ and $I^{2\omega}(\varphi,\psi+\pi/2)$ for the default settings of polarization and plane of incidence (PoI) in the LHP.
• Placing the cursor on the polar plots,and you can view a detailed description of the plots.

Now you have obtained your first numerical simulation of SHG polarimetry from a single interface. Next, let’s explore various polarization settings for the fundamental $\omega$ and SHG $2\omega$ waves.

• The definition of the incident polarization at frequency $\omega$ is given as, $E=(E_p,E_s)=E_0(\cos\varphi, \sin\psi e^{i\Delta\delta})$, where $E_p$ and $E_s$ are the p and s polarization components that are, respectively, parallel to and perpendicular to the PoI. The PoI is a plane formed by the incident wavevector, $\pmb{k}^{\omega}$, and the normal to the crystal surface. The polarization of the SHG wave (at frequency $2\omega$ is measured using a linear polarizer at an angle $\psi$, where $\psi=0^o$ and $90^o$ corresponds to the p-polarized and s-polarized SHG light, respectively.
• In the “Polarimetry Settings”, you can vary
  1. the incident angle ( $\theta^i$, in the range between $0-90^o$ ),
  2. the incident ( $\omega$ ) polarization direction ( $\varphi$ )
  3. output ( $2\omega$ ) polarization direction ( $\psi$ ), and
  4. the ellipticity of incident wave ( $\Delta\delta$ )
• In the “Incident Angle θⁱ(^o)”, set incident angle (at frequency $\omega$)  to $0^o$, and press “Update” to evaluate the changes. In the first row on the RHP, the polarized SHG response will change. The change of incident angle will also be revealed in the “Probing Geometry” plot in the second row of the RHP.
• Now set the “Output SHG Polarization” to “Fix Analyzer”. A “Fixed Analyzer Angle” option will appear, allowing you to select a specific output polarization direction. Keep the angle at $0^o$ for now (corresponding to p-polarization), and click “Update”. The SHG polar plots will change accordingly, and you can view the change to the polarization setting in the Polarization Relations plot in the second row of the RHP.
• Now click "3D Schematic" and then “Update”. The 3D Schematic option will allow you to visualize the Probing Geometry and Polarization Relations in 3D, and you can change the view direction by dragging the plots.
• You can also change the incident polarization. Try playing with the polarization settings and press “Update” to generate new plots.

Now you have some experience in setting up the polarization conditions for both fundamental and the SHG waves, in viewing your polarization settings using Probing Geometry and Polarization Relations plots, and generating different SHG polarimetry plots. Next, we will play with the crystal settings.

• In the “Crystal Structure Section”, the user has to specify the point group symmetry and the lattice parameters in order to properly determine the nonvanishing SHG coefficients as well the relevant crystal plane orientations. As an example, LiNbO₃ has the point group 3m and its lattice parameters have been well studied. Presently, leave the default settings as they are from selecting the Case Studies “LiNbO₃ (1120) MTI X-cut”.
• “Crystal Orientations” determines how the single crystal is oriented for the study with respect to the lab coordinates. It requires two user inputs:
  1. The orientation of the crystal plane of surface
  2. The orientation of the plane of incidence(PoI)
• These two pieces of information can be input in one of two ways: (a) By using Miller indices, (hkl), by clicking “Use Miller Indices (hkl)”. (b) By specifying the crystal physics axes with respect to the lab coordinates by clicking “Use crystal physics direction.” Note that $L_2$ lab axis is always fixed be perpendicular to the PoI.
• For the “LiNbO₃ (1120) MTI X-cut” case, the default orientation under “Use Miller Indices (hkl)” is set such that $(110)$ is the surface plane, and $[1 \overline{1} 0]$ is perpendicular to the PoI. Hovering your mouse over ?(hkl) describes the relation between 4-index and 3-index Miller indices for the trigonal and hexagonal systems.
• The “Use crystal physics direction” should be automatically populated for the preset “LiNbO₃ (1120) MTI X-cut”.  Otherwise, by hovering over the ?crystal physics specifies the relationship between the crystal physics and crystallographic axes for the selected point group. You can use this information to directly evaluate and input this information as well.
• Select 3D Schematic in the “Polarimetry Settings and press “Update”. You can then view the crystal orientations $(Z_1,Z_2,Z_3)$ in relation to the lab frame $(L_1,L_2,L_3)$ in the 3D Probing Geometry plot. You can click, hold and rotate the 3D schematics to get a 3D perspective.
• Note that for the “LiNbO₃ (1120) MTI X-cut” case, $Z_3$ is parallel to $L_1$ in 3D Probing Geometry plot.
• Next, change the h,k, and l under the “[hkl] -> Direction Perpendicular to the Plane of Incidence” to 0, 0, and 1. Press “Update” to update the calculation. The change means that you have rotated your crystal 90 degrees in-plane.
• Note that $Z_3$ is now parallel to $L_2$ in the 3D Probing Geometry plot.

Now you have finished the tutorial on the orientations. Detailed discussions about various coordinate systems and orientations can be found in the manual. Question marks next to “Use Miller Indices (hkl)” can help you quickly refer to the definition of $(Z_1,Z_2,Z_3)$ based on the point group you have selected.

In this section, you will need to provide either complex refractive index $\widetilde{n}$ or the complex relative dielectric permittivity $\widetilde\varepsilon$ for both the $\omega$ and $2\omega$ frequencies. Note, these tensors are property tensors specificed in the crystal physics coordinates and do not need to be changed when changing the crystal orientation.

• You can change the complex refractive index at $2\omega$ to complex values. You can change $(2.2529, 2.2529, 2.1604)$ to $(2.2529+I, 2.2529+I, 2.1604+2I)$. Press “Update” to update the calculation.
• Note, use capital $I$ to represent the imaginary part in Mathematica.
• Now the polar plots are updated with a real refractive index tensor at $\omega$ but a complex refractive tensor at $2\omega$. Next, try making the dielectric tensor at $\omega$ to be complex as well.

• Based on the point group symmetry you have selected in the “Crystal Structure”, the nonvanishing terms in the SHG tensor can be uniquely determined and provided in the Voigt notation¹.
• The default setting for the “LiNbO₃ (1120) MTI X-cut”provides the nonlinear coefficients of LiNbO₃ measured using the fundamental wavelength at 800nm. You can try to change those values and press “Update” to evaluate the influence on the SHG polar plots.

Now you have finished a quick tutorial on the input panels and performing the SHG simulation. ♯SHAARP not only provides simulations of the polarized SHG response but also features in generating analytical expressions for fitting experimentally polar plots for determining nonlinear SHG coefficients and the crystal symmetry.

Partial analytical expressions generate analytical expressions with only SHG coefficients as the unknown variables, providing a reliable way to experimentally determine SHG coefficients by fitting these expressions to the experimental polar plots. In this function, the numerical values of linear optical tensors, orientation, and the incident angle will be provided by the user and assumed to be known while the unknown SHG coefficients are left as variables. On the other hand, the Full analytical expressions provides the complete variables-based analytical expressions where the linear and nonlinear optical properties, as well as the various polarization and incidence angles of measurement are assumed to be variables. This method will provide a comprehensive expression for the wave mixing in the nonlinear medium.

• Taking LiNbO₃ as an example. Click LiNbO₃ (1120) MTI X-cut in the Case Study, in order to to use the default parameters.
• Click Partial Analytical Expressions in the Functionality section. Press “Update” to initiate the calculation.
• The RHP will display the final intensity expressions of $I^{2\omega}(\varphi,\psi)$, and $I^{2\omega}(\varphi,\psi+\pi/2)$.
• The Copy button to the left of the expressions will allow you directly copy the expression and paste it into another document such as a Mathematica notebook, or Microsoft files.

• Taking LiNbO₃ as an example. Click LiNbO₃ (1120) MTI X-cut in the Case Study, in order to to use the default parameters.
• Click Full Analytical Expressions in the Functionality section. Press “Update” to initiate the calculation. Note, this calculation could take up to a few minutes due to the amount of calculations involved. Please be patient! The progress bar shown at the top of the GUI helps keep track of the progress in the calculation.
• The RHP will display the equations for reflected electric fields at $2\omega$ frequency, $E^{2\omega}(\varphi,\psi)$ and $E^{2\omega}(\varphi,\psi+pi/2)$, and provide the complete set of equations in a separate notebook file, Full Analytical Expressions.nb,in the same directory of ♯SHAARP.si.
• The Copy button to the left of the expressions will allow you directly copy the expression and paste it into another document such as a Mathematica notebook, or Microsoft files.

Note, in the RHP, you can still use Copy to copy equations though the equations are not displayed.

• The analytical expressions are provided in sequence. For example, $E^{2\omega}(\varphi,\psi)$ is expressed using C^T,ee,2ω and E^T,e,2ω. The expressions for the C^T,ee,2ω and E^T,e,2ω will also be provided in the output expressions as a function of the material property coefficients and the measurement geometry.

Materials Properties Preset Values provides a way to store your input information so that you can review those input later. It works similarly to the buttons in the Case Study but allows users to customize input settings based on the need. This setting is particularly useful if you are varying some input settings to explore the changes towards the SHG response.

• Take LiNbO₃ as an example. Click LiNbO₃ (1120) MTI X-cut in the Case Study, to use the default parameters. Click SHG Simulations in the Functionality section. Press “Update” to initiate the calculation.
• Click “Preset 1” to store the input information including orientations, crystal structure, linear optical tensor values, and SHG tensor values. You can set labels for your presets and enter “LNO” in the “Preset 1 Label”. Use “Update” to secure the settings. Hovering the cursor on the pset buttons will provide more information about preset values.
• Now, you can click GaAs(111) in the Case Study and press “Update” to evaluate different materials. Go to “Preset 2” and use “Update” to make changes to the settings.
• By pressing “Preset 1”, the input information will direct you back to the settings for preset 1. You can also click “Preset 2” to return to saved settings for preset 2.
• If you click the “Clear Presets”, this process will erase all the saved settings for all four of the presets. You can redefine presets after clearing the definition.

• Detailed description of the method can be found here
• Some specific cases can be found here

1. Denev, S. A., Lummen, T. T. A., Barnes, E., Kumar, A. & Gopalan, V. Probing Ferroelectrics Using Optical Second Harmonic Generation. Journal of the American Ceramic Society 94, 2699–2727 (2011).
2. Boyd, R. W. & Prato, D. Nonlinear Optics. (Academic Press, 2008).
3. Bloembergen, N. & Pershan, P. S. Light Waves at the Boundary of Nonlinear Media. Phys. Rev. 128, 606–622 (1962).
4. Miller, R. C. & Nordland, W. A. Absolute Signs of Second-Harmonic Generation Coefficients of Piezoelectric Crystals. Phys. Rev. B 2, 4896–4902 (1970).
5. Jerphagnon, J. & Kurtz, S. K. Maker Fringes: A Detailed Comparison of Theory and Experiment for Isotropic and Uniaxial Crystals. Journal of Applied Physics 41, 1667–1681 (1970).
6. Herman, W. N. & Hayden, L. M. Maker fringes revisited: second-harmonic generation from birefringent or absorbing materials. J. Opt. Soc. Am. B, JOSAB 12, 416–427 (1995).
7. Shoji, I., Kondo, T., Kitamoto, A., Shirane, M. & Ito, R. Absolute scale of second-order nonlinear-optical coefficients. J. Opt. Soc. Am. B, JOSAB 14, 2268–2294 (1997).