# Nonlinear Optics Through the Field Tensor Formalism

Julien Duboisset, Benoît Boulanger, Sophie Brasselet, P. Segonds, Joseph Zyss

► To cite this version:

Julien Duboisset, Benoît Boulanger, Sophie Brasselet, P. Segonds, Joseph Zyss. Nonlinear Optics Through the Field Tensor Formalism. Laser and Photonics Reviews, 2024, pp.2400411. 10.1002/lpor.202400411 . hal-04770878

HAL Id: hal-04770878

https://hal.science/hal-04770878v1

Submitted on 7 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear Optics Through the Field Tensor Formalism

Julien Duboisset,* Benoît Boulanger,* Sophie Brasselet, Patricia Segonds, and Joseph Zyss

The “field tensor” is the tensor product of the electric fields of the interacting waves during a sum- or difference-frequency generation nonlinear optical interaction. It is therefore a tensor describing light interacting with matter, the latter being characterized by the “electric susceptibility tensor.” The contracted product of these two tensors of equal rank gives the light-matter interaction energy, whether or not propagation occurs. This notion having been explicitly or implicitly present from the early pioneering studies in nonlinear optics, its practical use has led to original developments in many highly topical theoretical or experimental situations, at the microscopic as well macroscopic level throughout a variety of coherent or non-coherent processes. The aim of this review article is to rigorously explain the field tensor formalism in the context of tensor algebra and nonlinear optics in terms of a general time-space multi-convolutional development, using spherical tensors, with components expressed in the frame of a common basis set of irreducible tensors, or Cartesian tensors. A wide variety of media are considered, including biological tissues and their imaging, artificially engineered by various combinations of optical and static electric fields, with the two extremes of all-optical and purely electric poling, and also bulk single crystals.

## 1. Introduction

Optical processes, be they coherent or incoherent, linear or nonlinear, given or engineered, convolve material properties with that of interacting fields, both incoming and outgoing. They are related to classical or quantum processes, depending on the scale of interest, related to electron displacements in matter under external electromagnetic influences. They are known as polarization of matter phenomena in the general sense, whereby material and field properties are mutually intertwined. The polarizability is either generally given as such by nature or at the outcome of engineering processes leading to artificial structures. It is generally considered at the core molecular or bond scale and then upscaled to the wavelength scale or above. Relations between these two scales, including geometric and physical considerations, are a crucial part of the prediction and analysis of optical properties. Geometric considerations

entail framework rotations allowing to move from the molecular or bond scale to that of interest for electromagnetic waves, namely the principal dielectric axes for a crystal and more generally, any adequately chosen laboratory frame provided it is an orthonormal frame where external parameters can be controlled. As for the interacting fields themselves, they come to activate the built-in material properties so as to trigger various optical processes of interest. The fact that material and field properties cannot be dissociated is well illustrated by the wave-vector $\vec{k}^{\omega} = \omega n^{\omega}/c \vec{u}$ of a propagating electric field $\vec{E}^{\omega}$ at the angular frequency $\omega$, where $\vec{u}$ stands for the unit vector, $n^{\omega}/c$ for the inverse of the phase velocity of the corresponding wave, $c$ the speed of light in vacuum ad $n^{\omega}$ the refractive index in the material, itself related to its linear polarizability. Clearly, polarizability via the refractive index and propagation of light via the wave-vector are indissociable properties. Note that in addition, the refractive index is not only depending on the angular frequency but also on the polarization of the field, which is itself depending on the refractive index, thus adding another bind between field(s) and matter. This will play an important role especially in the crystalline part of this review.

The dual nature of optical processes has led, from the inception of modern studies at the end of the 19th century, to the mathematically correct definition of a (multi-)linear convolution response including successive steps of formalization due to the slow emergence of a full-fledged tensor algebra toolbox commonly available to the optical community. The output of such

J. Duboisset, S. Brasselet
Aix Marseille Univ
CNRS
Centrale Med
Institut Fresnel
Marseille 13397, France
E-mail: julien.duboisset@fresnel.fr

B. Boulanger, P. Segonds
Université Grenoble Alpes
CNRS
Grenoble INP
Institut Néel
Grenoble 38000, France
E-mail: benoit.boulanger@neel.cnrs.fr

J. Zyss
Laboratoire Lumière
Matière et Interfaces (LuMin)
ENS Paris-Saclay
CentraleSupélec
CNRS
Université Paris-Saclay
Gif-sur-Yvette 91190, France

The ORCID identification number(s) for the author(s) of this article can be found under https://doi.org/10.1002/lpor.202400411

© 2024 The Author(s). Laser & Photonics Reviews published by Wiley-VCH GmbH. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

DOI: 10.1002/lpor.202400411
a response takes the form of an induced polarization of matter vector that further radiates an outgoing electric field via Maxwell equations, whereas the input is a tensor comprising all possible products of the input electric field components for a given rank that defines the order of the process. The convolution response, also known as the (non-)linear impulse response, stands upstream to the (hyper-)polarization properties of the medium.

In this review, we focus on the description of electric fields interacting with matter. This has been recognized as a scalar product between the susceptibility tensor and the field tensor, which is the main concern of the present study. Whereas this notion cannot claim any conceptual originality as such, having been explicitly or implicitly present from the early pioneering studies in nonlinear optics, its practical use in specific, however important situations, has led to original developments. These have either considerably simplified or even in some cases short-circuited tedious derivations as per the earlier crude methods that may have obscured the underlying physics. It has also led to previously incomplete insights on optical phenomena by way of new tensor selection rules allowing us to classify, predict or forbid phenomena as well as to open the way to simultaneous optimization of both the material and the probing electric field(s). Seemingly haphazard results at the outcome of calculations generally performed in a Cartesian basis, can then be seen from a more global and rational perspective.

This review will be organized in view of two main objectives: the first one is aiming at depicting and matching the symmetry properties of both the electric field and electric susceptibility tensors, mainly in terms of spherical tensors, with components expressed in the frame of a common basis set of irreducible tensors abiding to canonical transformation rules. This part will be mainly applicable to situations whereby propagation considerations can be ignored, such as for surfaces or thin films of dimensions much smaller than the wavelength. Various scenarios of nonlinear interactions will be considered in this frame: harmonic generation, statistically oriented media, nonlinear nano-plasmonic where outward shapes become essential as opposed to internal symmetry in the more classical case of dielectrics, nonlinear spectroscopy including coherent nonlinear Raman effects, multiphoton transfer of angular momentum and the currently forefront domain of high order harmonic generation. Special attention will be given in this framework to statistically oriented media, be they given as such in biological tissues and their imaging, or somehow artificially engineered by various combinations of optical and static electric fields, with the two extremes of all-optical and purely electric poling. The terminology and relevance of multipolar field tensor and that of the matching multipolar electric susceptibility terms will be discussed in this perspective.

The second line explores situations where propagation considerations prevail, in connection with thick media and phase-matching considerations, with emphasis put on single crystals. Then the field tensor will be indexed by the nature of eigen-field polarizations, referring to an ordinary (0) or extraordinary (e) refractive index in the case of uniaxial crystals, to a larger (+) or smaller (-) refractive index in biaxial crystals, whichever is relevant toward the fulfilment of phase-matching conditions. In these crystals, the field tensor proves to be a unique tool: on the one hand, to identify crystal classes for which some phase-

matching polarization configurations are forbidden, which corresponds to orthogonality between the field tensor and the electric susceptibility tensor under consideration; on the other hand, to determine the relative signs and amplitudes of the coefficients of the electric susceptibility tensor to be characterized.

We will close this review at a downstream level of the (multiple) Fourier transform of the convolution response which is in general the starting point of theoretical as well as experimental studies. Details of the derivation of the Fourier response from the upstream general convolution are to be found in textbooks.[1,2] It is nevertheless important to keep in mind the initial (multi-)linear convolution that accounts for memory effects to be further reflected in dephasing factors in the Fourier response. These are of particular importance for resonant responses where the polarization response is actually retarded by the finite lifetime of the addressed excited levels of the material, be they electronic, vibrational or both, as it will be considered in this review.

# 2. Field Tensor Formalism

## 2.1. Interaction Energy

Introduced by Maker in the context of second harmonic scattering in liquids,[3] the field tensor derives originally from the constitutive relation between the response of a medium to a multiple electric fields excitation given by:

$$ \vec{P}(\omega_n) = \epsilon_0 \chi^{(n-1)} \vec{E}(\vec{k}_1, \omega_1) \dots \vec{E}(\vec{k}_{n-1}, \omega_{n-1}) \quad (1) $$

where $\vec{P}$ is the induced polarization vector, $\epsilon_0$ the vacuum permittivity, $\chi^{(n-1)}$ the $(n-1)$-order electric susceptibility tensor of rank $n$, and $\vec{E}$ the electric field.

From this constitutive relation, Maker used the definition of the interaction energy written as[3,4]

$$ W = -\epsilon_0 \chi^{(n-1)} \vec{E}^*(\vec{k}_n, \omega_n) \vec{E}(\vec{k}_1, \omega_1) \dots \vec{E}(\vec{k}_{n-1}, \omega_{n-1}) = -\vec{E}^* \cdot \vec{P} \quad (2) $$

$$ = -\epsilon_0 \chi^{(n-1)} \cdot F^{(n-1)} \quad (3) $$

where “.” stands for the contraction between two tensors of the same rank $n$, that is, $\chi^{(n-1)}$ and $F^{(n-1)}$, the Field tensor being defined as:

$$ F^{(n-1)} = [\vec{E}^*(\vec{k}_n, \omega_n) \otimes \vec{E}(\vec{k}_1, \omega_1) \otimes \dots \otimes \vec{E}(\vec{k}_{n-1}, \omega_{n-1})] \quad (4) $$

where $\otimes$ stands for the tensor product between the interacting electric fields, $*$ is the complex conjugated, $\omega$ is the angular frequency and $\vec{k}$ the wave vector.

Note that the rank of the field tensor corresponds to the $n$ electric fields involved in the light-matter interaction, that is, the electric fields exciting the matter plus the electric fields emitted by the polarization of matter. This formulation is only valid for coherent optical processes, that is, when the emitted electric fields obey a phase relationship with the incoming electric fields, excluding for example processes based on fluorescence or luminescence. The
field tensor depends exclusively on the linear optical properties of the materials, that is, the wavelength dispersion of the refractive index, the linear absorption, and also on its orientation with respect to the material. The field tensor can be divided into two parts: the first one is a complex amplitude contribution, which reflects the amplitudes of the phase dependent electric fields, while the second one is the tensor contribution built-up from the unit electric field vectors of the different interacting waves. Then, in the context of plane wave propagating in the z direction of the laboratory frame (x, y, z), the field tensor can be rewritten as follows:

$$F^{(n-1)} = \|\vec{E}_1\| \dots \|\vec{E}_n\| F^{(n-1)} e^{i\Delta\vec{k}\cdot\vec{z}} \quad (5)$$

where $F^{(n-1)}$ is the unit field tensor defined as,

$$F^{(n-1)} = [\vec{e}^*(\vec{u}_n, \omega_n) \otimes \vec{e}(\vec{u}_1, \omega_1) \otimes \dots \otimes \vec{e}(\vec{u}_{n-1}, \omega_{n-1})] \quad (6)$$

$\vec{e}$ is the unit vector of the electric field $\vec{E}_i$ and $\vec{u}_i$ denotes the unit vector of the wave vector $\vec{k}_i$ of the wave $i$. Depending on the considered applications, the field tensor can be expressed in different coordinates and frames. Using Cartesian coordinates, the polarization states and the associated components of the wave vectors are expressed along the three axes x, y, and z of the laboratory frame. Then any field tensor coefficient can be written as:

$$F_{i_n i_1 \dots i_{n-1}}^{(n-1)} = e_{i_n}^* e_{i_1} \dots e_{i_{n-1}} \quad (7)$$

where $i_1 \dots i_n$ refer to the three laboratory frame axes x, y, z.

In this frame, the Cartesian coordinates of the field tensor coefficients $F_{i_n i_1 \dots i_{n-1}}^{(n-1)}$ is linked to the corresponding electric susceptibility tensor coefficients $\chi_{i_n i_1 \dots i_{n-1}}^{(n-1)}$ both expressed in the same frame with a large number of studies having addressed the symmetry properties of the susceptibility tensor to the symmetry elements of the studied material.[5,6] The Cartesian coordinates are also suited to describe the macroscopic interaction in crystals to be developed in Section 5.

Another frame can be used to describe the field tensor, known as the irreducible tensor base or spherical tensor base.[7,8] This base has been extensively used in atomic and molecular physics however less in optics and nonlinear optics.[9] The advantage of such a spherical approach is to identify invariant quantities through rotations of the reference frame and keep track of them via canonical transformation rules. Earlier work on irreducible tensors in the context of atomic physics led to the definition and development of a suited algebraic toolbox.[10]

Here, we remind the connection between irreducible spherical tensors, defined by $j$ and $m$ coordinates, and Cartesian tensors defined via Cartesian coordinates. First, we start with the elementary transformation of a Cartesian tensor of rank 1, namely a vector with three independent components that will transform into a canonical $(j, m)$ labelled coordinates of the Spherical tensor also with three independent components. The unit vector relation[3,10] is defined as

$$e_m^j = \sum_p (jm|p)e_p \quad (8)$$

with

$$(jm|p) = \begin{pmatrix} i/\sqrt{2} & 1/\sqrt{2} & 0 \\ 0 & 0 & -i \\ -i/\sqrt{2} & 1/\sqrt{2} & 0 \end{pmatrix} \quad (9)$$

with $j = 1, m = +1, 0, -1$ and $p = x, y, z$ and $i$ is the complex number.

Using the spherical tensor notation for a vector in the laboratory frame, and assuming a wave vector along the z-axis, which means an electric field orientation located in the x-y plane, the three basis vectors can be written:

$$\begin{align}
e_1^1 &= \frac{1}{\sqrt{2}}(ix + y) \\
e_0^1 &= -iz \\
e_{-1}^1 &= \frac{1}{\sqrt{2}}(-ix + y) \tag{10}
\end{align}$$

Each of them corresponds to a circularly invariant polarization state of the electric field, namely left- ($m = 1$), right-handed ($m = -1$) circular polarization state and axial ($m = 0$) polarization state. Note that the transformation of a vector into a spherical tensor is said to be irreducible since three independent components are transformed into three independent components in a unique way.

In linear optics, the field tensor is a coupling of two electric fields. Using any of the previous unit vectors, the irreducible spherical tensor is given by:

$$F_M^j = e_{m_1}^{*,1} \otimes e_{m_2}^1 \quad (11)$$

with $m_1, m_2 = (-1, 0, 1)$

In a more explicit way, the field tensor can be generated by way of the vector coupling method (Wigner or CB) to yield:

$$F_M^j = \sum_{p,q} \sum_{\begin{subarray}{c} -1 \le m_1 \le 1 \\ -1 \le m_2 \le 1 \end{subarray}} (j_1 m_1 j_2 m_2 | JM) e_{m_1}^{*,1} e_{m_2}^1 \quad (12)$$

where $(j_1 m_1 j_2 m_2 | JM)$ are the Clebsh-Gordan coefficients.

These coefficients are non-zero when $M = m_1 + m_2$ and $J \le j_1 + j_2$. Using the transformation between the Cartesian and the spherical coordinates basis, Equation (12) becomes:

$$F_M^j = \sum_{p,q} \sum_{\begin{subarray}{c} -1 \le m_1 \le 1 \\ -1 \le m_2 \le 1 \end{subarray}} (j_1 m_1 j_2 m_2 | JM) (j_1 m_1 | p) (j_2 m_2 | q) e_p e_q \quad (13)$$

This equation leads to a relation between the Cartesian components and the spherical components expressed in Table 1,[9]

In a full spherical picture, that is, when the electric fields are expanded over the circular basis and using the orthogonality of the
<table><caption>**Table 1.** Cartesian components as a function of the spherical components of a second-order tensor.</caption><thead><tr><th>J; &plusmn;Mpq</th><th>2;&plusmn;2</th><th>2;&plusmn;1</th><th>2; 0</th><th>1;&plusmn;1</th><th>1; 0</th><th>0;0</th></tr></thead><tbody><tr><td>xx</td><td>-1/2</td><td></td><td>1/&radic;6</td><td></td><td></td><td>1/&radic;3</td></tr><tr><td>yy</td><td>1/2</td><td></td><td>1/&radic;6</td><td></td><td></td><td>1/&radic;3</td></tr><tr><td>xy</td><td>&mpart;1/2</td><td></td><td></td><td>-1/&radic;2</td><td></td><td></td></tr><tr><td>yx</td><td>&mpart;1/2</td><td></td><td></td><td></td><td>1/&radic;2</td><td></td></tr><tr><td>zx</td><td></td><td>&mpart;1/2</td><td></td><td>1/2</td><td></td><td></td></tr><tr><td>xz</td><td></td><td>&mpart;1/2</td><td></td><td>-1/2</td><td></td><td></td></tr><tr><td>zy</td><td></td><td>-1/2</td><td></td><td>&mpart;1/2</td><td></td><td></td></tr><tr><td>yz</td><td></td><td>-1/2</td><td></td><td>&mpart;1/2</td><td></td><td></td></tr><tr><td>zz</td><td></td><td></td><td>-&radic;2/3</td><td></td><td></td><td>1/&radic;3</td></tr></tbody></table>

Clebsh-Gordan coefficients, the link between the electric fields and the tensor field can be expressed as:

$$e_{m_1}^{*,1} e_{m_2}^{1} = \sum_{0 \le J \le 2, -J \le M \le J} (j_1 m_1 j_2 m_2 |JM) F_M^J \quad (14)$$

with $m_1, m_2 = (-1, 0, 1)$. The combination of electromagnetic fields expressed in the spherical base allows to generate linear combination of field tensors of different orders, as given in **Table 2**.

The irreducible frame shows that the first-order field tensor built-up from Equation (14) has only spherical components of degrees $J = 0, 1$ and $2$ with the corresponding orders $M$ following $-J \le M \le J$. When the polarization fields are limited to in-plane polarizations, that is, along the $x$ and $y$ coordinates, the degrees $J = 0, 2$ can only be generated (Table 1). In order to further generate a $F_2^2$ field tensor, the excitation electric field must be of order $m = 1$, which corresponds to a left-handed circular polarization state, and the emitted electric field should be of order $m = 1$ as well. Since the complex conjugate of the emitted electric field has to be considered, it leads to a right-handed circular polarization state. Circular polarizations of the same handedness will generate several tensor fields of different degrees $J$, all with the same $M = 0$ order. Note that the generation of a purely isotropic field tensor $F_0^0$ is not possible, as seen in Table 2. Therefore, the incident polarization states ultimately govern the order $M$, which is

<table><caption>**Table 2.** Spherical tensor components generated from combinations of circular and axial polarization fields. $m = 1$ is the right-handed circular polarization, $m = -1$, the left-handed polarization and $m = 0$, axial/linear polarization.</caption><thead><tr><th>e<sub>m<sub>1</sub></sub> &otimes; e<sub>m<sub>2</sub></sub></th><th>F<sub>2</sub><sup>2</sup></th><th>F<sub>1</sub><sup>2</sup></th><th>F<sub>0</sub><sup>2</sup></th><th>F<sub>1</sub><sup>1</sup></th><th>F<sub>0</sub><sup>1</sup></th><th>F<sub>0</sub><sup>0</sup></th></tr></thead><tbody><tr><td>e<sub>1</sub></td><td>e<sub>1</sub></td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>e<sub>1</sub></td><td>e<sub>0</sub></td><td>1/&radic;2</td><td></td><td></td><td></td><td></td></tr><tr><td>e<sub>1</sub></td><td>e<sub>-1</sub></td><td></td><td>1/&radic;6</td><td></td><td></td><td>1/&radic;3</td></tr><tr><td>e<sub>-1</sub></td><td>e<sub>1</sub></td><td></td><td>1/&radic;6</td><td></td><td></td><td>1/&radic;3</td></tr><tr><td>e<sub>0</sub></td><td>e<sub>0</sub></td><td></td><td>&radic;2/3</td><td></td><td></td><td>1/&radic;3</td></tr><tr><td>e<sub>1</sub></td><td>e<sub>0</sub></td><td></td><td></td><td>1/&radic;2</td><td></td><td></td></tr><tr><td>e<sub>0</sub></td><td>e<sub>1</sub></td><td></td><td></td><td>0</td><td></td><td></td></tr><tr><td>e<sub>1</sub></td><td>e<sub>-1</sub></td><td></td><td></td><td></td><td>-1/&radic;2</td><td></td></tr><tr><td>e<sub>-1</sub></td><td>e<sub>1</sub></td><td></td><td></td><td></td><td>1/&radic;2</td><td></td></tr></tbody></table>

the order of symmetry of the projection of the field tensor spherical representation, onto the plane orthogonal to the field propagation direction $z$. $M$ is the sum of the $m$ values of each field, such that:

$$M = \sum_{n} m_{n} \quad (15)$$

The symmetry of the field tensor will be designated by this value $M$ in the following.

The building principle is the same for the second-order and the third-order field tensors[3,11,12] by adding tensorial products as:

$$F_M^J = e_{m_n}^{*,1} \otimes e_{m_1}^{1} \otimes \dots \otimes e_{m_{n-1}}^{1} \quad (16)$$

Note that irreducible tensors can be more easily visualized by using the associated spherical harmonic functions. Indeed, electric vectors can be expressed as a function of the spherical harmonics $Y_M^L(\theta, \phi)$ as:

$$F_M^L(\theta, \phi) = e_{m_1}^1(\theta, \phi) \times \dots \times e_{m_{n-1}}^1(\theta, \phi) \times e_{m_n}^{*,1}(\theta, \phi) \quad (17)$$

where $L = J$ and with $m_1, \dots, m_n = (-1, 0, 1)$.

The field tensor then becomes the algebraic product of spherical harmonics of order 1, representing each field. This construction makes it possible to go beyond the tensor description and take advantage of an analytical formalism.[13] **Figure 1** shows the tensor fields generated by several combinations of circular polarization states for different multiphotonic processes.

## 2.2. Electric Susceptibility Tensor

The use of the multipolar terminology has been introduced by way of the analogous behaviors of any arbitrarily defined symmetric tensor and that of the classically defined multipole charge tensor, both of the same rank and defined in the same frame. Any moment of order $n$ for an arbitrary distribution could also be used as a reference but charge distribution, including their most reduced form as point charges, has been considered more consistent with electronic charge displacement and even explicitly related to it in early models of nonlinearities.[14,15]

Indeed, a general momentum tensor or multipolar tensor for a $f(\vec{x})$ charge distribution can be expressed jointly in Cartesian and irreducible terms as:

$$M^{(n-1)} = \sum_{i_1, i_2, \dots, i_n} M_{(n-1)}^{i_1, i_2, \dots, i_n} x_{i_1} \otimes x_{i_2} \otimes \dots \otimes x_{i_n} = \sum_{-J \le m \le J} M_{(n-1),J}^{m} e_{J}^{m} \quad (18)$$

$$M_{(n-1)}^{i_1, i_2, \dots, i_n} = \iiint_{\mathbb{R}^{n-1}} f(\vec{x}) x_{i_1} x_{i_2} \cdots x_{i_n} dx_1 dx_2 \cdots dx_n \quad (19)$$

The relations between the Cartesian and irreducible components are the same as discussed earlier.

In general, the irreducible component of $J$ order of any symmetric susceptibility tensor has been referred to a multipole tensor, hence $J = 0$ a scalar charge, $J = 1$ a vector-like dipole, $J = 2$
<figure><img src="image_2.png" /><figcaption>Figure 1. Spherical tensor generated in the cases of second- and third-order processes: a) second-harmonic generation (SHG) energetic diagram (odd orders $F_{M=1}$ and $F_{M=3}$ tensor field are generated by combination of circular polarizations); b) third-harmonic generation (THG) energetic diagram (only $F_{M=2}$ and $F_{M=4}$ can be generated); c) Degenerated four-wave mixing (FWM) (orders $F_{M=0}$, $F_{M=2}$, $F_{M=4}$ can be generated, $F_{M=2}$ from two different possible configurations).</figcaption></figure>

a rank 2 quadrupole, a rank 3 $J = 3$ an octupole, $J = 4$ a rank 4 hexadecapole (or $2^4$-pole) and more generally an arbitrary $J$ to be labelled as a $2^J$-pole. Note that $J$ cannot be higher than the rank of the field tensor, while the rotational spectrum is greatly reduced for symmetric tensors, imposing identical parities for $J$ and $n$, namely $J$ taking only odd (respectively even) values for an odd (respectively even) $n$ rank.

This leads for example to quadratic (second-order) and cubic (third-order) electric susceptibility tensors as well as second- and third-order field tensors to be made-up respectively of two irreducible components, namely a $J = 1$ dipole and a $J = 3$ octupole for the $\chi^{(2)}$ susceptibility and the corresponding $F^{(2)}$ field tensor and three, namely a $J = 0$ scalar charge, a $J = 2$ quadrupole and a $J = 4$, $2^4$-pole for $\chi^{(3)}$ and the corresponding $F^{(3)}$.

Similarly as for a distribution of charges, the tensor formalism and its spherical decomposition has been exploited to express the characteristics of macroscopic electric susceptibility tensors originating from a collection of molecules in the focal volume using of a microscope, in particular in molecular and biologi-

cal media.[16] Supposing a molecular angular distribution $f(\Omega)$ with $\Omega = (\theta, \phi, \psi)$ the Euler set of angles determining the molecular orientation in 3D, the macroscopic electric susceptibility tensor of order $n$ deduced from the molecular electric susceptibility $\chi_{molec}^{(n-1)}$ can be expressed as:

$$ \chi^{(n-1)} = N \int \chi_{(n-1)}^{molec} f(\Omega) d\Omega \quad (20) $$

with $d\Omega = \sin \theta d\theta d\phi d\psi$ and $N$ is the molecular density.

The $\chi_{molec}^{(n-1)}$ molecular electric susceptibility tensor corresponds to the microscopic molecular electric susceptibility, denoted $\beta^{(n-1)}$ in increasing order of optical interactions. Its components can be decomposed onto symmetry elements when expressed in the macroscopic laboratory frame (indices $I_1, I_2 \dots I_n$),[3,16]

$$ \begin{aligned} \chi_{(n-1)}^{molec, I_1 \dots I_n}(\Omega) &= \sum_{i_1 \dots i_n} \chi_{(n-1)}^{molec, i_1 \dots i_n} (I_1 \cdot i_1) \dots (I_n \cdot i_n) (\Omega) \\ &= \sum_{i_1 \dots i_n} \sum_{-J \le m \le J} \chi_{(n)}^{molec, i_1 \dots i_n} c(i, I)_J^m Y_J^m(\Omega) \\ &= \sum_{-J \le m \le J} \chi_{(n-1), J}^{molec, m} Y_J^m(\Omega) \end{aligned} \quad (21) $$

with $(i_1 \dots i_n)$ the microscopic molecular frame and $c(i, I)_J^m$ coefficients of the projection of the $(I_1 \cdot i_1) \dots (I_n \cdot i_n)(\Omega)$ angular functions on the $Y_J^m(\Omega)$ spherical harmonics.[16]

On the other hand, the angular distribution $f(\Omega)$ can be decomposed onto an orthonormal basis of Wigner functions $D_{mm'}^J(\theta, \phi, \psi)$. Supposing that this function is of cylindrical symmetry, it is sufficient to use spherical harmonics $Y_J^m(\theta, \phi)$ with $f(\Omega) = \sum_{m,m'} f_m^J Y_{m'}^m(\Omega)$. Thanks to the orthogonality rule between spherical harmonics, the previous expression (20) can be reduced to:

$$ \chi^{(n-1)} = \sum_{J,m} \chi_{(n-1),J}^{molec,m} f_J^m \quad (22) $$

This leads to a general property applying beyond the case of cylindrical symmetry distributions: the macroscopic electric susceptibility tensor is the result of a symmetry coupling between the microscopic spherical components of the microscopic electric susceptibility tensor and those of the angular distribution of its constituents. As a result, the symmetry orders $J \le n$ of a molecular distribution have the ability to influence the formation of a macroscopic electric susceptibility tensor up to the order ($n-1$), only if they match with the symmetry orders $J$ of this tensor. In even order processes such as linear scattering ($n = 2$ for which $L = 0, 2$) and third order nonlinearities ($n = 4$, for which $L = 0, 2, 4$), the distribution symmetry orders influencing the optical response are therefore even orders, while in odd order processes such as second order nonlinearities ($n = 3$, for which $L = 1, 3$), these symmetry orders are odd, as expected in non-centrosymmetric media which are second-harmonic generation (SHG) active. A large amount of work has been dedicated to biological nonlinear imaging, where the analysis of macroscopic
electric susceptibility $\chi^{(n-1)}$ has allowed the retrieval of information on molecular distribution symmetry orders. This has been applied in SHG in collagen in tissues or in-vitro,[17-19] but also in third-order nonlinear processes such as Coherent Anti Stokes Raman Scattering.[12,20,21]

## 2.3. Light-Matter Interaction

It is convenient in nonlinear optics to introduce the effective co-efficient $\chi_{eff}$ that is calculated as a tensor contraction between the electric susceptibility tensor discussed before and the unit field tensor defined by Equation (6), expressed in an orthonormal frame using the Cartesian or spherical coordinates,[4,9,11,22-25] according to:

$$ \chi_{eff} = \sum_{I...K} \chi_{I...K} F_{I...K} = \sum_{J,M} (-1)^J \chi_M^J F_M^J \quad (23) $$

Of particular relevance in terms of simplification and emergence of tensorial selection rules, it is to be remembered that irreducible subspaces of different order $J$ being orthogonal, the generalized scalar product or contraction of two tensors of the same rank will not scramble different $J$ order.

The amplitude scattered by a single emitter can thus be expressed in terms of a scalar product between the field tensor and the electric susceptibility tensor, the latter one only depending on the emitter orientation described by the following spherical coordinates angles $(\theta, \varphi)$. For a given orientation of the emitter into the laboratory frame, the effective electric susceptibility is then:

$$ \chi_{eff} = \sum_{I...K} \sum_{i...k} (T_i^j(\theta, \varphi) ... T_K^k(\theta, \varphi)) \beta_{i...k} F_{I...K} \quad (24) $$

$$ = \sum_{J,M} \sum_{m} (-1)^J \mathcal{D}_{mM}^J(\theta, \varphi) \beta_m^J F_M^J \quad (25) $$

The quantity $T_i^j(\theta, \varphi)$ is the $iI$ components of the Euler rotation matrix, and $\mathcal{D}_{mM}^J(\theta, \varphi)$ the $mM$ components of the $J$-order Wigner rotation matrix.

At the end, the overall scattered amplitude is the sum of the amplitude scattered by each elementary source, with their respective orientation. If the phase between nonlinear emitters is not preserved during the integration time,[4] the total scattering cross section is $N$ times the cross section of one randomly oriented molecule. In this case, the intensity scattered $I_{inc}$ is found to be incoherent and expressed as:

$$ I_{inc} = NI_0 \langle \chi_{eff} \chi_{eff}^* \rangle \quad (26) $$

where $I_0$ is a constant and $\langle \rangle$ represents the average over all the position and orientation of emitters.[3,26,27]

At the opposite situation of coherent processes, the phase relationship between emitters is conserved, which confers the property to the process to scale nonlinearly with the number of emitters. The coherent intensity $I_c$ of a $n$-order nonlinear process can be then written as[4,13,22]

$$ I_c \propto \langle \chi_{eff} \rangle \langle \chi_{eff}^* \rangle \quad (27) $$

In what follows, Expressions (21) and (22) will be exploited under different conditions that emphasize the role of the field tensor formalism in the understanding of nonlinear optical interactions in molecular environments.

## 3. Microscopic Interaction

The field tensor formalism is a privileged tool for the study of the nonlinear interactions of sum- and difference-frequency generations in the context of microscopy. The nonlinear interaction is significant only inside the focusing volume which does not exceed several wavelengths in length. It is usually assumed that the electric fields are not disturbed by the sample itself, and propagate without any distortion nor depletion. In this context, the field tensor can be reduced to its polarization part $F_M^J$, given by Equation (16), whereas the amplitudes and the phase-matching term are assumed to remain constant. The case where the beams propagate in the nonlinear medium over several millimeters or centimeters length into a non-isotropic medium is discussed in section 5.

### 3.1. Raman Scattering

Raman spectroscopy is currently used to address a very wide range of issues in chemistry and physics. It is now recognized as a valuable technique in biology, being label-free and non-invasive. Individual Raman bands have several characteristic properties such as energy, line shapes, intensity and also polarization states.[28] The polarization of a Raman band is characteristic of the non-zero elements of the associated scattering tensor and allows the assignment of its vibrational symmetry. This method is used qualitatively to distinguish between fully symmetric vibrations modes (that generally scatter highly polarized light, and non-symmetric modes. The molecular electric susceptibility tensor associated to Raman scattering is determined by the symmetry properties of vibrations, which are well documented. The use of polarized light is of great interest to probe these symmetry properties in isotropic media. In oriented media however, the average orientation of molecules makes it more difficult to directly recover the information as the symmetry rules influence the readout at both microscopic and macroscopic scales. Several strategies have been developed to determine the symmetry of vibrations at the molecular scale using circular or linear polarized light. The most elegant and effective method is to concentrate on the irreducible molecular tensor $\chi_M^J$.[29] For each symmetry group, it is necessary to perform the appropriate transformation to obtain the associated isotropic averages, and finally probe it with the field tensor spherical components $F_M^J$.[4,12] The Raman effective coefficient, which embeds both the first-order electric susceptibility and field tensors, is defined as:

$$ \chi_{eff}^R(\omega_r) = \chi_R^{(1)} \cdot F^{(1)} \quad (28) $$

where $F^{(1)} = \vec{e}^*(\omega - \Omega)\vec{e}(\omega)$ and $\chi_R^{(1)}$ is the macroscopic Raman susceptibility tensor at the angular frequency $\Omega$, which contains the molecular symmetry information about the vibrating bonds.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Principle of circularly polarized Raman spectrum. Left, measurement principle with same and opposite circular polarization and field tensors generated in both case.[33] Right, spectrum of CCl₄, with two antisymmetric vibrations and one symmetric, Adapted with permission.[34] Copyright 2024, ACS Publications.</figcaption></figure>

Raman scattering in solutions being an incoherent process, the intensity of Raman scattering follows Equation (21):

$$ 
\begin{aligned}
I_M &= I_0 \chi_{eff}^R \chi_{eff}^{R*} \\
&= I_0 \sum_{J' M' m'} J M m \left( \chi_m^J \chi_{m'}^{*J'} \right) \left( F_M^J F_{M'}^{*J'} \right) \mathcal{D}_{m M}^J \mathcal{D}_{m' M'}^{*J'} \\
&= I_0 \sum_J \frac{1}{2J+1} \left( \chi_M^J \right)^2 \left( F_M^J \right)^2
\end{aligned} 
\quad (29) $$

Equation (29) emphasizes the consequence of orthogonality rules (where the term $\frac{1}{2J+1}$ is a normalization factor) between spherical harmonics: $-2 \le M \le 2$ as the orders of the field tensor range between 0 and 2 by construction as seen in Equation (14) and Table 2. Supposing that only circular polarizations are used in the construction of the field tensor, then M is equal to either 0 or $\pm 2$ (Table 2). Circular polarizations allow therefore to quantify independently the electric susceptibility tensor invariants of order M.[28,30]

In the spectrum obtained from circular Raman spectroscopy, vibrational bands emerge with different amplitudes directly related to the non-zero components of the scattering tensor.[31] The $I_{RR}$ bands generated from right circular polarizations only at both the detection and excitation, also denoted $I_+$, correspond to the symmetric bands ($M=0$), whereas $I_{RL}$ (generated from excitation and detection polarizations of opposite handedness), also denoted $I_-$, corresponds to anti-symmetric bands ($M=2$). An elegant way to quantify the distinction between these different symmetry bonds is to plot the difference $I_+ - I_-$,[32] as depicted in **Figure 2**. The magnitude of the difference $I_+ - I_-$ has also been

used to distinguish totally symmetric from non-totally symmetric modes.[33]

## 3.2. Harmonic Generation

The symmetry coupling described in the previous Section expands to optical effects of higher orders, such as second- and third-order nonlinear coherent processes. 3-wave interactions involving three angular frequencies ($\omega_1, \omega_2, \omega_3$) are governed by the second-order electric susceptibility tensor $\chi^{(2)}$, which covers sum-frequency generation effects SFG ($\omega_3 = \omega_1 + \omega_2$) including second-harmonic generation SHG ($2\omega = \omega + \omega$) and difference-frequency generation DFG ($\omega_1 = \omega_3 - \omega_2$) or DFG ($\omega_2 = \omega_3 - \omega_1$). Four-wave processes involving four angular frequencies ($\omega_1, \omega_2, \omega_3, \omega_4$) are governed by the third-order electric susceptibility tensor $\chi^{(3)}$, and dealing with sum-frequency generation SFG ($\omega_4 = \omega_1 + \omega_2 + \omega_3$), including third-harmonic generation THG ($3\omega = \omega + \omega + \omega$), difference-frequency generation DFG ($\omega_1 = \omega_4 - \omega_2 - \omega_3$) or DFG ($\omega_2 = \omega_4 - \omega_1 - \omega_3$) or DFG ($\omega_3 = \omega_4 - \omega_1 - \omega_2$), and also four-wave mixing FWM ($\omega_4 = -\omega_1 + \omega_2 + \omega_3$) or FWM ($-\omega_1 = \omega_4 - \omega_2 - \omega_3$) or FWM ($\omega_2 = \omega_4 + \omega_1 - \omega_3$) or FWM ($\omega_3 = \omega_4 + \omega_1 - \omega_2$).

### 3.2.1. Second-Harmonic Generation

The lowest-order electric susceptibility $\chi^{(2)}$ is a rank-3 polar tensor which cancels in media with inversion symmetry. It exhibits, in the most general cases, 27 tensor elements, however this number may be drastically reduced by symmetry constraints.[6] The
field tensor associated to this nonlinear process is the tensorial product of the three fields, that is, two of them at $\omega$ the other one being complex conjugated at $2\omega$. [26]

From the irreducible formalism described in Section 2.1, the combination of left- or right-handed circular polarizations generates $F_{M=1}^{(2)}$ and $F_{M=3}^{(2)}$ field tensors, as illustrated in Figure 1. Similarly, $\chi_M^{(2),J}$ only exhibits $M=1,3$ components due to its odd number of involved Cartesian indices, so that the resulting effective electric susceptibility $\chi_{eff}^{(2)}$ involves coupling of odd M values. This result is in adequation with the well-known property of coherent SHG to vanish in isotropic or centrosymmetric materials. Indeed, the isotropic symmetry cannot match the symmetry of the field tensor necessary to “read” SHG processes. The extensive use of the irreducible tensor formalism in molecular crystals has allowed to emphasize the interest of dipolar ($M=1$) and octupolar ($M=3$) symmetries to enhance material-field coupling.[26,27,35] In particular, octupolar materials have the ability to lead to SHG response that are insensitive to the incident polarization state, making them orientation-insensitive, in contrast to materials of dipolar symmetry.[26]

More recently, SHG has permitted to address the symmetry of metamaterials and nano-objects.[36] Behind the complexity of such metallic nano-objects where the local field is driven by Plasmon resonances,[37] some simple symmetry rules emerge. In many cases, the local nonlinear source indeed takes its origin in the symmetry breaking at the metallic interface.[38,39] This implies that the electric susceptibility symmetry is the same as the symmetry of the nanoparticle imposed by its geometry, when both the excitation and emission energies are far from electronic resonances. In the case of a pure dipolar material, a right-handed circular polarization generates harmonic photons with the same polarization. Conversely, a pure octupolar (or a threefolded symmetry) nanostructure will change the direction of rotation of the harmonic photons. In Figure 3a, nanotriangle cavities show a clear threefold symmetry which efficiently matches with a $F_{M=3}$ field tensor, whereas circular nanocavities do not exhibit any SHG.[40,41]

### 3.2.2. Third-Harmonic Generation

Using circular polarization states, the available field tensors for THG are limited to $F_{M=4}^{(3)}$ and $F_{M=2}^{(3)}$. This result also confirms the well-known property of THG, to vanish in isotropic materials when a circular polarization is used.[43] This property of THG that is not sensitive to isotropic distributions under circularly polarized light has also been used in bioimaging to select tissue symmetries and add orientation selectivity.[44] In the plasmonic context, Figure 3b clearly illustrates the sensitivity of THG to the symmetry of different nanoparticles. Nanorods, of second-order symmetry, interact with the circularly polarized fields combination producing $F_{M=2}$ (Right Circular Polarization – Right Circular Polarization) and nano-crosses, of fourth-order symmetry, interact effectively with the $F_{M=4}$ tensor field component (Right Circular Polarization – Left Circular Polarization).[42] This rank-3 process shows its complementarity with SHG, even orders being accessible while odd orders are absent.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Symmetry resolved imaging. a) Field tensor of order $M=1$ and $M=3$ addressed to nano-cavites; Left: nanotriangles show strong up conversion of the signal in counter hand circular polarization; Right: the isotropic symmetry of the nanoholes does not match the field symmetry and no SHG is allowed; scale bar 2 µm. Adapted with permission.[40] Copyright 2024 by the American Physical Society b) HG from linear, three-fold cross and fourfold cross; the RCP-RCP probes the $M=2$ order, mainly in the nanowire; the RCP-LCP probe $M=4$ in the fourfold cross; the three-fold nano-cross are not excited whatever the RCP/LCP configuration; scale bar 400 nm. Adapted with permission.[42] Copyright 2024 by the American Physical Society c) FWM on salt micro-crystal of cubic crystal system; the isotropic order is present in the surrounded medium composed of oil and from the crystal. The order 2 is absent, and the order 4 is only present into the cubic crystal, which highlights the symmetry sensitivity of the technique; Adapted with permission.[11] Copyright 2024 by the American Physical Society. scale bar 10 µm.</figcaption></figure>

### 3.2.3. Four-Wave-Mixing

Similar to THG, FWM processes give access to specific even symmetry orders up to order four. An access to the isotropic order ($M=0$) is only possible when FWM involves non-degenerated
angular frequencies. In this case, isotropic order can be ad-
dressed when all the electric fields display the same circular po-
larization. This behavior differs from THG thanks to different an-
gular frequencies in the energetic path. FWM images of Sodium
Chloride (NaCl) salt microcrystals are depicted in Figure 3c. The
different sets of circular polarizations lead to different informa-
tion: the F<sub>M=0</sub><sup>(3)</sup> field tensor interacts with the isotropic electric sus-
ceptibility from the liquid (oil) and that from the salt crystal; the
F<sub>M=2</sub><sup>(3)</sup> tensor does lead to any signal since there is no twofold sym-
metry in the sample; However, the F<sub>M=4</sub><sup>(3)</sup> field tensor efficiently
matches the electric susceptibility tensor of the NaCl crystal that
belongs to the symmetry O<sub>n</sub> group and has a strong cubic elec-
tric susceptibility.[45] Furthermore, the oil background is invisible
with this last configuration. By using field tensors of pure sym-
metry, it is therefore possible to reveal the symmetry of the elec-
tric susceptibility tensor in one single configuration, leading to a
symmetry-resolved imaging capability.[13]

3.2.4. *High-Harmonic Generation*

Since the advent of ultra-short pulse lasers, HHG in gas has been
at the center of a great interest. It provides a unique light source
in the UV and X-ray ranges, however HHG is strongly affected by
symmetry issues.[46,47] The full tensor description used as for low
harmonic generation, that is, SHG or THG, cannot be exploited
in this configuration, since the number of tensor coefficients dra-
matically increases. In this context, symmetries have been used
to explain the presence or absence of some harmonics in the large
spectrum of the emitted light. Even if the irreducible formal-
ism is not strictly used in this field of research, the use of circu-
lar polarizations has emerged to monitor HHG. HHG has been
generated in a molecular sample,[48–50] and more often in rare
gases, which exhibit an isotropic distribution of molecules. As
in THG whereby isotropic sample cannot be addressed when all
the exciting beams are circularly polarized, circular HHG is for-
bidden in such samples.[48] To circumvent this issue, two beams
at different wavelengths, whose polarization can be monitored
independently, have been focused into the sample.[51–53] In this
scheme, the circular polarization of each beam is adjusted in or-
der to generate a field tensor containing isotropic components.
Hence, the field tensor generated can sustain an isotropic order
that matches the symmetry of the sample. However, the tensor
formalism is questionable in the present context where it is as-
sociated with an electric susceptibility, if not both invalid, in the
case of strong fields much above the internal atomic field, such
as in HHG related phenomena where one is playing a role at
the threshold of ionization. In this context, the perturbative ex-
pansion as such is not correct and may need to be completed by
thermal and squeezed states of light.[54]

**3.3. Coherent Molecular Spectroscopy**

3.3.1. *Coherent Raman Scattering*

Coherent Raman scattering microscopy is based on a third-order
nonlinear process in which a vibrational resonance is reached
through angular frequencies mixing. Because of such a sensi-

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Symmetry resolved molecular spectroscopy. CARS spectrum of a zeolite crystal; two circular polarization configurations probe $M = 0$ and $M = 4$ symmetries; the spectrum shows the symmetry signature of each vibration and its assignment. Adapted with permission.[12] Copyright 2024 by the American Physical Society.</figcaption></figure>

tivity to vibration bonds, this provides a label-free, 3D, chemi-
cally selective and real-time imaging method.[55] In coherent anti-
Stokes Raman scattering (CARS), two beams of different angular
frequencies (pump $\omega_p$ and Stokes $\omega_s$) interact with the sample to
excite a vibrational resonance at the angular frequency $\Omega = \omega_p$
$- \omega_s$. The symmetry of the vibration bond being generally differ-
ent from the crystal symmetry of the corresponding macroscopic
sample, the probe at $\omega_p$ is interacting with the excited vibration
at $\Omega$ to generate a new "blue shifted" photon at $\omega_{as}$ called anti-
Stokes, as shown in **Figure 4**. This CARS process occurs over-
all through a four-wave mixing process with an emission at the
anti-Stokes angular frequency $\omega_{as} = 2\omega_p - \omega_s$. A CARS recorded
spectrum is typically formed by tuning the Stokes or pump fre-
quency, emphasizing the resonance when $\omega_p - \omega_s = \Omega$ depicted
Figure 4.

The symmetry of the intermediate vibrational level probed in
the CARS path has a crucial importance for the symmetry of the
tensor field used. This is illustrated in the example described
in[13] where a Zeolite crystal belonging to the cubic crystal system
was investigated. This molecular crystal has a large number of
atomic vibrations spanning a broad variety of symmetries.[56]
Off any vibrational resonance, the response from the crystal
is mainly isotropic. However, the symmetry orders of the two
vibrations depicted Figure 4 differ drastically from that of the
crystal symmetry. The T<sub>2g</sub> resonance is an anti-symmetric vi-
bration appearing in the spectral response of the CARS process
when the field tensor probes the symmetry order M = 4, using
the appropriate circular polarizations mixture. The Eg resonance
is a two-time degenerated vibration possessing an isotropic
order M = 0 as well as M = 4 order contribution. The M = 0
component of this vibration interacts with the M = 0 component
from the non-resonant response of the crystal, leading to a Fano
profile (e.g., originating from an interference process), as seen in
Figure 4.[12] In this context, the symmetry resolved spectroscopy,
which uses circularly polarized beams, allows the symmetry of
vibrations to be directly identified, bringing a clear added-value
onto the sample. Usually, the samples are less ordered than
crystals, especially in biology. In these cases, the symmetry of
molecular assemblies must consider the molecular distribution
in addition to the symmetry of the vibration, similarly to the
description in Section 2.2. In order to reveal the symmetry of the
sample, it is then necessary to introduce an angular distribution
function representing the probability to find a molecule pointing
in a certain direction. The macroscopic electric susceptibility
tensor in the laboratory frame is the result of the combination of
both molecular electric susceptibilities and of their distribution which can be probed with the same field tensor configuration.

Beyond this organization selectivity, the field tensor formalism has been used to suppress the isotropic background in CARS images and spectra, thus enhancing the contrast by one to two orders of magnitude.[57] It also allows to retrieve quantitative information on the molecular organization, without any post-processing and independently of the sample orientation in the transverse plane.[13]

### 3.3.2. Spinning Molecules

Whereas the CARS process probes the vibrational level of an assembly of molecules, other types of energetic levels can be probed, such as the rotational ones.

A way to probe rotational levels is the rotational Doppler shift (RDS) which is related to the transfer of angular momentum between an anisotropic rotating body and circularly polarized light.[58] In order to populate the rotational level, several methods exist based on optical or physical interactions. For example, molecules have been rotated by optical TeraHz pulses and probed by circular THG.[59] In materials with isotropic symmetry, like gases, THG with circular polarization is forbidden, as seen previously. However, when another process aligns the molecules, and therefore populates a rotational level, the process is allowed again. RDS can be considered as a multiphoton effect where the first two photons are coming from a first excitation to prepare the alignment of molecules. Then a THG or HHG can occur.[60] Figure 5 illustrates the multiphoton process using circular polarizations. In the case of isotropic materials, such as randomly oriented molecules, the electric susceptibility belongs to the isotropic order which is inaccessible with a field tensor of order 2 or 4. When the molecules are oriented with a two-photon process, the electric susceptibility tensor of the sample belongs to the order $M = \pm 2$, or higher even orders. From there, the molecule can convert three photons at $\omega$ onto one photon at $3\omega$. The scheme of Figure 5 describes the THG as a fifth-order nonlinear process. The THG signal is blue shifted for molecules rotating in the same direction as circular polarization of the electric field, and red shifted for counter-rotating molecules.

As a way to cancel a nonlinear process into an isotropic medium, the circular THG offers a clear picture of rotational Doppler Effect by way of a multiphotonic process.

## 4. Multi-Poling Field Tensor

### 4.1. Induced Poling processes

Imprinting order in an initially isotropic media is of crucial importance to allow for rank-3 tensors such as attached to SHG and other three-wave-mixing phenomena. In general, a sequence of two processes is required: a first one to print a single dipole, or more generally as will be further developed, to multipole an initially amorphous or isotropic media such as a liquid, polymer or ferroelectric medium, and a second one to probe the emergence and eventually decay of a previously or simultaneously imprinted

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Multiphotonic picture of the rotational Doppler shift. Circular THG is forbidden in static gas but becomes available when molecules rotate due to 2-photon orientation. Adapted with permission.[59] Copyright 2024 by the American Physical Society.</figcaption></figure>

(multi)polar order. Tensor considerations have played from the beginning a central role in poling processes, with irreducible expressions shown to considerable rationalize derivations and lending themselves to generalizations that may have been obscured by the earlier purely Cartesian expressions. In that respect, idiosyncratic loops had locked molecular engineering of materials to the sole pool of dipolar molecules until the mid-80's, whereas broader tensorial considerations have allowed to enlarge the domain to octupolar molecules for SHG or higher multipolar orders for two-photon absorption. This has triggered to this day an active search for multipolar methods of both electric, optical or mixed nature toward their orientation. A basic challenge in the domain had originally been the impossibility to upscale and measure the octupolar component of a molecule as it cannot be coupled to a dipolar field due to tensorial orthogonality of the $J = 3$ and $J = 1$ irreducible spaces.

Two main generic configurations have been investigated in-depth, both theoretically and experimentally: one relates to equilibrium "writing" conditions whereby a Maxwell-Boltzmann statistic can be applied[27]; the other one relates to out-of-equilibrium dynamical processes where transient states are participating by way of quantum transitions at the molecular level and iterative cycling processes, until an equilibrium is reached.[27,61]

Both situations feature a sequential combination of write- and read-multipolar field tensors, in conjunction with the nonlinear electric susceptibility at stake. This write- and read- field tensors and associated experiments are the equivalent of pump-probe experiments in time-resolved spectroscopy. This similarity has been purposely used in[62] where a first write sequence imprints a
non-centrosymmetric grating in a viscous fluid, to be followed by a probe field tensor diffracted by the grating, the whole scheme pertaining to a six-wave-mixing configuration. A wealth of studies has been dedicated to oriented liquids and polymers for SHG, the former under the label of EFISH, for electric field induced second-harmonic.[27,26] EFISH has remained throughout decades the gold standard for the measurement of the quadratic nonlinear tensor of dipolar molecules. Dipolar orientation of guest-host polymers has been a domain of intensive investigations in the perspective of Pockels electro-optic modulation and their possible large-scale implication in broad-band optical telecommunication systems. In both cases, the write tensor is made-up of a single static or low angular frequency field that is generated by planar capacitor types of electrodes sandwiching the medium. In such situations whereby the field-molecule coupling is weak compared to intermolecular interactions, the write tensor reduces to a single constant $\vec{E} = E^{J=1}$ vector component that can only couple to the corresponding $\beta^{J=1}$ component of the molecular quadratic hyperpolarizability $\beta$.

There are two important exceptions to this configuration: that of planar interdigitated electrodes for waveguided applications as in the case of periodically-poled LiNbO₃.[63] In this medium, the field distribution is more complex, in particular at domain interfaces as seen by using Pockels linear electro-optic microscopy (PLEOM) imaging[64] requiring the consideration of multipolar J orders higher than 1. Another important configuration consists in breaking both the weak-field and homogenous conditions of poling, thus involving a writing-field tensor consisting higher multipolar orders beyond J = 1. Such conditions are met in nanoscale experiments where a non-homogeneous write-field is applied over very short distances, possibly down to the nanometer scale via a scanning tunneling microscope (STM) tip leading to strong field molecule interactions that invalidate the weak field assumption.[65] Finally, polymer poling in a liquid crystal environment would also preclude a weak field type coupling assumption and require the consideration of higher order components in the write-tensor.

As introduced in Section 2.2, the macroscopic tensor $\chi^{(n-1)}$ is obtained from the statistical superposition of molecular electric susceptibilities $\beta^{(n-1)}$ averaged over the orientational molecular distribution function $f(\Omega)$, that is:

$$ \langle \chi^{(n-1)} \rangle_{\Omega} = \frac{N}{kT} \int R_{\Omega}(\beta^{(n-1)}) f(\Omega) d\Omega \quad (30) $$

where $R_{\Omega}$ is a rotation operator and $f(\Omega)$ originates from the Boltzmann distribution

$$ f(\Omega) = \exp \left[ -\frac{W(\Omega)}{kT} \right] \quad (31) $$

where T, associated to the Boltzmann constant, is the temperature and the potential $W(\Omega) = W_0 - T.E$ occurring from the interaction between the writing-field E and the molecular tensor T.T is a general coupling tensor best defined in the molecular frame with J components eventually beyond the sole J = 1 dipolar order.

The writing-field tensor E is defined along the suited axis of the laboratory frame, thus requiring rotation of T to bring it in

the same frame, so that:

$$ f(\Omega) = \exp \left[ -\frac{W(\Omega)}{kT} \right] \approx 1 - \frac{W_0}{kT} + \frac{1}{kT} \sum_{m,J} (-1)^m R_{\Omega}(T) E_J^m \quad (32) $$

then

$$ \langle \chi^{(n-1)} \rangle_{\Omega} \cdot F = \frac{N}{kT} \sum_{J=0 \text{ or } 1}^{n} \frac{1}{2J+1} (T^J \cdot \beta^{(n-1)J}) E^J \cdot F^J \quad (33) $$

In this last general expression, the contracted product of the irreducible components $E^J$ and $F^J$ of the write-field tensor with the read-out electric field tensor is weighted by the contracted product at the same J order of the coupling tensor T with the molecular hyperpolarizability $\beta^{(n-1)}$ tensor. Summation over J is limited by parity considerations on the rank n as mentioned earlier in the case of symmetric tensors. The detailed derivation for this expression is based on a fundamental tensor permutation scheme valid within each J-labelled irreducible sub-space,[27] that is:

$$ (T^J \cdot E^J) \beta^{(n-1)J} = (T^J \cdot \beta^{(n-1)J}) E^J \quad (34) $$

In the special cases of quadratic processes with respectively n = 3, it comes:

$$ \chi^{(2)} = \sum_{J=1,3} \chi^{(2)J} = \frac{N}{kT} \sum_{J=1,3} \frac{1}{2J+1} (T^J \cdot \beta^J) E^J \quad (35) $$

In order to fully upscale at the macroscopic level the full rotational spectrum of $\beta$, these expressions show that the coupling term needs to match the J = 1 and 3 rotational spectrum of the nonlinear susceptibility either. It therefore appears from Equation. (35) that a purely dipolar write-field tensor, namely $T = T^1$ while allowing to orient the corresponding J = 1 vector component $\beta^J$ will not allow for the corresponding orientation of the J = 3 octupolar component, thus requiring more elaborate poling schemes. This can be ensured by way of strong coupling conditions or via and adequately designed set of multipolar, in the case octupolar electrodes both at the nanoscale of an octupoling and more generally a multipoling scheme.[66,67]

In the case of third-order processes with n = 4, it comes:

$$ \chi^{(3)} = \sum_{J=0,2,4} \chi^{(3)J} = \frac{N}{kT} \sum_{J=0,2,4} \frac{1}{2J+1} (T^J \cdot \gamma^J) E^J \quad (36) $$

Likewise, orientational upscaling of the full rotational spectrum of the cubic susceptibility $\gamma$ requires the involvement of a coupling tensor T comprising the same J = 0, 2, 4 orders as for $\gamma$.

The above expressions (35) and (36) can be specialized to the more classical case of dipolar coupling in both EFISH and poled polymer cases, with the poled molecular dipole moment for the coupling, namely $T = T^1 = \vec{\mu}$. Dynamical studies of build-up and relaxation of poled polymers have been performed by way of alternating electric write- and optical read-tensors within adequately defined duty cycles.[68]
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Imprint of a statistical order in a DR1-MMA functional polymer for various write-read field tensor configurations. Adapted with permission.[27] Copyright 2024 Optica Publishing Group. Polar plot $|e^{2\omega}|(\phi)$ of the second-harmonic signal read out by linearly polarized fundamental beam (with the polarization angle $\phi$) in the case of grafted DR1-MMA spin-coated films. Squares correspond to experimental values; solid curves to theoretical fit. a) Analyzer at $2\omega$ along X, $E^{2\omega}$ and $E^\omega$ writing fields linearly polarized along X. b) Analyzer at $2\omega$ along Y, $E^{2\omega}$ and $E^\omega$ writing fields linearly polarized along Y. c) Analyzer at $2\omega$ along X, $E^{2\omega}$ and $E^\omega$ writing fields linearly polarized along respectively Y and X. d) Analyzer at $2\omega$ along Y, $E^{2\omega}$ and $E^\omega$ writing fields linearly polarized along respectively Y and X.</figcaption></figure>

## 4.2. All-Optical Poling Processes

All-optical processes involve in general a sequential or inter-
twinned cooperation of write and read field tensors as it will be
derailed later for different configurations. We therefore revisit
and generalize our earlier definition of the field tensor as a mul-
tiple field tensor so as to accommodate such more general con-
figurations, as developed in Appendix B.

The possibility of inducing a permanent $\chi^{(2)}$ by light in a
centrosymmetric material has been originally discovered by
serendipity in silica fibers.[69] It used the irradiation of a single
laser beam and was subsequently optimized by simultaneous
seeding of harmonic light at the doubled frequency.[40] The origin
of the ensuing SHG process has been phenomenologically ac-
counted for by the non-zero average of the write-field tensor over
a period $T = 2\pi/\omega$, namely $\langle \tilde{E}^{2\omega} \otimes \tilde{E}^\omega \otimes \tilde{E}^\omega \rangle_T \neq 0$.[70] The quan-
tum origin of this symmetry breaking process had been shown
by R. Glauber back in the 70's[71] to be governed by interferences
between one- and two-photon pathways, then in the context of
statistical correlation studies in quantum optics. In photosen-
sitive polymers abiding to molecular reorientation under cycles
of randomizing cis-trans-photoisomerization, like PMMA-DR1
shown in **Figure 6**, it has been first demonstrated by[41,42]
in a dipolar scheme. Further it has been extended by[27] and

subsequent papers for the full molecular and orientational or-
ders, including the otherwise elusive $J = 3$ octupolar component
shown to be orientable under an adequately defined all-optical oc-
tupoling scheme. Microscopically, the photo-induced excitation
of molecules, coupled to possible transient photo-isomerization
processes and other possible intramolecular or intermolecular
mechanisms, results in an angularly selective randomization
process leading to an orientationally selective hole-burning with
the possibility to engineer its multipolar symmetry at will. This
dynamical process has been studied in[72,73] based on alternation
of write- and read-tensor exposure of the sample.

The introduction of the irreducible tensor formalism high-
lights the rotational symmetry of the core-molecules via their
nonlinear hyperpolarizability tensor. It is in conjunction with the
orientating mechanism via the adequate transition interference
tensor that is in general amenable to an equivalent hyperpolar-
izability tensor and last of both the write- and the read-tensor
fields.[14,21,44] Here it has been assumed that an orientational hole-
burning process is responsible for the imprint via partial ran-
domization of the molecules. The depth of the orientational hole
burnt in the initially isotropic distribution of photo-isomerizable
nonlinear molecules scales, with the probability of activating the
reorienting trans-cis transition. A simple semi-classical quantum
model of joint one- and n-photon transition probability $P_{01}(\Omega)$
within a two-level quantum molecular model accounts for the emergence of the non-linearity $\beta^{(n-1)}$ appearing again, however in relation with a multiphoton transition:

$$f(\Omega) \propto P_{01}(\Omega) \propto \beta^{(n-1)} \cdot E \qquad (37)$$

$$\langle \chi^{(n-1)} \rangle_{\Omega} \propto \int R_{\Omega}(\beta^{(n-1)}) f(\Omega) d\Omega \propto \left[ \int R_{\Omega}(\beta^{(n-1)}) \beta^{(n-1)} d\Omega \right] \cdot E \qquad (38)$$

with the writing field E now occurring from multiple optical fields:

$$E = (\vec{E}^{n\omega})^* \otimes \vec{E}^{\omega} \cdots \otimes \vec{E}^{\omega} \qquad (39)$$

The writing field E is then a field tensor rank n comprising the (n-1) harmonic and n-1 times the fundamental field. Finally, along similar irreducible tensor expressions and their canonical orientational transformation properties it comes:

$$\langle \chi^{(n-1)} \rangle_{\Omega} \cdot F \propto \sum_{j=0}^{n} \frac{1}{2J+1} (\beta^{(n-1)J} \cdot \beta^{(n-1)J}) E^j \cdot F^j \qquad (40)$$

Equation (40) looks, at least at a formal level, strikingly similar to the corresponding one for multipoling at equilibrium under a Maxwell-Boltzman statistic. However, it is of a fundamentally different physical origin. It is most striking to notice here the same $\beta^{(n-1)}$ tensor showing-up as a nonlinear electric susceptibility and also as a multiphoton transition rate for parallel mutually coherent interfering one- and n-photon pathways. The E write-field tensor provides a powerful control on the orientation of the medium by way of its rotational spectrum.

In a transverse planar configuration with beams perpendicular to a planar sample and in-plane transverse field polarizations, the double integer (J, m) labelling can be reduced to $m = \pm J$ very much like general $Y_m^J(\theta, \phi)$ spherical harmonics are amenable to $e^{im\phi}$ complex exponentials. Splitting the fundamental beam at frequency $\omega$ in n circularly polarized beams of individually controlled handedness labelled by $\epsilon_i\omega$, with $\epsilon_i = \pm 1$ with $1 \le i \le n$ and $\epsilon_0 = \pm 1$ for the $n^{th}$- harmonic, allows to span the full rotational spectrum for the write-tensor and therefore promote any $J = \pm m$ component of the molecular hyperpolarizability via the following expression of the write field:

$$\begin{aligned} E &= (E^{\epsilon_n n \omega})^* \otimes E^{\epsilon_1 \omega} \cdots \otimes E^{\epsilon_{n-1} \omega} \\ &= E_1 \cdots E_n e^{i(-\epsilon_n + \epsilon_1 + \cdots + \epsilon_{n-1})} \phi \vec{c}_*^{\epsilon_n} \otimes \vec{c}_{\epsilon_1} \cdots \otimes \vec{c}_{\epsilon_{n-1}} \\ &\quad - n \le \epsilon_1 + \cdots + \epsilon_n \le n \\ \vec{c}_{\epsilon_i} &= \frac{1}{\sqrt{2}} (\vec{x} + \epsilon_i \vec{y}) \end{aligned} \qquad (41)$$

where the $\vec{c}_{\epsilon_i}$ stand for the normalized right or left-handed circularly polarized complex unit vectors. The more general case, where the transition is governed by the interference between p- and q- multiphoton pathways, would introduce other degrees of control. Its derivation would be technically tedious but follows essentially the same path as above for $p = 1$ and $q = n$.

The following simpler expression applies to the case of quadratic nonlinear properties, sustained by a molecular hyperpolarizability tensor with both $J = 1$ dipolar and $J = 3$ irreducible components, that is:

$$\langle \chi^{(n-1)} \rangle_{\Omega} \cdot F \propto \sum_{J=1,3} \frac{1}{2J+1} \|\beta^J\|^2 E^J \cdot F^J \qquad (42)$$

where F is the field tensor of Equation (4). The $\beta^J \cdot \beta^J$ contracted product is expressed here as the square of the tensor norm $\beta^J$ norm. In order to single out the dipolar (respectively octupolar) component of $\beta$, it is necessary to make use of the corresponding dipolar (respectively octupolar) write- and read-field tensors. In a planar configuration with in-plane write- and read-field tensors, the write-field tensor can be expressed on a basis of complex exponential with integer exponents as following:

$$E^J = (\vec{E}^{\epsilon_3 2\omega})^* \otimes \vec{E}^{\epsilon_1 \omega} \otimes \vec{E}^{\epsilon_2 \omega} = E_1 E_2 E_3 e^{i(-\epsilon_3 + \epsilon_1 + \epsilon_2)} \phi \vec{c}_*^{\epsilon_3} \otimes \vec{c}_{\epsilon_1} \otimes \vec{c}_{\epsilon_2} \qquad (43)$$

The writing field tensor can be engineered by different combinations of circular polarizations of the harmonic and fundamental beams to sustain an adjustable balance of dipolar or octupolar symmetries, namely $\epsilon_1 = \epsilon_2 = \mp 1$, $\epsilon_3 = \pm 1$ for dipolar symmetry and $\epsilon_3 = \mp 1$, $\epsilon_1 = \epsilon_2 = \pm 1$ for octupolar symmetry.[74] A reading field tensor and analyzer with matching parities are to be further used to probe the dipolar and octupolar components of the imprinted $\chi^{(2)}$. Note that there exists the possibility of chiral orientation for the octupolar component by inverting the parity of the fields contributing to the write-field tensor.

As an example, a sample of DR1-methyl methacrylate was excited with collinear circular or counter circular ($\omega, 2\omega$) polarizations, as shown in Figure 6.[75] The $E^1$ field tensor component produces in the collinear ($\omega, 2\omega$) circular configuration imprints a dipolar statistical order resulting in a macroscopic dipolar response of the sample. The read-field tensor then reveals the dipolar response that can be identified by its two lobes on the polarization plot (Figure 6a). In the counter circular ($\omega, 2\omega$) configuration, the $F_3$ “writing” field tensor imprints an octupolar symmetry order into the sample, as seen in Figure 6c via its isotropic polarization response. In the general topic of all-optic poling, the purely octupolar configuration is of particular interest, as it is the unique configuration which is capable of generating a polarization-independent harmonic response.

## 5. Macroscopic Interaction

The field tensor formalism is also a privileged tool for the study of the second- and third-order nonlinear interactions of sum- and difference-frequency generations at the macroscopic scale in crystals. The configuration of polarization is of prime importance as in the case of the microscopic scale, but it becomes necessary to take into propagative considerations when the interacting waves propagate in the nonlinear medium over several millimeters or centimeters, a configuration which is of prime interest when phase-matching conditions are fulfilled. It is therefore our aim in the following paragraphs to illustrate the role and potential of the field tensor to investigate such configurations.
## 5.1. Optical Classes and Eigenmodes of Polarization in Crystal Optics

In this context, the field tensor is built from the unit electric field vectors corresponding to the eigenmodes of polarization of the different interacting waves. These polarization states are imposed by the direction of propagation $\vec{u}$ that is considered. Situations differ depending on the crystal of interest that can be long to the isotropic optical class or to the two anisotropic optical classes, labeled uniaxial and biaxial. It is based on relative values between the three principal refractive indices $n_x(\omega) = \sqrt{1 + \chi_{xx}^{(1)}(\omega)}$, $n_y(\omega) = \sqrt{1 + \chi_{yy}^{(1)}(\omega)}$ and $n_z(\omega) = \sqrt{1 + \chi_{zz}^{(1)}(\omega)}$ where $\omega$ refers to the angular frequency and $(x, y, z)$ stands for the dielectric frame, that is the frame that diagonalizes the first-order electric susceptibility tensor $\chi^{(1)}(\omega)$ that is a rank-3 polar tensor. The optical classes are derived from the application of Neumann's principle to $\chi^{(1)}$. This principle stipulates that if matter remains invariant according to several orientation symmetry operations, then any physical property of matter has to be invariant according to the same orientation symmetry operations. As a consequence, the tensor of interest, that is $\chi^{(1)}$, has to be invariant according to all the orientation symmetry operations of the point group of matter. Doing this for the 32 point groups leads to three possible forms for the matrix representing $\chi^{(1)}$, with the following symmetries when the matrices are expressed in the dielectric frame: $\chi_{xx}^{(1)} = \chi_{yy}^{(1)} = \chi_{zz}^{(1)}$ for the isotropic optical class, which corresponds to crystals belonging to the cubic crystal classes, but also to gas, liquids and glasses; $\chi_{xx}^{(1)} \neq \chi_{yy}^{(1)} \neq \chi_{zz}^{(1)}$ for the uniaxial anisotropic optical class to which the trigonal, tetragonal and hexagonal crystals belong; and $\chi_{xx}^{(1)} \neq \chi_{yy}^{(1)} \neq \chi_{zz}^{(1)}$ for the biaxial anisotropic optical class corresponding to triclinic, monoclinic and orthorhombic crystals. Then it comes for the refractive indices: $n_x(\omega) = n_y(\omega) = n_z(\omega) \equiv n(\omega)$ for the isotropic optical class; $n_x(\omega) = n_y(\omega) = n_o(\omega) \neq n_z(\omega) \equiv n_e(\omega)$ for the uniaxial anisotropic optical class, where $n_o(\omega)$ and $n_e(\omega)$ are the ordinary and extraordinary principal refractive indices, respectively; and $n_x(\omega) \neq n_y(\omega) \neq n_z(\omega)$ for the biaxial anisotropic optical class. Let $\vec{e}(\omega, \vec{u})$ be the unit electric field vector describing the linear polarization state of a wave propagating in a direction $\vec{u}(\theta, \varphi)$, where $(\theta, \varphi)$ are the angles of spherical coordinates linked to the Cartesian coordinates in the dielectric frame.

The calculation of $\vec{e}(\omega, \vec{u})$ can be done from solving the wave propagation equation projected on the three axes of the dielectric frame.[6] It can be shown that any direction of polarization is allowed in media of the isotropic optical class. The case of the two anisotropic optical classes, for which there are only two possible directions of polarization imposed by the direction of propagation $\vec{u}(\theta, \varphi)$, is completely different. In the cases of any direction of propagation in a uniaxial crystal, the two eigenmodes of polarization are given in terms of ordinary and extraordinary modes, respectively $\vec{e}^0(\omega, \vec{u})$ and $\vec{e}^e(\omega, \vec{u})$ and are expressed in the dielectric frame as:

$$ \vec{e}^0(\omega, \vec{u}) = \begin{pmatrix} e_x^0(\omega, \vec{u}) \\ e_y^0(\omega, \vec{u}) \\ e_z^0(\omega, \vec{u}) \end{pmatrix} = \begin{pmatrix} -\sin(\varphi) \\ \cos(\varphi) \\ 0 \end{pmatrix} \quad (44) $$

and

$$ \vec{e}^e(\omega, \vec{u}) = \begin{pmatrix} e_x^e(\omega, \vec{u}) \\ e_y^e(\omega, \vec{u}) \\ e_z^e(\omega, \vec{u}) \end{pmatrix} = \begin{pmatrix} -\cos(\varphi) \cos(\theta \pm \rho^e(\omega, \vec{u})) \\ -\sin(\varphi) \cos(\theta \pm \rho^e(\omega, \vec{u})) \\ \sin(\theta \pm \rho^e(\omega, \vec{u})) \end{pmatrix} \quad (45) $$

with

$$ \rho^e(\omega, \vec{u}) = \mathrm{Arcos} \left[ \frac{n_{-2}^0(\omega) \cos^2(\theta) + n_{-2}^e(\omega) \sin^2(\theta)}{\sqrt{n_{-4}^0(\omega) \cos^2(\theta) + n_{-4}^e(\omega) \sin^2(\theta)}} \right] \quad (46) $$

The sign ($\pm$ or $+$) in Equation (45) is ($+$) for a negative uniaxial crystal, and ($-$) for a positive one. The quantity $\rho^e(\omega, \vec{u})$ is the spatial walk-off. Note that $\rho^e(\omega, \vec{u}) = 0$ in the $xy$-plane.

Equation (44) shows that the ordinary polarization is orthogonal to the $z$-axis that is called the optical axis of the uniaxial crystal. The extraordinary polarization is then orthogonal to the ordinary polarization, and so parallel to a plane including the optical axis according to Equation (45).

The ordinary and extraordinary eigenmodes of polarization are also relevant in the case of a propagation in the principal planes of the dielectric frame of biaxial crystals. But it is necessary to replace $(n_o, n_e)$ by $(n_x, n_z)$ in the $xz$-plane ($\varphi = 0$), by $(n_y, n_z)$ in the $yz$-plane ($\varphi = \pi/2$), and by $(n_y, n_x)$ in the $xy$-plane ($\theta = \pi/2$) in Equations (45) and (46) for the extraordinary polarization. Equation (44), that describes the ordinary polarization, can be also used in the $xz$- and $yz$-planes of biaxial crystals provided the previous modifications are done, but not in the $xy$-plane where it has to be replaced by:[6]

$$ \left\{ \begin{aligned} \vec{e}^0(\omega, \vec{u}) &= \begin{pmatrix} -\sin(\varphi \mp \rho^0(\omega, \vec{u})) \\ \cos(\varphi \mp \rho^0(\omega, \vec{u})) \\ 0 \end{pmatrix} \\ \rho^0(\omega, \vec{u}) &= \mathrm{Arcos} \left[ \frac{n_y^{-2}(\omega) \cos^2(\varphi) + n_x^{-2}(\omega) \sin^2(\varphi)}{\sqrt{n_y^{-4}(\omega) \cos^2(\varphi) + n_x^{-4}(\omega) \sin^2(\varphi)}} \right] \end{aligned} \right. \quad (47) $$

The sign ($-$ or $+$) in Equation (47) is ($-$) for a negative uniaxial crystal, and ($+$) for a positive one.

For a propagation out of the principal planes of biaxial crystals, it is necessary to use the following equation:

$$ e_i^{\pm}(\omega, \vec{u}) = \frac{[n^{\pm}(\omega, \vec{u})]^2}{[n^{\pm}(\omega, \vec{u})]^2 - [n_i(\omega)]^2} u_i [\vec{u} \cdot \vec{e}^{\pm}(\omega, \vec{u})] \quad (48) $$

where $i = x, y, z$.

$n^{\pm}(\omega, \vec{u})$ refers to the two possible values that the refractive index can take in the direction of propagation $\vec{u}$ that is considered. It is expressed by the two solutions of the index surface equation,[6] that is,:

$$ n^{\pm}(\omega, \vec{u}) = \sqrt{\frac{2}{-B \mp \sqrt{B^2 - 4C}}} \quad (49) $$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Index surfaces and ordinary (o) and extraordinary (e) eigenmodes of polarizations in the cases of positive uniaxial (left) and biaxial (right) crystals. OA denotes the optical axis and (x, y, z) is the dielectric frame. V is the angle between the optical axis of the biaxial crystal and the z-axis.</figcaption></figure>

with

$$ 
\begin{aligned}
B &= -u_x^2 (b + c) - u_y^2 (a + c) - u_z^2 (a + b) \\
C &= u_x^2 bc + u_y^2 ac + u_z^2 ab \\
a &= n_x^{-2} (\omega), \quad b = n_y^{-2} (\omega), \quad c = n_z^{-2} (\omega)
\end{aligned} 
\quad (50) 
$$

Equations (48–50) allow us to establish that $\vec{e}^+(\omega, \vec{u})$ and $\vec{e}^-(\omega, \vec{u})$ are never orthogonal out of the principal planes of biaxial crystals. Note that the quantity $n^+(\omega, \vec{u}) - n^-(\omega, \vec{u})$ is the birefringence in the direction $\vec{u}$ and that the relation of order between the principal refractive indices leads to the definition of the optical sign of uniaxial and biaxial crystals: the sign is positive if $n_0(\omega) < n_e(\omega)$ for uniaxial crystals, and if $n_x(\omega) < n_y(\omega) < n_z(\omega)$ for biaxial crystals. The inequalities are reversed for a negative optical sign. Since the eigenmodes are linear polarizations, then they are described by real quantities.

As an example, the index surface with the associated eigenmodes of polarization of positive uniaxial and biaxial crystals are shown in **Figure 7**.

## 5.2. Symmetries of Field Tensors

The relations between the field tensor coefficients at different angular frequencies are built up from the commutativity of the algebra product of the unit electric field vectors Cartesian coordinates that compose the field tensors. Below are detailed their expression for all the possible 3-wave and 4-wave nonlinear optical processes, coming from Equation (6), where all the electric fields are real quantities, that is, $\vec{e}^*(\vec{u}_n, \omega_n) = \vec{e}(\vec{u}_n, \omega_n)$, since they correspond to linear polarizations described by Equations (44)-(48):

* for SFG/DFG 3-wave mixing,[23]

$$ 
\begin{aligned}
F_{ijk} (\omega_3 = \omega_1 + \omega_2, \vec{u}) &= e_i (\omega_3, \vec{u}) e_j (\omega_1, \vec{u}) e_k (\omega_2, \vec{u}) \\
&= F_{jik} (\omega_1 = \omega_3 - \omega_2, \vec{u}) = e_j (\omega_1, \vec{u}) e_i (\omega_3, \vec{u}) e_k (\omega_2, \vec{u}) \\
&= F_{kij} (\omega_2 = \omega_3 - \omega_1, \vec{u}) = e_k (\omega_2, \vec{u}) e_i (\omega_3, \vec{u}) e_j (\omega_1, \vec{u})
\end{aligned} 
\quad (51) 
$$

* for SFG/DFG 4-wave mixing,[24]

$$ 
\begin{aligned}
F_{ijkl} (\omega_4 = \omega_1 + \omega_2 + \omega_3, \vec{u}) &= e_i (\omega_4, \vec{u}) e_j (\omega_1, \vec{u}) e_k (\omega_2, \vec{u}) e_l (\omega_3, \vec{u}) \\
&= F_{jikl} (\omega_1 = \omega_4 - \omega_2 - \omega_3, \vec{u}) = e_j (\omega_1, \vec{u}) e_i (\omega_4, \vec{u}) e_k (\omega_2, \vec{u}) e_l (\omega_3, \vec{u}) \\
&= F_{kijl} (\omega_2 = \omega_4 - \omega_1 - \omega_3, \vec{u}) = e_k (\omega_2, \vec{u}) e_i (\omega_4, \vec{u}) e_j (\omega_1, \vec{u}) e_l (\omega_3, \vec{u}) \\
&= F_{lijk} (\omega_3 = \omega_4 - \omega_1 - \omega_2, \vec{u}) = e_l (\omega_3, \vec{u}) e_i (\omega_4, \vec{u}) e_j (\omega_1, \vec{u}) e_k (\omega_2, \vec{u})
\end{aligned} 
\quad (52) 
$$

* for FWM 4-wave mixing:

$$ 
\begin{aligned}
F_{ijkl} (\omega_4 = -\omega_1 + \omega_2 + \omega_3, \vec{u}) &= e_i (\omega_4, \vec{u}) e_j (\omega_1, \vec{u}) e_k (\omega_2, \vec{u}) e_l (\omega_3, \vec{u}) \\
&= F_{jikl} (-\omega_1 = \omega_4 - \omega_2 - \omega_3, \vec{u}) = e_j (\omega_1, \vec{u}) e_i (\omega_4, \vec{u}) e_k (\omega_2, \vec{u}) e_l (\omega_3, \vec{u}) \\
&= F_{kijl} (\omega_2 = \omega_4 + \omega_1 - \omega_3, \vec{u}) = e_k (\omega_2, \vec{u}) e_i (\omega_4, \vec{u}) e_j (\omega_1, \vec{u}) e_l (\omega_3, \vec{u}) \\
&= F_{lijk} (\omega_3 = \omega_4 + \omega_1 - \omega_2, \vec{u}) = e_l (\omega_3, \vec{u}) e_i (\omega_4, \vec{u}) e_j (\omega_1, \vec{u}) e_k (\omega_2, \vec{u})
\end{aligned} 
\quad (53) 
$$

Equations (51–53) show that the field tensor coefficients remain unchanged by concomitant permutations of the unit electric field vectors coordinates at the different angular frequencies and the corresponding Cartesian indices. These relations are equivalent to thermodynamic symmetry regarding the electric susceptibility tensor in the case of a lossless medium, that is, when the electric susceptibility coefficients are all real numbers.[76] It is illustrated by the following equalities in the case of a 3-wave mixing process:

$$ \chi_{ijk} (\omega_3 = \omega_1 + \omega_2) = \chi_{jik} (\omega_1 = \omega_3 - \omega_2) = \chi_{kij} (\omega_2 = \omega_3 - \omega_1) \quad (54) $$

The effective coefficients of SFG and DFG can be then derived from Equations (51) and (54) for the 3-wave mixing processes, leading to:

$$ 
\begin{aligned}
\chi_{\text{eff}}^{\text{SFG}} (\omega_3 = \omega_1 + \omega_2, \vec{u}) &= \chi^{(2)} (\omega_3 = \omega_1 + \omega_2). F^{(2)} (\omega_3 = \omega_1 + \omega_2, \vec{u}) \\
&= \chi_{\text{eff}}^{\text{DFG}} (\omega_1 = \omega_3 - \omega_2, \vec{u}) = \chi^{(2)} (\omega_1 = \omega_3 - \omega_2). F^{(2)} (\omega_1 = \omega_3 - \omega_2, \vec{u}) \\
&= \chi_{\text{eff}}^{\text{DFG}} (\omega_2 = \omega_3 - \omega_1, \vec{u}) = \chi^{(2)} (\omega_2 = \omega_3 - \omega_1). F^{(2)} (\omega_2 = \omega_3 - \omega_1, \vec{u})
\end{aligned} 
\quad (55) 
$$
Similar relations can be established for the 4-wave processes.

As seen previously from Equations (44–50), there are two possible eigenmodes of polarization for each wave, which gives a total of $2^3 = 8$ possible configurations of polarization in any direction of propagation for SFG or DFG 3-wave mixing processes: (eee), (ooo), (eo), (oee), (eoe), (eeo), (ooe) and (oeo) in uniaxial crystals or in the principal planes of biaxial crystals, and the same in terms of (+) and (-) polarization modes out of the principal planes of biaxial crystals. Similarly, there are $2^4 = 16$ possible configurations in the case of a SFG and DFG 4-wave mixing or FWM: (oooo), (eeee), (eooo), (oeeee), etc. There are four kinds of symmetry directly depending on these configurations of polarization that are built from ordinary (o) and extraordinary (e) eigenmodes; they are due to: i) cancellations of Cartesian coordinates of the unit electric field, ii) orthogonality between ordinary and extraordinary polarizations, iii) independence of polarization directions occurring at different angular frequencies, and iv) equality between angular frequencies.[6] These symmetries will determine the shape of the (3 × 9) or (3 × 27) matrices describing the field tensor of 3-wave or 4-wave mixing processes, respectively. Below we demonstrate the rationale that underlies the construction of the field tensor matrix based on a few examples applying as well to principal planes of positive uniaxial or biaxial crystals, in particular with a (eo) configuration.

i. According to Equations (44–47) and to Figure 7, the ordinary polarization always lies in the xy-plane, so that $F_{yz}^{eoo} = F_{zjk}^{eoo} = F_{zz}^{eoo} = 0$, which leads to the cancellation of 15 field tensor components in this case. This symmetry depends also on the direction of propagation. By keeping the same example of a propagation in the xy-plane, Equation (45) and Figure 7 show that the extraordinary wave is polarized along the z-axis, then: $F_{xjk}^{eoo} = F_{yjk}^{eoo} = 0$ leading to the annihilation of 8 more field tensor coefficients. This type of relations never applies out of the principal planes of biaxial crystals since all the Cartesian coordinates of the unit electric fields are non-zero.

ii. Ordinary and extraordinary polarizations are always perpendicular so that for example $F_{xxk}^{eoo} + F_{yyk}^{eoo} + F_{zzk}^{eoo} = F_{xjx}^{eoo} + F_{yjy}^{eoo} + F_{zjz}^{eoo} = 0$, which leads to the following relations in addition to those mentioned in (i): $F_{xxy}^{eoo} = -F_{xxy}^{eoo} = -F_{yyx}^{eoo}$ and $F_{xyy}^{eoo} = -F_{xyy}^{eoo} = -F_{xxy}^{eoo}$. Note that these relations do not apply out of the principal planes of biaxial crystals where the two eigenmodes of polarization are never orthogonal.

iii. Equation (44) and Figure 7 show that the ordinary polarization of a uniaxial crystal, or in the xz- and yz-planes of a biaxial crystal, never depends on the frequency, whereas it is the case for a direction of propagation at $\theta = \pi/2$, that is, in the xy-plane, according to Equations (45–47). Again using the example of the configuration (eo), a permutation relation can be found over the two last Cartesian indices, that is, $F_{ijk}^{eoo} = F_{ikj}^{eoo}$, which leads to the equality $F_{xxy}^{eoo} = F_{zyx}^{eoo}$ in addition to the relationships provided by (i) and (ii). Note also that (iii) leads to the full permutation symmetry for the (ooo) field tensor, that is, $F_{ijk}^{eoo} = F_{ikj}^{eoo} = F_{jik}^{eoo} = F_{kjk}^{eoo} = F_{kji}^{eoo}$ for any direction of propagation of a uniaxial crystal or in the xz- and yz-planes of a biaxial crystal.

iv. Equalities between the frequencies of interacting waves of same polarization can provide new permutation rules. It

is the case of SHG ($2\omega = \omega + \omega$), THG ($3\omega = \omega + \omega + \omega$) or the optical Kerr effect, that is, ($\omega = -\omega + \omega + \omega$), which corresponds to a full degenerated FWM. The relations found in (iiii) do not provide new equalities between the field tensor coefficients in the cases of the previous examples with (eo) and (oo) due to (iii). But it is not always the case, in particularly for a propagation out of the principal planes of biaxial crystals where (iiii) is the only allowed permutation. For example it comes $F_{ijk}^{++} (2\omega = \omega + \omega, \bar{u}) = F_{ikj}^{--} (2\omega = \omega + \omega, \bar{u})$ and similarly $F_{ijkl}^{++} (3\omega = \omega + \omega + \omega, \bar{u}) = F_{iljk}^{+-} (3\omega = \omega + \omega + \omega, \bar{u})$.

The symmetries of both the field tensor and the electric susceptibility tensor can lead to situations where their contraction is nil, that is:

$$ \chi_{eff} = \chi^{(n-1)}.F^{(n-1)} = 0 \qquad (56) $$

with $n = (3, 4)$.

This orthogonality between the field tensor and the electric susceptibility tensor occurs when second-order or third-order processes are achieved in specific crystal classes for given configurations of the input polarizations. It happens in the following situations involving uniaxial crystals only: $D_4$ (422) and $D_6$ (622) for the configurations of polarization (eo), (oeo) and (ooe); $C_{4v}$ (4 mm) and $C_{6v}$ (6 mm) for (oee), (eoe) and (eeo); $D_6$ (622), $D_{6h}$ (6mmm), $D_{3h}$ (62m) and $C_{6v}$ (6 mm) for (eo00), (oe00), (ooe0), (oo0e) and also (oeee), (eoee), (eeoe), (eeee).[6,23,24] Note that the corresponding conversion efficiency is nil even if phase-matching or quasi-phase-matching is achieved. This orthogonality criterion is then of prime importance for choosing the right crystal and the well-suited configuration of polarization.

### 5.3. Determination of the Absolute Magnitudes and Relative Signs of the Nonlinear Tensor Coefficients

The field tensor formalism is of particular interest for the characterization of the second- and third-order electric susceptibility tensors of crystals, the goal being to determine both the magnitude and the relative signs of the non-zero coefficients of these tensors from the measurement of frequency conversion efficiencies.

The best situation is to perform such a measurement when the nonlinear process is phase-matched, which allows the conversion efficiency to be enhanced, providing a good signal to noise ratio. Two examples are considered here, allowing us to show two possible strategies.

The first example is that of the famous nonlinear optical crystal, KTiOPO₄ (KTP) that is used nowadays for numerous applications in the visible and far infrared. It crystallizes in the orthorhombic crystal class $C_{2v}$ (mm2), it belongs to the biaxial optical class, and its optical sign is positive.[77]

There are 7 non-zero and independent coefficients of the second-order electric susceptibility tensor for this crystal class.[20]

Kleinman's assumption can be used when both the absorption and the wavelength dispersion of the principal refractive indices are weak, which is often relevant in the used wavelength ranges. In this case, the second-order electric susceptibility
<figure><img src="image_9.png" /><figcaption>Figure 8. Field tensor coefficients calculated as function of the phase-matching angle for Type II SHG ($\lambda_\omega = 1.064$ µm; $\lambda_{2\omega} = 0.532$ µm) in KTP. Adapted with permission.[78] Copyright 2024 Optica Publishing Group.</figcaption></figure>

tensor becomes symmetric over the three Cartesian indices, which gives for $C_{2\omega}$: $\chi_{xxz} = \chi_{xzx} = \chi_{zxx} \neq \chi_{yyz} = \chi_{yz\gamma} = \chi_{zyy} \neq \chi_{zzz}$, corresponding respectively to $\chi_{15} = \chi_{31} \neq \chi_{24} = \chi_{32} \neq \chi_{33}$ using the contracted notation. This leads to three independent coefficients to determine.

It is interesting at this stage to consider the corresponding effective coefficients expressed on the basis of the field tensor. For example, a phase-matched Type II SHG configuration of polarization in the positive biaxial crystal KTP, leads to:

$$ \chi_{eff} (2\omega = \omega + \omega, \vec{u}) = \chi_{15} F_{15} + \chi_{24} F_{24} + \chi_{31} F_{31} + \chi_{32} F_{32} + \chi_{33} F_{33} \quad (57) $$

with

$$ \left\{ \begin{array}{l} F_{15} = \\ \left[ e_x^- (2\omega, \vec{u}) \cdot e_x^+ (\omega, \vec{u}) \cdot e_z^- (\omega, \vec{u}) + e_x^- (2\omega, \vec{u}) \cdot e_z^+ (\omega, \vec{u}) \cdot e_x^- (\omega, \vec{u}) \right] \\ F_{24} = \\ \left[ e_y^- (2\omega, \vec{u}) \cdot e_y^+ (\omega, \vec{u}) \cdot e_z^- (\omega, \vec{u}) + e_y^- (2\omega, \vec{u}) \cdot e_z^+ (\omega, \vec{u}) \cdot e_y^- (\omega, \vec{u}) \right] \\ F_{31} = e_z^- (2\omega, \vec{u}) \cdot e_x^+ (\omega, \vec{u}) \cdot e_x^- (\omega, \vec{u}) \\ F_{32} = e_z^- (2\omega, \vec{u}) \cdot e_y^+ (\omega, \vec{u}) \cdot e_y^- (\omega, \vec{u}) \\ F_{33} = e_z^- (2\omega, \vec{u}) \cdot e_z^+ (\omega, \vec{u}) \cdot e_z^- (\omega, \vec{u}) \end{array} \right. \quad (58) $$

The Cartesian coordinates $e_x^{+/-}$ and $e_x^{0/}$ stand for a propagation out-of and in the principal planes, respectively. They can be calculated using Equations (44)-(49) according to the direction of propagation. The Type II SHG phase-matching cone lying from the xy-plane to the yz-plane,[78] the corresponding variation of the field tensor coefficients are plotted as a function of the corresponding phase-matching angle $\varphi$ in Figure 8 for ($\lambda_\omega = 1.064$ µm; $\lambda_{2\omega} = 0.532$ µm).

Figure 8 shows first that $\chi_{15}(0.532 \mu m)$ is the only excited nonlinear coefficient in the yz-plane, that is, at $\varphi = 90^\circ$. The measure-

<figure><img src="image_10.png" /><figcaption>Figure 9. Type II SHG ($\lambda_\omega = 1.064$ µm; $\lambda_{2\omega} = 0.532$ µm) in KTP. Calculated and measured second-harmonic power as function of the phase-matching angle; the coefficients $d_{15}$ and $d_{24}$ correspond to $\chi_{15}/2$ and $\chi_{24}/2$, respectively. Adapted with permission.[72] Copyright 2024 Optica Publishing Group.</figcaption></figure>

ment of the corresponding Type II SHG conversion efficiency, which depends on the square of the effective coefficient, can directly lead to the magnitude of $\chi_{15}(0.532 \mu m)$: it had been found that $|\chi_{15}(0.532 \mu m)| = 2.8$ pm/V.[78] From this data combined with such a measurement performed in the xy-plane at $\varphi = 23^\circ$, Figure 8 shows that it is then possible to access to two possible values for the magnitude of $\chi_{24}(0.532 \mu m)$ according to its relative sign with $\chi_{15}(0.532 \mu m)$: $|\chi_{24}(0.532 \mu m)| = 5.3$ pm/V in the hypothesis of same sign, and $|\chi_{24}(0.532 \mu m)| = 6.5$ pm/V in the hypothesis of an opposite sign. The choice between these two values can be then done by measuring the Type II SHG conversion efficiency for different directions of propagation ranging between the xy- and yz-planes. Actually, the weights of excitation of the two nonlinear coefficients, that is, $F_{15}$ and $F_{24}$, vary antagonistically while the three other field tensor coefficients $F_{31}$, $F_{32}$ and $F_{33}$ are negligible, as shown in Figure 9.

Figure 9 well indicates a clear difference between the two sign assumptions, with or without a zero crossing according to the phase-matching angle. The corresponding measurements had been carried out using a KTP crystal cut as a sphere. It undoubtedly showed an agreement between the experimental data and the calculated curve that never vanishes to zero, which means an identical sign for $\chi_{15}$ and $\chi_{24}$. It was then possible to determine the suited magnitude of $\chi_{24}$, that is, $|\chi_{24}(0.532 \mu m)| = 5.3$ pm/V. And as mentioned above, under Kleinman's assumption, the magnitudes of $\chi_{31}$ and $\chi_{32}$ are the same than those of $\chi_{15}$ and $\chi_{24}$, respectively. Finally, using Type I SHG and the sphere experimental setup, it had also been possible to determine $\chi_{33}$, which is not detailed here.[78]

The second example is that of the positive biaxial crystal BaGa$_4$Se$_7$ (BGSe) belonging to the monoclinic crystal class $C_s$ ($m$). This crystal is a promising material for mid-infrared
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Type I SHG as function of the phase-matching angle in the xz-plane of BGSe. a) Calculated tuning curve of the fundamental wavelength. b) Calculated and measured square of the effective coefficient. The coefficients $d_{21}$ and $d_{23}$ correspond to $\chi_{21}/2$ and $\chi_{23}/2$, respectively. Adapted with permission.[80] Copyright 2024 Optica Publishing Group.</figcaption></figure>

generation.[79] However, the characterization of the coefficients of its second-order electric susceptibility tensor is much more complicated than in the case of KTP. Indeed, BGSe has six nonzero and independent components under Kleinman's assumption, that is, using the contracted notation and the proper dielectric frame: $\chi_{15} = \chi_{31}$, $\chi_{16} = \chi_{21}$, $\chi_{22}$, $\chi_{23} = \chi_{34}$, $\chi_{24} = \chi_{32}$ and $\chi_{33}$. The field tensor analysis shows that it is necessary to combine non-phase-matched and phase-matched measurements on the one hand, and to consider the wavelength evolution of the phase-matching angle in different principal planes in order to access to the different relative signs on the other hand.[80] The case of the relative signs between $\chi_{21}$ and $\chi_{23}$ is detailed hereafter. Type I SHG phase-matched in the xz-plane is suited for that purpose because the effective coefficient only depends on these two coefficients whose amplitude had been previously determined, that is: $|\chi_{21}(0.532 \mu m)| = 10.6 \text{ pm/V.}^{[80]}$ and $|\chi_{23}(0.532 \mu m)| = 28.4 \text{ pm/V.}^{[81]}$ Since a principal plane is considered, the field tensor coefficients can be expressed in terms of ordinary and extraordinary waves, which gives in the xz-plane at $V < \theta < \pi/2$:

$$ \chi_{\text{eff}}(2\omega = \omega + \omega, \vec{u}) = \chi_{21}(2\omega = \omega + \omega) \cdot F_{21}(2\omega = \omega + \omega, \vec{u}) + \chi_{23}(2\omega = \omega + \omega) \cdot F_{23}(2\omega = \omega + \omega, \vec{u}) \quad (59) $$

with

$$ \left\{ \begin{array}{l} F_{21}(2\omega = \omega + \omega, \vec{u}) = e_y^o(2\omega, \vec{u}) e_x^e(\omega, \vec{u}) e_x^e(\omega, \vec{u}) \\ F_{23}(2\omega = \omega + \omega, \vec{u}) = e_y^o(2\omega, \vec{u}) e_z^e(\omega, \vec{u}) e_z^e(\omega, \vec{u}) \end{array} \right. \quad (60) $$

where the expressions of $e_y^o$ and $e_x^e$ can easily be found using Equations (44) and (45) at $\varphi = 0$.

The corresponding calculations as a function of wavelength shown in Figure 10 are performed in the two hypotheses of relative

sign using the numerical amplitudes of $\chi_{21}$ and $\chi_{23}$ given above and also Miller's rule.[82]

The corresponding experimental data given in Figure 10 undoubtedly showed an agreement with the curve that starts from zero, which means different signs for $\chi_{21}$ and $\chi_{23}$.[80]

## 6. Beyond the Conventional Susceptibility and Field Expansion: Nonlinear Response as a Multi-Variable Time-Space Convolution of Material and Field Tensors

Modeling of nonlinear optical phenomena under weak field condition is sustained by of nonlinear electric susceptibility tensor $\chi^{(n)}(\omega_1, \omega_2 \dots \omega_n)$ with coefficients that take the familiar and ubiquitous expressions $\chi_{i_1 i_2 \dots i_n}(\omega_1, \omega_2 \dots \omega_n)$. These expressions stand for the tensor coefficients of the electric susceptibility of rank $n$, order $n-1$, at the angular frequencies $(\omega_1, \omega_2 \dots \omega_n)$.[83]

Knowledge of these expressions provide *inter alia* the Fourier spectrum of the $(n-1)^{\text{th}}$ order response of the medium under an input with space-time Fourier spectrum $E_{ij}(\vec{k}_j, \omega_j)$ for $1 \le j \le n$, this expression standing for the $(\vec{k}_j, \omega_j)$ Fourier components of a vectorial input's spatial component along an axis $i_j = (x, y, z)$. The wide use of the expression

$$ \begin{array}{l} P_{(n)}^{i_1}(\vec{k}_1; \omega_1) = \chi_{(n)}^{i_1 i_2 \dots i_n}(\vec{k}_1, \vec{k}_2 \dots \vec{k}_n; \omega_1, \omega_2 \dots \omega_n) \\ E_{i_2}(\vec{k}_2, \omega_2) E_{i_3}(\vec{k}_3, \omega_3) \dots E_{i_n}(\vec{k}_n, \omega_n) \end{array} \quad (61) $$

as a starting point tends to be oblivious of both its more physically intuitive and mathematically grounded upstream roots as well as of their implications. It is preferable to start from our usual perception of space and time and not from more elaborate time and/or space frequencies that are the less intuitive Fourier transform of the former. One should therefore try and express nonlinear phenomena by way of some upstream expression of the response function where only time and space variable are taking part, whereas the above expression appear as downstream consequences by way of an adapted multiple time-space Fourier transformation. This point of view is particularly important toward the purpose of setting firm foundations to the field tensor such as would pertain for example to exposure to multiple time-varying field in nonlinear media mediated by non-local delaying interactions.

For the stake of simplicity, we shall start here from a *purely time-dependent* input-output pair, the spatial dependence disappearing in the case of a purely local phenomenon as it will be further detailed. It will then be straightforward to extend this purely time dependent approach to the more general combined time-space one. Moreover, for the sake of simplification and with no prejudice to the bottom-line issues, we shall first restrict ourselves to a scalar isotropic medium and further generalize to an anisotropic one where tensor inputs and outputs are to be considered.

It is generally admitted in physics textbooks, such as in nonlinear optics, that an arbitrary nonlinear response in time can be expanded into a sum of multilinear responses. While being taken
as granted for a starting point, this assumption requires justification based on reasonably rigorous mathematical foundations.

Let us consider a couple of input and output, or action and response time dependent functions $E(t)$, $P(t)$, it is indeed generally admitted that the output can be expressed as a sum of multiple integrals of the form:

$$P[t, E(.)] = \sum_{n=1}^{\infty} P^n[t, E(.)] \\ = \iint \dots \int \tilde{\mathbb{R}}^n(t; \tau_1, \tau_2, \dots, \tau_n) E(\tau_1) E(\tau_2) \dots E(\tau_n) d\tau_1 \tau_2 \dots \tau_n \qquad (62)$$

The response $P$ at time $t$ is a scalar (to be further generalized to a vector or a tensor) which depends on the input $E$ as a whole, that is at all times and not just at any given time, hence the notations $E(.)$ and $P[t, E(.)]$ that have been chosen here.

We propose to call this important property as the basic theorem for a weak nonlinear response. We will later clarify the label “weak” by way of the convergence condition imposed on the development of the nonlinear response function $P[t, E(.)]$.

The special cases of linear (e.g., $n = 1$) and quadratic (e.g., $n = 2$) contributions of the induced polarization of matter, the general response takes the more familiar forms:

$$P^1[t, E(.)] = \int \tilde{\mathbb{R}}^1(t; \tau) E(\tau) d\tau \qquad (63)$$

$$P^2[t, E(.)] = \iint \tilde{\mathbb{R}}^2(t; \tau_1, \tau_2) E(\tau_1) E(\tau_2) d\tau_1 d\tau_2 \qquad (64)$$

We justify in **Appendix A** those expressions, before proceeding to their generalization and subsequent multivariable Fourier transform, to serve as a sound basis for the definition of frequency dependent linear and nonlinear electric susceptibility tensors.

Based now on Appendix A with the basic theorem for a scalar isotropic purely time-dependent response of a medium, we are now in a position to generalize to the more general case of an anisotropic space-time-dependent response which is at this stage a matter of technicality entailing no difficulty.

We now use in the following vectors notations for the space variable $\vec{x}$, the induced polarization of matter $\vec{P}$, the electric field $\vec{E}$ and the nonlinear functional $\vec{\tilde{F}}$:

$$\vec{P}(t; \vec{x}) = \vec{\tilde{F}}[\vec{E}(t, \vec{x})] \qquad (65)$$

The vector functional $\vec{\tilde{F}}$ relates an input electric field to an output response which stands for the induced polarization in nonlinear optics. This notation features a $(t; \vec{x})$ dependence on both sides may be somewhat misleading. It is left here simply to remind that the left and right hand-side vector functions are $(t; \vec{x})$ dependant while in no way meant to single-out their values at $(t; \vec{x})$.

The vector functional $\vec{\tilde{F}}$ is assumed to be “well behaved” in the sense of fulfilling conditions for further derivations already

detailed in Appendix A as well as nonlinear in the most general sense that is:

$$\vec{P}(t; \vec{x}) = \vec{\tilde{F}}[\vec{E}_1(t, \vec{x}) + \vec{E}_2(t, \vec{x})] \neq \vec{\tilde{F}}[\vec{E}_1(t, \vec{x})] + \vec{\tilde{F}}[\vec{E}_2(t, \vec{x})] \qquad (66)$$

As for vectors of electric fields, we assume that they fulfil conditions met for well-behaved distributions in the sense of generalized functions encompassing such distributions as that of Heaviside or Dirac as extreme cases.[84] They can be further Fourier transformed in the sense again of well-behaved functions that can be considered as regular distributions in the general case.

At this stage, the notion of nonlinearity does not entail any given order such as attached to individual nonlinear optical processes, namely SHG quadratic, FWM cubic or other as it will appear later in this paragraph.

A more general functional can be expressed in terms of **n** input-fields, such as in:

$$\vec{P}(t; \vec{x}) = \vec{\tilde{F}}[E_1(t, \vec{x}), \dots, E_n(t, \vec{x})] \qquad (67)$$

By way of a technical extension of the derivation in **Appendix A**, the more general nonlinear vector functional $\vec{\tilde{F}}$ can be expanded in terms of a sum of multiple convolution vector functions of increasing order $p$ provided perturbation conditions are satisfied. This basically amounts to assume that all participating electric fields are much smaller than the atomic electric fields within the material of interest.

Under this condition, each $p$-order time-space response can be expressed as a multiple convolution of order $p$ with the full tensor contraction (hereafter as $\text{cont}(\tilde{\mathbb{F}}^{(p)}(.), \tilde{\mathbb{F}}^{(p)}(.))$) over the $p+1$ Cartesian indices of the time-space field tensor with the multiple impulse function expressed as follows:

$$\vec{P}(t; \vec{x}) \cdot \vec{u} = \int d\tau_1 \dots \int d\tau_p \int d\vec{\rho}_1 \dots \int d\vec{\rho}_p \text{cont}[(\tilde{\mathbb{F}}^{(p)}(\tau_1 \dots \tau_p; \vec{\rho}_1 \dots \vec{\rho}_p) \dots \tilde{\mathbb{F}}^{(p)}(t - \tau_1 \dots t - \tau_p; \vec{x} - \vec{\rho}_1 \dots \vec{x} - \vec{\rho}_p)] \qquad (68)$$

with the read-field tensor expressed in space-time variables:

$$\tilde{\mathbb{F}}^{(p)}(t_1 \dots t_p; \vec{x}_1 \dots \vec{x}_p) = \vec{u} \otimes \vec{E}(t_1, \vec{x}_1) \otimes \dots \otimes \vec{E}(t_p, \vec{x}_p) \qquad (69)$$

and the impulse response $\tilde{\mathbb{F}}^{(p)} = \tilde{\chi}^{(p)}(\tau_1 \dots \tau_p; \vec{\rho}_1 \dots \vec{\rho}_p)$

The different time-space variables are now unambiguously designating the application points of the related functions as would be demanded for nonlocal interactions such as for nonlinear plasmonic propagative media. In accordance with classical conventions, we will henceforth switch to the $\tilde{\chi}^{(p)} = \tilde{\mathbb{F}}^{(p)}$ notation for a nonlinear electric susceptibility of order $p$ (expressed here in time-space variables leaving the final $\chi^{(p)}$ notation to represent that same electric susceptibility now expressed with respect to $(\omega, \vec{k})$ frequency variables following a multi-variable Fourier transform.

The following equation

$$\vec{P}(t; \vec{x}) \cdot \vec{u} = \int d\tau_1 \dots \int d\tau_p \int d\vec{\rho}_1 \dots \int d\vec{\rho}_p \dots \\ \tilde{\chi}^{(p)}(t - \tau_1 \dots t - \tau_p; \vec{x} - \vec{\rho}_1 \dots \vec{x} - \vec{\rho}_p) \vec{u} \otimes \vec{E}(\tau_1, \vec{\rho}_1) \otimes \dots \otimes \vec{E}(\tau_p, \vec{\rho}_p) \qquad (70)$$
can be alternatively written, due to the fundamental property of convolutions as:

And corresponds more specifically to a Dirac tensor field polarized along $u_i = \{x, y \text{ or } z\}$ for $1 \le i \le p$ and of unit amplitudes. The fully index expression for the multiple-convolution writes:

$$ (\vec{P})_{u} (t, \vec{r}) = \left\{ \sum_{w, u_1, \dots, u_p} \int d\tau_1 \int d\vec{\rho}_1, \dots \int d\tau_p \int d\vec{\rho}_p \dots \right. \\ \left. \dots (\vec{u})_{w} [\vec{\chi}_{w, u_1, \dots, u_p}^{(p)} (\tau_1, \vec{r}_1, \dots, \tau_p, \vec{r}_p) \dots E_{u_1}^{(p)} (t - \tau_1, \vec{r} - \vec{\rho}_1) \dots E_{u_p} (t - \tau_p, \vec{r} - \vec{\rho}_p) \dots \right\} \tag{76} $$

$$ \begin{gathered} \vec{P}(t; \vec{x}) \cdot \vec{u} = \int d\tau_1 \dots \int d\tau_p \int d\vec{\rho}_1 \dots \int d\vec{\rho}_p \dots \\ \vec{\chi}^{(p)}(\tau_1 \dots \tau_p; \vec{\rho}_1 \dots \vec{\rho}_p) \vec{u} \otimes \vec{E}(t - \tau_1, \vec{x} - \vec{\rho}_1) \otimes \dots \otimes \vec{E}(t - \tau_p, \vec{x} - \vec{\rho}_p) \end{gathered} \tag{71} $$

The impulse response terminology originates from the special case of a $\vec{F}_\delta$ field-tensor made-up of pure space-time Dirac impulses (e.g., fully instantaneous and localized):

$$ \vec{F}_\delta = \delta^{(p)} (\tau_1 - t_1 \cdots \tau_p - t_p; \vec{\rho}_1 - \vec{r}_1 \cdots \vec{\rho}_p - \vec{r}_p) \tag{72} $$

With all the properties of a multi-variable Dirac distribution where the tensorial indices are omitted for the sake of simplicity, it comes:

$$ \begin{gathered} \delta^{(p)} (\tau_1 - t_1 \cdots \tau_p - t_p; \vec{\rho}_1 - \vec{r}_1 \cdots \vec{\rho}_p - \vec{r}_p) = 0 \quad \text{for} \\ \tau_1 \neq t_1 \text{ or } \cdots \neq t_p, \vec{\rho}_1 \neq \vec{r}_1 \text{ or } \cdots \neq \vec{\rho}_p = \vec{r}_p \end{gathered} \tag{73} $$

$\delta^{(p)}(\tau_1 - t_1 \cdots \tau_p - t_p; \vec{\rho}_1 - \vec{r}_1 \cdots \vec{\rho}_p - \vec{r}_p)$ diverging for $\tau_1 = t_1$ or $\cdots = t_p; \vec{\rho}_1 = \vec{r}_1$ or $\cdots = \vec{\rho}_p = \vec{r}_p$ with the following more rigorous definition of a Dirac distribution:

$$ \begin{gathered} \int dt_1 d\vec{r}_1 \cdots dt_p d\vec{r}_p \delta^{(p)} (\tau_1 - t_1 \tau_p - t_p; \vec{\rho}_1 - \vec{r}_1 \cdots \vec{\rho}_p - \vec{r}_p) \\ \varphi (t_1 \cdots t_p; \vec{r}_1 \cdots \vec{r}_p) = \cdots \int dt_1 d\vec{r}_1 \cdots dt_p d\vec{r}_p \delta^{(p)} (t_1 \cdots t_p; \vec{r}_1 \cdots \vec{r}_p) \\ \varphi (\tau_1 - t_1 \cdots \tau_p - t_p; \vec{\rho}_1 - \vec{r}_1 \cdots \vec{\rho}_p - \vec{r}_p) = \varphi (\tau_1 \cdots \tau_p; \vec{\rho}_1 \cdots \vec{\rho}_p) \end{gathered} \tag{74} $$

valid for any infinitely differentiable $\varphi(t_1 \dots t_p; \vec{r}_1, \dots \vec{r}_p)$ function of $4p$ variables (e.g., $p$ in time and $3p$ in space).

In the present context, taking $\varphi(t_1 \dots t_p; \vec{r}_1 \dots \vec{r}_p) = \vec{\chi}^{(p)}(t_1 \dots t_p; \vec{r}_1 \dots \vec{r}_p)$ leads to identify $\vec{\chi}^{(p)}(t_1 \dots t_p; \vec{r}_1 \dots \vec{r}_p)$ with the response of the NL system of order $p$ to the multiple Dirac function input as a set of $p$ time-space vector field distributions of instantaneous durations at times $t_i$ and point-wise localization at positions $\vec{r}_i$ for $1 \le i \le p$, hence the terminology of impulse response.

The full tensorial expression of the multiple Dirac function is self-understandably as follows:

$$ (\vec{u}) w \delta_{u_1, \dots, u_p}^{(p)} (t_1 \cdots t_p; \vec{r}_1 \cdots \vec{r}_p) \tag{75} $$

Its multiple Fourier transform leads to the usual expressions including time-space frequency dispersion whereby a convolution in space time simplifies into a product in the frequency space:

$$ \begin{gathered} (\vec{P})_u (\omega, \vec{k}) = \sum_{w, u_1, \dots, u_p} (\vec{u}) w \chi_{w, u_1, \dots, u_p}^{(p)} (\omega_1, \vec{k}_1, \dots, \omega_p^r, \vec{k}_p) \cdots \\ E_{u_1} (\omega_1, \vec{k}_1) \cdots E_{u_p} (\omega_p, \vec{k}_p) \end{gathered} \tag{77} $$

The Fourier transform of the electric field tensor is simply the tensor product of the individual fields, namely:

$$ E_{u_1} (\omega_1, \vec{k}_1) = \mathcal{F} [E_{u_1} (t, \vec{r})] = \iint E_{u_1} (t, \vec{r}) e^{i(\omega_1 t + \vec{r} \cdot \vec{k}_1)} dt d^3 \vec{r} \tag{78} $$

The multiple Fourier transform leading to the usual time and space frequency dispersed expression of the electric susceptibility is as follows:

$$ \begin{align} \chi_{w, u_1, \dots, u_p}^{(p)} (\omega_1, \vec{k}_1, \dots, \omega_p^r, \vec{k}_p) &= \mathcal{F}[\vec{\chi}_{w, u_1, \dots, u_p}^{(p)} (\tau_1, \vec{r}_1, \dots, \tau_p, \vec{r}_p)] \nonumber \\ &= \cdots \chi_{w, u_1, \dots, u_p}^{(p)} (\omega_1, \vec{k}_1, \dots, \omega_p^r, \vec{k}_p) = \mathcal{F}[\vec{\chi}_{w, u_1, \dots, u_p}^{(p)} (\tau_1, \vec{r}_1, \dots, \tau_p, \vec{r}_p)] \nonumber \\ &= \cdots \chi_{w, u_1, \dots, u_p}^{(p)} (\omega_1, \vec{k}_1, \dots, \omega_p^r, \vec{k}_p) = \mathcal{F}[\vec{\chi}_{w, u_1, \dots, u_p}^{(p)} (\tau_1, \vec{r}_1, \dots, \tau_p, \vec{r}_p)] \nonumber \\ &= \cdots = \int d\tau_1 \iiint d^3 \vec{r}_1 \cdots \int d\tau_p \iiint d^3 \vec{r}_p \vec{\chi}_{w, u_1, \dots, u_p}^{(p)} (\tau_1, \vec{r}_1, \dots, \tau_p, \vec{r}_p) e^{i(\omega_1 \tau_1 + \cdots + \omega_p \tau_p + \vec{k}_1 \cdot \vec{r}_1 + \cdots + \vec{k}_p \cdot \vec{r}_p)} \tag{79} \end{align} $$

A $\sqrt{2\pi}^{-4p}$ multiplying factor can be introduced to symmetrize the direct and inverse Fourier transform expressions.

Note that all integrations over the $\tau_i$ delays have to be performed over positive values corresponding to the past, namely $t - \tau_i \ge 0$. In order to allow for the implementation of classical rules of Fourier transform. Therefore, for the sake of simplifying notations, all $\tau_i$- Fourier transforms are made to implicitly containing a $Y(\tau_i)$ Heaviside distribution equal to unity for $\tau_i$ positive and cancels otherwise, namely by introduction of

$$ \int Y(\tau_1) d\tau_1 \cdots \int Y(\tau_p) d\tau_p \tag{80} $$
in the above expression of the frequency dependent electric susceptibility as well as electric fields.

Note also that the multivariable transform of the field-tensor is the product of the Fourier transform of individual electric fields, namely:

$$ 
\begin{aligned}
F(\omega_1 \cdots \omega_p; \vec{k}_1 \cdots \vec{k}_p) &= \mathcal{F}[E_{u_1}(t_1, \vec{r}_1)] \otimes \cdots \otimes \mathcal{F}[E_{u_p}(t_p, \vec{r}_p)] \\
&= \bar{u} \otimes \bar{E}(\omega_1, k_1) \otimes \cdots \otimes \bar{E}(\omega_p, k_p) \quad (81)
\end{aligned}
$$

One will have recognized here the time space expression of the field tensor, namely:

$$ F(t_1 \cdots t_p, \vec{r}_1 \cdots \vec{r}_p) = \mathcal{F}[E_{u_1}(t_1, \vec{r}_1)] \otimes \cdots \otimes \mathcal{F}[E_{u_p}(t_p, \vec{r}_p)] \quad (82) $$

At this stage, we can safely express the nonlinear response of order n in a proper context and discuss its physical meaning and practical implications.

A number of configurations involve a sequence of two incoming optical perturbations with different time and space features, leading to an outgoing beam as a result of nonlinear interactions in a given medium. Such a situation is met in the all-pervading pump and probe configuration toward time-resolved experiments or write and read fields in the different context of optical poling. These configurations require the extension of the above framework from a single electric field $\vec{E}(\vec{r}, t)$ to a couple of distinct field perturbations that we conventionally designate respectively as $\vec{E}^w(\vec{r}, t)$ and $\vec{E}^r(\vec{r}, t)$ each one featuring its own time-space distribution. We explore the extension of the previous single-field frame in Appendix B which provides full expressions of the extended multilinear operators and the associated impulse response that require independent integration over the different time and space variables related to the two interacting fields.

We show that the combined write-read process can be expressed by a single multiple convolution response:

$$ \vec{P}(t, \vec{r}) = \tilde{\mathcal{L}}_{(p,q)}^{(r,w)} \circ \tilde{\mathcal{F}}_{(p,q)}^{(r,w)} \quad (83) $$

which expresses via the ⊙ symbol connecting the generalized multilinear response function $\tilde{\mathcal{L}}_{(p,q)}^{(r,w)}$ and the $\tilde{\mathcal{F}}_{(p,q)}^{(r,w)}$ field tensor via the combination of a multiple convolution over $4(p+q)$ spacetime set of variables together with a full tensor contraction over Cartesian indices. The $(p+q)$ order nonlinear effect involves here p-times the write field and q-times the read field.

## 7. Conclusion

In this review, we have explored the practical as well as conceptual relevance of the field tensor throughout a variety of coherent and non-coherent processes. To recall but a few, we have shown that the field tensor comprehends the full extension of polarization and phase properties of the intervening fields. Indeed, the field tensor has been shown to uniquely comprise all the optical physics ingredients that are required to analyze and discuss in a synthetic way a broad variety of nonlinear phenomena. However, as expected, the field tensor was never found useful on its

own, but only by way of an intimate one-to-one conjunction with the corresponding material based (non-)linear electric susceptibility tensor of the same rank. Throughout this review, we have made use of the all-pervading algebraic operation underlying the field-matter tensor coupling known as tensor contraction which is a generalization of the usual scalar product of vectors. Its properties are particularly striking and useful when expressed in the irreducible tensor formalism. We have basically explored two directions with different results for each.

Whenever symmetry considerations predominate and conversely propagation is less important, such as in thin films or nanoscale samples, the key issue is to express the field tensor and the electric susceptibility tensor in the same tensor reference frame. This can of course be readily done in Cartesian coordinates with a subsequent reduction taking into considerations symmetry. However, it results in complex expressions that scramble components and fail to properly reflect the core invariances that are attached to different symmetry types. By contrast, we have shown that the irreducible tensor formalism allows to decouple different rotational orders corresponding to combinations of circular polarization states with different parities. These irreducible tensor components make-up for a symmetry eigenstate of the field tensor very much like adequate o/e and ± combination account for the propagation tensor eigenstates in the context of crystalline nonlinear optics. The power of the irreducible tensor formalism is showing up in the designation of specific combination of circular polarization states that are capable of singling-out specific irreducible components of the nonlinear electric susceptibility of a material and its environment. Indeed, a one-to-one coupling between irreducible components of the same J order for both the electric susceptibility and field tensors is evidenced both theoretically and experimentally, the latter in a nano-plasmonic context. The case of statistically oriented media with built-in features or an externally applied order by means of electrical and/or optical fields has been reviewed and shown to lead to the definition of two field tensors. The first one is a “write” field tensor which is meant to impart a given rotational order to a plastic medium, and the second one a “read-out” or “probing” field tensor, both made of adequate combinations of circular polarization states. It is possible to imprint any combination of irreducible orders onto a photo-polymerizable material, with so-far unexploited consequences in terms of high capacity tensorial optical data storage with protected access due to the need of complex nonlinear combinations for the read-out.

Extensions of this formalism to such processes as coherent-anti-Stokes Raman spectroscopy, harmonic light scattering, high-harmonic generation, transfer or photon angular momentum to fast-rotating molecules and others have been addressed, all of them pointing-out at the “filtering” potential of the field tensor in order to characterize the nonlinear properties, including intermediate vibrational resonances, of a given medium.

Second- and third-order phase-matched nonlinear propagation phenomena are also discussed, making central the use of the field tensor by way of combinations of the polarization eigenstates of the propagation equations in both uniaxial and biaxial crystals. This allows to discuss in the proper frame the emergence of birefringence-based phase-matching and the existence as well as of systematic expressions of the corresponding effective
nonlinear coefficients leading in particular to the determination of the absolute values and relative signs of the coefficients of the nonlinear electric susceptibility tensors. This is developed in the prototypical case of some available inorganic nonlinear crystals of broad application interest, in particular toward three- and four-wave parametric emission that are currently playing a major role in the production of entangled photons and related perspectives in future quantum information and computation schemes.

Finally, we revisit and justify the expression of the nonlinear response in the perturbation regime limit in terms of a general time-space multi-convolutional development that had been so far taken as a starting basis with no rigorous justification. Based on this firm ground, usual expressions of the susceptibilities and of the tensor field on a time-space Fourier series can be derived allowing to retrieve the usual form used in the literature. We furthermore extend the previous considerations to the case of a sequence of two perturbating fields, such as in pump-prove or write-read configurations that are respectively attached to time-resolved nonlinear spectroscopy and to light induced structural modifications of photosensitive materials.

Some more advanced issues have been purposely left aside, due to their nature of an ongoing research or of their more speculative nature. Among the current works where the formalism associated to the field tensor is most probably of high importance is quantum nonlinear optics, aiming at the generation and frequency conversion of quantum states of light. It benefits today from considerable technological progresses where efficient interactions between photons in nonlinear media are made possible, with applications in quantum computing or sensitive metrology. It however requests to introduce quantum expressions of the field tensor, in which non-classical states of light are involved through tensorial forms for the efficient production of entangled pairs or triplets of photons. Another interest for the use of field tensors is the extension to situations involving chirality, either from the optical side through spin or orbital angular momentum, or in chiral molecular or artificial media. In optical chirality, fields carry orbital angular momentum as a spatial property, such as occurring from a phase vortex, extending the field tensor along propagation pathways based on tensorial extensions of the nonlinear Maxwell-Helmholtz equations and others. It is also remarkable that singularities appearing in the field tensor under specific situations, such as non-paraxial phase/polarization-shaped beams or beams propagating along a crystal optical axis for example, can lead to novel situations, like original phase-matching scenarios based on the Berry geometric phase.[85,86] At last, in chiral molecular media or nanostructures, electric susceptibilities involve additional magnetic dipole or quadrupolar electric interactions that complicate the material susceptibility tensor to which the field tensor must be coupled.[7,87] These components have been evidenced and accounted for in the modelling of nonlinear chirality from surfaces,[88] and in biological molecules such as collagen.[89] Of high interest is also the interaction between complex vortex beams and chiral materials, for which mostly linear optical processes have been studied so far.[90]

Over the last decades, the field of modern optics has been prone to constant renewals that forbid to limit ourselves to the sole extrapolation of current developments. In any scenario, it is likely that the field tensor will remain a useful tool albeit in evolving formats along developments to come.

## Appendix A: Proof of the Expression of the Weak Nonlinear Response as a Sum of Multilinear Responses

We proceed to show that the decomposition into a sum of multilinear kernel[91] functions as per Equation 62 is legitimate. Let us consider N sampling values in the time interval [-A, A] to be further extended, when needed, to the whole-time space spanned by $E(t)$ from $-\infty$ to $\infty$ by allowing $A \to \infty$ under certain divergence conditions. We now divide the [-A, A] interval in N time smaller time intervals $[t_i, t_{i+1}]$ defined by $\Delta\tau = t_{i+1} - t_i = 2A/N$, corresponding to N equidistant sampling times $t_i$ with $1 \le i \le N$ spanning the [-A, A] interval starting from $t_1 = -A$ to $t_N = A$.

Given the N discrete sampling times $t_1, t_2, \dots, t_N$ as defined above, we now define a corresponding basis of rectangular sampling functions[92] centered at times $t_i$, $1 \le i \le N$ given by

$$ \Delta_i(t) = \frac{1}{\Delta\tau} \qquad (A1) $$

that take a constant $1/\Delta\tau$ value over the $[t_i - \Delta\tau/2, t_i + \Delta\tau/2]$ interval and cancel elsewhere. As $\Delta\tau \to 0^{[93]}$ one can easily show that the sampling function $\Delta_i(t)$ tend toward Dirac functions $\delta(t - t_i)$ also centered at $t_i$.

Indeed, for any well-behaved[94] function it is straightforward to prove that

$$ \begin{aligned} \mathrm{Int} = \int_{-\infty}^{\infty} f(t) \Delta_i(t) dt &= \frac{1}{\Delta\tau} \int_{t_i - \Delta\tau/2}^{t_i + \Delta\tau/2} f(t) dt \to f(t_i) = \int_{-\infty}^{\infty} f(t) \delta(t - t_i) dt \end{aligned} \qquad (A2) $$

for $\Delta\tau \to 0$, making use here of the definition of integral Int as the primitive $F(t)$ of its integrand $f(t)$, that is $dF/dt = f$ and then

$$ \begin{aligned} \mathrm{Int} &= \frac{1}{\Delta\tau} \int_{t_i - \frac{\Delta\tau}{2}}^{t_i + \frac{\Delta\tau}{2}} f(t) dt = \frac{1}{\Delta\tau} \left[ F\left(t_i + \frac{\Delta\tau}{2}\right) - F\left(t_i - \frac{\Delta\tau}{2}\right) \right] \\ \Rightarrow \mathrm{Int} &= \frac{1}{\Delta\tau} \left[ \left( \frac{\Delta\tau}{2} \right) \frac{dF}{dt}(t_i) + \left( \frac{\Delta\tau}{2} \right) \frac{dF}{dt}(t_i) \right] + 0(\Delta\tau) \\ \Rightarrow \mathrm{Int} &= \frac{1}{\Delta\tau} \left[ \left( \frac{\Delta\tau}{2} \right) f(t_i) + \left( \frac{\Delta\tau}{2} \right) f(t_i) \right] + 0(\Delta\tau) = f(t_i) + 0(t) \end{aligned} \qquad (A3) $$

Hence $\mathrm{Int} \to f(t_i)$ as $\Delta\tau \to 0$ and the convergence property

$$ \Delta_i(t) \to \delta(t - t_i) \qquad (A4) $$

or $\Delta_i(\cdot) \to \delta(\cdot - t_i) = \delta_i(\cdot)$ is now duly established.

Introducing $F(t)$ as a sampled version of $E(t)$:

$$ F(t) = \Delta\tau \sum_{i=1}^{N} E(t_i) \Delta_i(t) \qquad (A5) $$

Let us now proceed to show that as $\Delta\tau \to 0$

$$ F(t) \to \int_{-\infty}^{\infty} E(\tau) \delta(\tau - t) d\tau = E(t) \qquad (A6) $$

We shall make use of a well-known definition of the Riemann integral as a limit of sampling sums, sometimes referred to as Darboux sums, as
the sampling period $\Delta\tau$ tends to zero, namely for a well-bahaved⁵ (w.b.) function $f(t)$:

$$ \sum_{i=1}^{N} f(t_i) \Delta\tau \rightarrow \int_{-A}^{A} f(\tau) d\tau \qquad (A7) $$

as $\Delta\tau \rightarrow 0$ or $N \rightarrow \infty$

Let us now extend the previous definition and consider a sequence of functions $f(N, t)$ with $(N, t)$ converging to $f(t)$ as $N \rightarrow \infty$.⁹⁵ It can straightforwardly be shown that

$$ \sum_{i=1}^{N} f(N, t_i) \Delta\tau \rightarrow \int_{-A}^{A} f(\tau) d\tau \qquad (A8) $$

as $\Delta\tau \rightarrow 0$ or $N \rightarrow \infty$

Indeed,

$$ \sum_{i=1}^{N} f(N, t_i) \Delta\tau = \sum_{i=1}^{N} [f(N, t_i) - f(t_i)] \Delta\tau + \sum_{i=1}^{N} f(t_i) \Delta\tau \qquad (A9) $$

from strong convergence conditions, for any arbitrarily small $\epsilon$ value, there is an integer $N_0$ such that for any $N > N_0$, and any $i$, $|f(N, t_i) - f(t_i)| < \epsilon$. Then for $N > N_0$,

$$ \sum_{i=1}^{N} |f(N, t_i) - f(t_i)| \Delta\tau < 2\epsilon A \qquad (A10) $$

and therefore,

$$ \left| \sum_{i=1}^{N} [f(N, t_i) - f(t_i)] \Delta\tau \right| = \left| \sum_{i=1}^{N} [f(N, t_i) - f(t_i)] \Delta\tau \right| \le \sum_{i=1}^{N} |f(N, t_i) - f(t_i)| \Delta\tau < 2\epsilon A \qquad (A11) $$

$$ \sum_{i=1}^{N} f(t_i) \Delta\tau - 2\epsilon A \le \sum_{i=1}^{N} f(N, t_i) \Delta\tau \le \sum_{i=1}^{N} f(t_i) \Delta\tau + 2\epsilon A \qquad (A12) $$

Now for $N > N_1$, as per the Riemann integrability condition for $f$:

$$ \int_{-A}^{A} f(\tau) d\tau - \epsilon \le \sum_{i=1}^{N} f(t_i) \Delta\tau \le \int_{-A}^{A} f(\tau) d\tau + \epsilon \qquad (A13) $$

In conclusion, for $N > \sup(N_1, N_0)$ and arbitrary small $\eta = \epsilon(1 + 2A)$,

$$ \int_{-A}^{A} f(\tau) d\tau - \eta \le \sum_{i=1}^{N} f(N, t_i) \Delta\tau \le \int_{-A}^{A} f(\tau) d\tau + \eta \qquad (A14) $$

and hence the now proven property is:

$$ \sum_{i=1}^{N} f(N, t_i) \Delta\tau \rightarrow \int_{-A}^{A} f(\tau) d\tau \qquad (A15) $$

as $N \rightarrow \infty$ or $\tau \rightarrow 0$

Back to the main course of the demonstration, as $\Delta\tau \rightarrow 0$ or $N \rightarrow \infty$, $F$ converges to a Riemann integral, by virtue of the above-mentioned extension of the basic definition of the Riemann integral as:

$$ F \rightarrow \int_{-A}^{A} E(\tau) \delta(\tau - \tau) d\tau = E(t) \qquad (A16) $$

Let us start with a linear approximation for the response function and proceed to establish Equation 63:

$$ P^1[t, E(.)] = \int \tilde{R}^1(t; \tau) E(\tau) d\tau \qquad (A17) $$

$$ P^1[t, E(.)] = \lim_{N \to \infty} \{ P^1[t, F(.)] \} = \lim_{N \to \infty} \sum_{i=1}^{N} E(t_i) P^1[t, \Delta_i(.)] \Delta\tau \qquad (A18) $$

where we have made use of the linearity of $P^1[t, E(.)]$.

As $N \rightarrow \infty$, we know that $\Delta_i(\cdot) \rightarrow \delta(\cdot - t_i)$, hence $P^1[t, \Delta_i(.)] \rightarrow P^1[t, \delta(\cdot - t_i)]$ and, by virtue of the extended definition of the Riemann integral,

$$ P^1[t, E(.)] = \lim_{N \to \infty} \sum_{i=1}^{N} E(t_i) P^1[t, \Delta_i(.)] \Delta\tau = \int_{-A}^{A} E(\tau) P^1[t, \delta(\cdot - \tau)] d\tau \qquad (A19) $$

$$ P^1[t, E(.)] = \int_{-A}^{A} E(\tau) R^1(t, \tau) d\tau \qquad (A20) $$

where we have introduced the simplified notation $\tilde{R}^1(t, \tau) = P^1[t, \delta(\cdot - \tau)]$.

The interpretation of $R^1(t, u)$ appears now clearly as the linear impulse response of the system under a Dirac excitation at time $\tau$.

Most importantly, knowledge of the linear impulse response $\tilde{R}^1(t, \tau)$ allows the reconstruction of the output response of the linear system to any w.b. $E(.)$ input as from expression above.

Let us now extend this result to nonlinear responses beyond the linear approximation, leading us to introduce the multilinear impulse response $\tilde{R}^n(t; \tau_1, \tau_2, \dots, \tau_n)$ at order $n$ according to Equation 62 namely:

$$ P^n[t, E(.)] = \iiint \dots \int \tilde{R}^n(t; \tau_1, \tau_2, \dots, \tau_n) E(\tau_1) E(\tau_2) \dots E(\tau_n) d\tau_1 \dots d\tau_n \qquad (A21) $$

$$ P[t, E(.)] = \lim_{N \to \infty} \{P[t, F(.)]\} = \lim_{N \to \infty} P \left[ t, \sum_{i=1}^{N} E(t_i) \Delta_i(.) \Delta\tau \right] \qquad (A22) $$

$$ P \left[ t; \sum_{i=1}^{N} E_i \Delta_i(.) \Delta\tau \right] \sim P(t; 0) + \sum_{i=1}^{N} E_i \Delta\tau \left( \frac{\partial P[t; u_i \Delta_i(.)]}{\partial u_i} \right)_{u=0} \\ + \sum_{i=1}^{N} \sum_{j=1}^{N} E_i E_j \Delta\tau^2 \left( \frac{\partial^2 P[t; u_i \Delta_i(.) + u_j \Delta_j(.)]}{\partial u_i \partial u_j} \right)_{u=0} \\ + \sum_{i=1}^{N} \sum_{j=1}^{N} \sum_{k=1}^{N} E_i E_j E_k \Delta\tau^3 \left( \frac{\partial^3 P[t; u_i \Delta_i(.) + u_j \Delta_j(.) + u_k \Delta_k(.)]}{\partial u_i \partial u_j \partial u_k} \right)_{u=0} \dots \qquad (A23) $$
Note that we used here the following property, expressed here in a simplified form for an n-variable w.b. function f which cancels at the origin hereafter developed at first order:

$$f(x_1, \dots, x_n) = \sum_{j=1}^{n} x_j \left( \frac{\partial f}{\partial x_j} \right)_{x_1, \dots, x_n=0} = \sum_{j=1}^{n} x_j \frac{df(0, \dots, 0, x_j, 0 \dots 0)}{dx_j} \quad (A24)$$

This property can be easily established on a function of two variables and then straightforwardly generalized to n:

$$f(x, y) = x \left( \frac{\partial f(x, y)}{\partial x} \right)_{x,y=0} + y \left( \frac{\partial f(x, y)}{\partial y} \right)_{x,y=0} = x \left( \frac{\partial f(u, 0)}{\partial u} \right)_{u=0} \\ + y \left( \frac{\partial f(0, u)}{\partial u} \right)_{u=0} \quad (A25)$$

$$P \left[ t; \sum_{i=1}^{N} E_i \Delta_i(.) \Delta \tau \right] \rightarrow P[t; E(.)] = P(t; 0) \quad (A32)$$

$$+ \int_{-A}^{A} d\tau E(\tau) \left( \frac{\partial P[t; u\delta(\cdot, -\tau)]}{\partial u} \right)_{u=0} \quad (A33)$$

$$+ \int_{-A}^{A} \int_{-A}^{A} E(\tau_1)E(\tau_2)d\tau_1 d\tau_2 \left( \frac{\partial^2 P[t; u_1\delta(\cdot, -\tau_1) + u_2\delta(\cdot, -\tau_2)]}{\partial u_1 \partial u_2} \right)_{u_1, u_2=0} \\ + \int_{-A}^{A} \int_{-A}^{A} \int_{-A}^{A} E(\tau_1)E(\tau_2)E(\tau_3)d\tau_1 d\tau_2 d\tau_3 \left( \frac{\partial^3 P[t; u_1\delta(\cdot, -\tau_1) + u_2\delta(\cdot, -\tau_2) + u_3\delta(\cdot, -\tau_3)]}{\partial u_1 \partial u_2 \partial u_3} \right)_{u_1, u_2, u_3=0} \dots \quad (A34)$$

Indeed

$$\left(\frac{\partial f(x, y)}{\partial x}\right)_{x,y=0} = \lim_{\delta x \to 0} \left(\frac{f(x + \delta x, y) - f(x, y)}{\delta x}\right)_{x,y=0} \\ = \lim_{\delta x \to 0} \frac{f(\delta x, 0) - f(0, 0)}{\delta x} = \left(\frac{df(u, 0)}{du}\right)_{u=0} \quad (A26)$$

As $N \to \infty$,

$$\left(\frac{\partial P[t; u_i \Delta_i(.)]}{\partial u_i}\right)_{u=0} \rightarrow \left(\frac{\partial P[t; u_i \delta_i(.)]}{\partial u_i}\right)_{u=0} \quad (A27)$$

$$\left(\frac{\partial^2 P[t; u_i \Delta_i(.) + u_j \Delta_j(. )]}{\partial u_i \partial u_j}\right)_{u=0} \rightarrow \left(\frac{\partial^2 P[t; u_i \delta_i(.) + u_j \delta_j(.)]}{\partial u_i \partial u_j}\right)_{u=0} \quad (A28)$$

$$\left(\frac{\partial^3 P[t; u_i \Delta_i(.) + u_j \Delta_j(.) + u_k \Delta_k(.)}{\partial u_i \partial u_j \partial u_k}\right)_{u=0} \\ \rightarrow \left(\frac{\partial^3 P[t; u_i \delta_i(.) + u_j \delta_j(.) + u_k \delta_k(.)}{\partial u_i \partial u_j \partial u_k}\right)_{u=0} \quad (A29)$$

We make now use of the generalized form of the Riemann integral convergence, for n variable functions, according to

$$\sum_{j=1}^{N} \sum_{i=1}^{N} f(x_i, x_j) \Delta x^2 \rightarrow \int_{-A}^{A} \int_{-A}^{A} f(x, y) \, dxdy \quad (A30)$$

and likewise for higher orders, thereby introducing multiple integral expressions.

In the present case, this allows to express accordingly the relevant induced polarization limits according to:

$$P \left[ t; \sum_{i=1}^{N} E_i \Delta_i(.) \Delta \tau \right] \rightarrow P[t; E(.)] \quad (A31)$$

Also

Now calling:

$$\bar{R}^1(t; \tau) = \left( \frac{\partial P[t; u\delta(\cdot, -\tau)]}{\partial u} \right)_{u=0} \quad (A35)$$

$$\bar{R}^2(t; \tau_1, \tau_2) = \left( \frac{\partial^2 P[t; u_1\delta(\cdot, -\tau_1) + u_2\delta(\cdot, -\tau_2)]}{\partial u_1 \partial u_2} \right)_{u_1, u_2=0} \quad (A36)$$

$$\bar{R}^3(t; \tau_1, \tau_2, \tau_3) = \left( \frac{\partial^3 P[t; u_1\delta(\cdot, -\tau_1) + u_2\delta(\cdot, -\tau_2) + u_3\delta(\cdot, -\tau_3)]}{\partial u_1 \partial u_2 \partial u_3} \right)_{u_1, u_2, u_3=0} \quad (A37)$$

One ends-up with the desired expressions now set on firm grounds:

$$P[t; E(.)] = P[t; 0] + P^1[t; E(.)] + P^2[t; E(.)] + P^3[t; E(.)] + \dots \quad (A38)$$

$$P^1[t; E(.)] = \int \bar{R}^1(t; \tau) E(\tau) d\tau \quad (A39)$$

$$P^2[t; E(.)] = \iint \bar{R}^2(t; \tau_1, \tau_2) E(\tau_1) E(\tau_2) d\tau_1 d\tau_2 \quad (A40)$$
$$P^3[t, E(\cdot)] = \iint \bar{R}^3(t; \tau_1, \tau_2, \tau_3) E(\tau_1) E(\tau_2) E(\tau_3) d\tau_1 d\tau_2 d\tau_3 \text{ etc } \quad \text{(A41)}$$

$$P^1[t, E(\cdot)] = \int \left( \frac{\partial P[t; u\delta(\cdot, -\tau)]}{\partial u} \right)_{u=0} E(\tau) d\tau \quad \text{(A42)}$$

$$P^2[t, E(\cdot)] = \iint \left( \frac{\partial^2 P[t; u_1 \delta(\cdot, -\tau_1) + u_2 \delta(\cdot, -\tau_2)]}{\partial u_1 \partial u_2} \right)_{u_1, u_2=0} E(\tau_1) E(\tau_2) d\tau_1 d\tau_2 \quad \text{(A43)}$$

$$P^3[t, E(\cdot)] = \iint \left( \frac{\partial^3 P[t; u_1 \delta(\cdot, -\tau_1) + u_2 \delta(\cdot, -\tau_2) + u_3 \delta(\cdot, -\tau_3)]}{\partial u_1 \partial u_2 \partial u_3} \right)_{u_1, u_2, u_3=0} E(\tau_1) E(\tau_2) E(\tau_3) d\tau_1 d\tau_2 d\tau_3 \text{ etc } \quad \text{(A44)}$$

$$P^3[\vec{r}, t; E(\cdot)] = \iint \bar{R}^3(\vec{r}, t; \vec{\rho}_1, \vec{\rho}_2, \vec{\rho}_3; \tau_1, \tau_2, \tau_3) E(\vec{\rho}_1, \tau_1) E(\vec{\rho}_2, \tau_2) E(\vec{\rho}_3, \tau_3) d\vec{\rho}_1 d\vec{\rho}_2 d\vec{\rho}_3 d\tau_1 d\tau_2 d\tau_3 \text{ etc } \quad \text{(A51)}$$

$$\bar{R}^1(\vec{r}, t; \vec{\rho}, \tau) = \left( \frac{\partial P[\vec{r}, t; u\delta(\cdot, -\vec{\rho}, -\tau)]}{\partial u} \right)_{u=0} \quad \text{(A52)}$$

$$\bar{R}^2(\vec{r}, t; \vec{\rho}_1, \vec{\rho}_2, \tau_1, \tau_2) = \left( \frac{\partial^2 P[\vec{r}, t; u_1 \delta(\cdot, -\vec{\rho}_1, -\tau_1) + u_2 \delta(\cdot, -\vec{\rho}_2, -\tau_2)]}{\partial u_1 \partial u_2} \right)_{u_1, u_2=0} \quad \text{(A53)}$$

$$\bar{R}^3(\vec{r}, t; \vec{\rho}_1, \vec{\rho}_2, \vec{\rho}_3; \tau_1, \tau_2, \tau_3) = \left( \frac{\partial^3 P[\vec{r}, t; u_1 \delta(\cdot, -\vec{\rho}_1, -\tau_1) + u_2 \delta(\cdot, -\vec{\rho}_2, -\tau_2) + u_3 \delta(\cdot, -\vec{\rho}_3, -\tau_3)]}{\partial u_1 \partial u_2 \partial u_3} \right)_{u_1, u_2, u_3=0} \quad \text{(A54)}$$

The above result may be reached by way of a shortcut allowing to bypass the introduction of the $\Delta_i$ functions, then directly considering the sole $\delta_i$ Dirac distributions.

Indeed, the $E(\cdot)$ input can be expressed as

$$E(\cdot) = \int E(\tau) \delta(\cdot, -\tau) d\tau \quad \text{(A45)}$$

where $\delta(\cdot, -\tau)$ stands for the Dirac distribution at $t = \tau$.

Then

$$P[t; E(\cdot)] = P\left[t; \int E(\tau) \delta(\cdot, -\tau) d\tau\right] \quad \text{(A46)}$$

Applying the definition of the Riemann integral, here in a reverse mode:

$$\int E(\tau) \delta(\cdot, -\tau) d\tau = \lim_{N \to \infty} P \left[ t, \sum_{i=1}^{N} E(\tau_i) \delta(\cdot, -\tau_i) \Delta \tau \right] \quad \text{(A47)}$$

One can then expand the above $P$ expression in terms of the $E(\tau_i)\Delta\tau$ values by way of a Taylor expansion, using a similar step as in the above derivation, and finally transform the discrete summations into integrals that have the same expressions as obtained by the first longer derivation.

The derivations above can be straightforwardly generalized to time and space, $E(\cdot)$ standing now for a function of space and time locally given by $E(\vec{r}, t)$:

$$P[\vec{r}, t; E(\cdot)] = P[\vec{r}, t; 0] + P^1[\vec{r}, t; E(\cdot)] + P^2[\vec{r}, t; E(\cdot)] + P^3[\vec{r}, t; E(\cdot)] + \dots \quad \text{(A48)}$$

$$P^1[\vec{r}, t; E(\cdot)] = \int \bar{R}^1(\vec{r}, t; \vec{\rho}, \tau) E(\vec{\rho}, \tau) d\vec{\rho} d\tau \quad \text{(A49)}$$

$$P^2[\vec{r}, t; E(\cdot)] = \iint \bar{R}^2(\vec{r}, t; \vec{\rho}_1, \vec{\rho}_2, \tau_1, \tau_2) E(\vec{\rho}_1, \tau_1) E(\vec{\rho}_2, \tau_2) d\vec{\rho}_1 d\vec{\rho}_2 d\tau_1 d\tau_2 \quad \text{(A50)}$$

where $\delta(\cdot, -\vec{\rho}, -\tau)$ is the multivariable Dirac distribution at pace time position $(\vec{\rho}, \tau)$ defined for any well behaved $E(\vec{r}, t)$ by:

$$\langle E(\vec{r}, t) | \delta(\cdot, -\vec{\rho}, -\tau) \rangle = \iiint E(\vec{r}, t) \delta(\vec{r} - \vec{\rho}, t - \tau) d\vec{\rho} d\tau = E(\vec{\rho}, \tau) \quad \text{(A55)}$$

Derivatives over $u$ such as

$$\frac{\partial P[\vec{r}, t; u\delta(\cdot, -\vec{\rho}, -\tau)]}{\partial u} \quad \text{(A56)}$$

have to be understood as distribution derivatives[84]

## Appendix B: Multiple Field Tensor

In the general case of a two-step, possibly inter-twinned, of write and read configurations, both assumed to be nonlinear, the upstream functional can be expressed:

$$\vec{P}(t; \vec{x}) = \vec{\bar{F}}_{(p,q)}^{(r,w)} [\vec{E}_1^{(r)}(t, \vec{x}) \dots \vec{E}_p^{(r)}(t, \vec{x}) \dots \vec{E}_q^{(w)}] \quad \text{(B1)}$$

where $\vec{E}_i^{(r)}, 1 \le i \le p$ (resp. $\vec{E}_q^{(w)}, 1 \le i \le q$) stand for read- (resp. write-) electric field vectors. One can then follow the derivation above extended to an increased number of input electric fields, leading to a generalized response function ultimately expressed in the Fourier domain to be further developed.

We therefore assume that the write-read process is referred to a single but dual impulse response dependent on 4q write process and 4p time-space variables related to the probe process, namely:

$$\vec{E}_{(p,q)}^{(r,w)} (\tau_1^r, \tau_1^r, \dots, \tau_p^r, \tau_p^r, \tau_1^w, \tau_1^w, \dots, \tau_q^w, \tau_q^w) \quad \text{(B2)}$$

The write-read tensor is correspondingly expressed a:

$$\bar{F}_{(p,q)}^{(r,w)} = \bar{F}_{(p)}^{(r)} \otimes \bar{E}_{(q)}^{(w)} = \bar{u}^{(r)} \otimes \bar{E}_{1}^{(r)}(\tau_{1}^{r}, \tau_{1}^{r}) \otimes \dots \otimes \bar{E}_{p}^{(r)}(\tau_{p}^{r}, \tau_{p}^{r}) \otimes \\ \bar{E}_{1}^{(w)}(\tau_{1}^{w}, \tau_{1}^{w}) \otimes \dots \otimes \bar{E}_{q}^{(w)}(\tau_{q}^{w}, \tau_{q}^{w}) \quad \text{(B3)}$$
The combined write-read process can then be condensed in a single multiple convolution response:

$$ \vec{P}(t, \vec{r}) = \vec{\mathcal{L}}_{(p,q)}^{(r,w)} \odot \vec{\mathcal{F}}_{(p,q)}^{(r,w)} \quad (B4) $$

where the ⊙ symbol stands for the combination of a multiple convolution of the two tensorial functions over the 4(p + q) space-time set of variables together with a full tensor contraction (e.g., generalized scalar product in the Euclidian space) over Cartesian indices. The Cartesian components of the impulse response tensor $\vec{\mathcal{L}}_{(p,q)}^{(r,w)}$ of rank p + q + 1 are the $3^{p+q+1}$ complex numbers $\vec{\mathcal{L}}_{w_1...u_1...u_p, \nu_1...\nu_q}$ with its Cartesian indices w, u₁... uₚ, ν₁... νq spanning independently the x, y, z directions in the chosen Cartesian reference frame.

Likewise for the dual write-read field tensor $\vec{\mathcal{F}}_{(p,q)}^{(r,w)}$ which is correspondingly made-up of the matching components $\vec{\mathcal{F}}_{w_1...u_1...u_p, \nu_1...\nu_q}^{(r,w)} = u_\omega E_{u_1}^{(r)} ... E_{u_p}^{(r)} E_{\nu_1}^{(w)} ... E_{\nu_q}^{(w)}$ and

trated by the all-optical poling configuration which entails a sequence of write and read fields, with a rich literature devoted to the domain (Refs.).

Upstream to these considerations stands a response functional, depending on write and read electric field vectors expressed in time-space, namely:

$$ \vec{P}(t; \vec{x}) = \vec{\mathcal{F}} \left[ \vec{E}^r(t, \vec{x}), \vec{E}^w(t, \vec{x}) \right] = \vec{\mathcal{F}}[\vec{E}^r(t, \vec{x}), \vec{E}^w(t, \vec{x})] + \vec{E}^{wf}(t, \vec{x}) + \vec{E}^{wh}(t, \vec{x}) \quad (B8) $$

Here the write electric field is the combination of interfering fundamental and harmonic input electric fields that we prefer to express at this upstream stage in this unusual time-space form, so as to leave room for further elaboration in terms of short pulse.

Likewise, the seminal work of Baranova and Zeld'ovitch[70] is based on time-dependent expressions, upstream to the more common frequency-dependent expression $(E^{2\omega})^n (E^{\omega})^2$.

One can then follow the general write-read derivation detailed above in the case of AOP configurations, ending with

$$ (\vec{P})_u (2\omega^{(r)}, 2\vec{k}^r) = \sum_{w_1 u_1...u_p, \nu_1...\nu_q} \left\{ \begin{aligned} & (\vec{u})_w \mathcal{L}_{w_1 u_1 u_2, \nu_1 \nu_2 \nu_3}^{(r,w)} (\omega^{(r)}, \vec{k}^r; \omega^{(w)}, \vec{k}^w) \\ & \dots E_{u_1}^{(r)} (\omega^{(r)}, \vec{k}^r) E_{u_2}^{(r)} (\omega^{(r)}, \vec{k}^r) E_{\nu_1}^{(w)} (2\omega^{(w)}, 2\vec{k}^w) E_{\nu_2}^{(w)} (\omega^{(w)}, \vec{k}^w) E_{\nu_3}^{(w)} (\omega^{(w)}, \vec{k}^w) \end{aligned} \right\} \quad (B9) $$

$$ (\vec{P})_u = \vec{P} \cdot \vec{u} = \text{cont} \left( \vec{\mathcal{L}}_{(p,q)}^{(r,w)} \vec{\mathcal{F}}_{(p,q)}^{(r,w)} \right) = \sum_{w_1 u_1...u_p, \nu_1...\nu_q} (\vec{\mathcal{L}}_{w_1 u_1...u_p, \nu_1...\nu_q}^{(r,w)})^* \left( u_\omega E_{u_1}^{(r)} ... E_{u_p}^{(r)} E_{\nu_1}^{(w)} ... E_{\nu_q}^{(w)} \right) \quad (B5) $$

where notations are self-understandingly cont(A, B) for the full contraction of the same rank Cartesian tensors A and B and the star stands for the multiple convolutions, namely

The field tensor can be expressed in condensed form in the Fourier space as:

$$ \begin{aligned} \mathcal{F}_{\text{AOP}}^{(w,r)} = & \dots \vec{u} \otimes E^{(r)}(\omega^{(r)}, \vec{k}^r) \otimes E^{(r)}(\omega^{(r)}, \vec{k}^r) \otimes (E^{(w)})^* (2\omega^{(w)}, 2\vec{k}^w) \otimes \\ & E^{(w)}(\omega^{(w)}, \vec{k}^w) \otimes E^{(w)}(\omega^{(w)}, \vec{k}^w) \end{aligned} \quad (B10) $$

$$ (\vec{P})_u (t, \vec{r}) = \left\{ \begin{aligned} & \sum_{w_1 u_1...u_p, \nu_1...\nu_q} \int d\tau_1^r \int d\tau_1^r \dots \int d\tau_p^r \int d\tau_p^r \int d\tau_1^w \int d\tau_1^w \dots \int d\tau_q^w \int d\tau_q^w \dots \\ & \{ (\vec{u})_w [\vec{\mathcal{L}}_{w_1 u_1...u_p, \nu_1...\nu_q}^{(r,w)} (\tau_1^r, \vec{r}_1^r, \dots, \tau_p^r, \vec{r}_p^r; \tau_1^w, \vec{r}_1^w, \dots, \tau_q^w, \vec{r}_q^w) \dots \\ & \dots E_{u_1}^{(r)} (t - \tau_1^r, \vec{r} - \vec{r}_1^r) \dots E_{u_p}^{(r)} (t - \tau_p^r, \vec{r} - \vec{r}_p^r) \dots \\ & \dots E_{\nu_1}^{(w)} (t - \tau_1^w, \vec{r} - \vec{r}_1^w) \dots E_{\nu_q}^{(w)} (t - \tau_q^w, \vec{r} - \vec{r}_q^w) \} \end{aligned} \right\} \quad (B6) $$

Leading to the more usual form after a multiple Fourier transform:

$$ (\vec{P})_u (\omega, \vec{k}) = \sum_{w_1 u_1...u_p, \nu_1...\nu_q} \left\{ \begin{aligned} & (\vec{u})_w \mathcal{L}_{w_1 u_1...u_p, \nu_1...\nu_q}^{(r,w)} (\omega_1^r, \vec{k}_1^r, \dots, \omega_p^r, \vec{k}_p^r; \omega_1^w, \vec{k}_1^w, \dots, \omega_q^w, \vec{k}_q^w) \dots \\ & \dots E_{u_1}^{(r)} (\omega_1^r, \vec{k}_1^r) \dots E_{u_p}^{(r)} (\omega_p^r, \vec{k}_p^r) E_{\nu_1}^{(w)} (\omega_1^w, \vec{k}_1^w) \dots E_{\nu_q}^{(w)} (\omega_q^w, \vec{k}_q^w) \end{aligned} \right\} \quad (B7) $$

The relevance of combining both read-and write-field tensors under a single tensor together with the corresponding response function is illus-

In the commonly encountered configuration whereby the write and read beam are generated in a fully collinear configuration by a single
fundamental laser source with an additional frequency doubling contribution to the write electric field, the above expression simplifies into:

$$ T_{AOP}^{(w,r)} = \dots \bar{u} \otimes E^{(r)}(\omega, \vec{k}) \otimes E^{(r)}(\omega, \vec{k}) \otimes [E^{(w)}(2\omega, 2\vec{k})]^* \otimes \\ E^{(w)}(\omega, \vec{k}) \otimes E^{(w)}(\omega, \vec{k}) \qquad (B11) $$

In the above, we leave aside a discussion of the relevant phase factors contributing significantly to the different electric fields above, herein implicitly contained in their amplitudes. These phase factors play a crucial role in the read process, as evidenced by propagative phase matching considerations to be developed in the part of this paper devoted to crystals. In the case of AOP processes and related ones, external phases have been identified early on as playing a crucial role in the formation of nonlinear gratings and allowing for the control of their diffraction efficiency by way of phase control devices.

Indeed, a rank-6 $\chi^{(5)}$ electric susceptibility tensor has been introduced and studied by means of nonlinear ellipsometry,[96] followed by the rich literature of AOP. Moreover, this work explores the particularly noteworthy case of a liquid where orientational relaxation of nonlinear solute molecules plays a central role, moreover in a AOP induced holographic grating, both pointing at the relevance of a time-space frame to start with.

In many studies, this dual write- and read- configuration is interwinned with an appropriate duty cycle depending on both, the relaxation constant of the medium such as pertain to its viscosity and related physical-chemical parameters, and on the relaxation times of the excited states in the spectrum that contribute to the NL write and read processes for the involved atoms or molecules. Note that in general the write beams must be tuned to the excited state of the molecule, so as to further trigger the Weigert effect, namely the diazo-bond photo-isomerization responsible for the statistical molecular randomization, whereas the read fundamental beam is not necessarily tuned to the excited state and furthermore, in cases involving a significant propagation, should not be resonant so as to avoid spurious reabsorption effects.

## Conflict of Interest

The authors declare no conflict of interest.

## Keywords

crystal optics, imaging, microscopy, nonlinear optics, symmetries, tensor

Received: March 25, 2024

Revised: July 11, 2024

Published online:

[1] P. N. Butcher, D. Cotter, *The Elements of Nonlinear Optics*, Cambridge University Press, Cambridge, 1990.

[2] R. W. Boyd, *Nonlinear Optics*, Elsevier, Amsterdam, 2008.

[3] P. D. Maker, *Phys. Rev. A* **1970**, 1, 923.

[4] R. Bonneville, D. S. Chemla, *Phys. Rev. A* **1978**, 17, 2046.

[5] R. W. Boyd, D. Prato, *Nonlinear Optics*, Third Edition, Academic Press, Amsterdam, 2008.

[6] B. Boulanger, J. Zyss, in *International Tables for Crystallography Volume D: Physical Properties of Crystals*, (Ed: A. Authier), Springer Netherlands, Dordrecht, 2014, pp. 181.

[7] D. L. Andrews, *Symmetry* **2020**, 12, 1466.

[8] D. L. Andrews, *J. Chem. Phys.* **2023**, 158, 034101.

[9] J. Jerphagnon, D. Chemla, R. Bonneville, *Adv. Phys.* **1978**, 27, 609.

[10] U. Fano, G. Racah, *Irreducible Tensorial Sets*, Academic Press, Cambridge, 1959.

[11] J. Duboisset, H. Rigneault, S. Brasselet, *Phys. Rev. A* **2014**, 90, 063827.

[12] C. Cleff, H. Rigneault, S. Brasselet, J. Duboisset, *Phys. Rev. A* **2017**, 96, 013851.

[13] C. Cleff, A. Gasecka, P. Ferrand, H. Rigneault, S. Brasselet, J. Duboisset, *Nat. Commun.* **2016**, 7.

[14] F. N. H. Robinson, *Bell Syst. Tech. J.* **1967**, 46, 913.

[15] H. Chen, M. Liu, T. Yan, *Commun. Theor. Phys.* **2020**, 72, 075503.

[16] S. Brasselet, *Adv. Opt. Photonics* **2011**, 3, 205.

[17] D. Aït-Belkacem, M. Guilbert, M. Roche, J. Duboisset, P. Ferrand, G. Sockalingum, P. Jeannesson, S. Brasselet, *J. Biomed. Opt.* **2012**, 17, 080506.

[18] J. Duboisset, D. Aït-Belkacem, M. Roche, H. Rigneault, S. Brasselet, *Phys. Rev. A* **2012**, 85, 043829.

[19] K. Tilbury, C.-H. Lien, S.-J. Chen, P. J. Campagnola, *Biophys. J.* **2014**, 106, 354.

[20] F.-Z. Bioud, P. Gasecka, P. Ferrand, H. Rigneault, J. Duboisset, S. Brasselet, *Phys. Rev. A* **2014**, 89, 013836.

[21] P. Gasecka, A. Jaouen, F.-Z. Bioud, H. B. de Aguiar, J. Duboisset, P. Ferrand, H. Rigneault, N. K. Balla, F. Debarbieux, S. Brasselet, *Bio-phys. J.* **2017**, 113, 1520.

[22] B. Boulanger, G. Marnier, *Opt. Commun.* **1990**, 79, 102.

[23] B. Boulanger, G. Marnier, *J. Phys.: Condens. Matter* **1991**, 3, 8327.

[24] B. Boulanger, J. P. Fève, G. Marnier, *Phys. Rev. E* **1993**, 48, 4730.

[25] G. Lerner, O. Neufeld, L. Hareli, G. Shoulga, E. Bordo, A. Fleischer, D. Podolsky, A. Bahabad, O. Cohen, *Sci. Adv.* **2023**, 9, eade0953.

[26] J. Zyss, I. Ledoux, *Chem. Rev.* **1994**, 94, 77.

[27] S. Brasselet, J. Zyss, *J. Opt. Soc. Am., B: Opt. Phys.* **1998**, 15, 257.

[28] G. Placzek, *The Rayleigh and Raman Scattering*, United States Atomic Energy Commission, Division Of Technical Information, Oak Ridge, Tennessee, 1962.

[29] B. F. Gächter, *J. Mol. Spectrosc.* **1976**, 63, 1.

[30] J. Nestro, T. G. Spiro, *J. Raman Spectrosc.* **1973**, 1, 539.

[31] L. Hecht, L. A. Nafie, *Mol. Phys.* **1991**, 72, 441.

[32] R. Clark, A. J. McCaffery, *J. Phys. Chem.* **1977**, 81, 1918.

[33] A. J. McCaffery, R. A. Shatwell, *Rev. Sci. Instrum.* **1976**, 47, 247.

[34] R. Clark, S. R. Jeyes, A. J. McCaffery, R. A. Shatwell, *J. Am. Chem. Soc.* **1974**, 96, 5586.

[35] J. Zyss, I. Ledoux, J.-F. Nicoud, in *Molecular Nonlinear Optics* (Ed: J. Zyss), Academic Press, San Diego, 1994.

[36] M. Kauranen, A. V. Zayats, *Nat. Photonics* **2012**, 6, 737.

[37] P. Dombi, Z. Pápa, J. Vogelsang, S. V. Yalunin, M. Sivis, G. Herink, S. Schäfer, P. Groß, C. Ropers, C. Lienau, *Rev. Mod. Phys.* **2020**, 92, 025003.

[38] J. Butet, P.-F. Brevet, O. J. F. Martin, *ACS Nano* **2015**, 9, 10545.

[39] G. Bachelier, J. Butet, I. Russier-Antoine, C. Jonin, E. Benichou, P.-F. Brevet, *Phys. Rev. B* **2010**, 82, 235403.

[40] K. Konishi, T. Higuchi, J. Li, J. Larsson, S. Ishii, M. Kuwata-Gonokami, *Phys. Rev. Lett.* **2014**, 112, 135502.

[41] A. Salomon, Y. Prior, M. Fedoruk, J. Feldmann, R. Kolkowski, J. Zyss, *J. Opt.* **2014**, 16, 114012.

[42] S. Chen, G. Li, F. Zeuner, W. H. Wong, E. Y. B. Pun, T. Zentgraf, K. W. Cheah, S. Zhang, *Phys. Rev. Lett.* **2014**, 113, 033901.

[43] D. Oron, E. Tal, Y. Silberberg, *Opt. Lett.* **2003**, 28, 2315.

[44] N. Olivier, F. Aptel, K. Plamann, M.-C. Schanne-Klein, E. Beaurepaire, *Opt. Express* **2010**, OE 18, 5028.

[45] Z. Özgo, S. Kielich, *Physica B* **1976**, 81, 151.

[46] N. Ben-Tal, N. Moiseyev, A. Beswick, *J. Phys. B: At., Mol. Opt. Phys.* **1993**, 26, 3017.

[47] O. E. Alon, V. Averbukh, N. Moiseyev, *Phys. Rev. Lett.* **1998**, 80, 3743.

[48] C. L. Tang, H. Rabin, *Phys. Rev. B* **1971**, 3, 4025.

[49] N. Saito, P. Xia, F. Lu, T. Kanai, J. Itatani, N. Ishii, *Optica* **2017**, 4, 1333.

[50] A. Fleischer, O. Kfir, T. Diskin, P. Sidorenko, O. Cohen, *Nat. Photonics* **2014**, 8, 543.

[51] D. B. Milošević, W. Becker, R. Kopold, *Phys. Rev. A* **2000**, 61, 063403.
[52] O. Kfir, P. Grychtol, E. Turgut, R. Knut, D. Zusin, D. Popmintchev, T. Popmintchev, H. Nembach, J. M. Shaw, A. Fleischer, H. Kapteyn, M. Murnane, O. Cohen, *Nat. Photonics* **2015**, *9*, 99.

[53] F. Mauger, A. D. Bandrauk, T. Uzer, *J. Phys. B: At., Mol. Opt. Phys.* **2016**, *49*, 10LT01.

[54] A. Gorlach, M. E. Tzur, M. Birk, M. Krüger, N. Rivera, O. Cohen, I. Kaminer, *Nat. Phys.* **2023**, *19*, 1689.

[55] A. Zumbusch, G. R. Holtom, X. S. Xie, *Phys. Rev. Lett.* **1999**, *82*, 4142.

[56] M. Baertsch, P. Bornhauser, G. Calzaferri, R. Imhof, *J. Phys. Chem.* **1994**, *98*, 2817.

[57] P. K. Upputuri, J. Lin, L. Gong, X.-Y. Liu, H. Wang, Z. Huang, *Opt. Lett.* **2013**, *38*, 1262.

[58] B. A. Garetz, S. Arnold, *Opt. Commun.* **1979**, *31*, 1.

[59] O. Faucher, E. Prost, E. Hertz, F. Billard, B. Lavorel, A. A. Milner, V. A. Milner, J. Zyss, I. S. Averbukh, *Phys. Rev. A* **2016**, *94*, 051402.

[60] E. Skantzakis, S. Chatziathanasiou, P. A. Carpeggiani, G. Sansone, A. Nayak, D. Gray, P. Tzallas, D. Charalambidis, E. Hertz, O. Faucher, *Sci. Rep.* **2016**, *6*, 39295.

[61] R. Loucifsaibi, K. Nakatani, J. Delaire, M. Dumont, Z. Sekkat, *Chem. Anal. Biol. Fate: Polynucl. Aromat. Hydrocarbons, Int. Symp.*, 5th **1993**, *5*, 229.

[62] A.-C. Etilé, C. Fiorini, F. Charra, J.-M. Nunzi, *Phys. Rev. A* **1997**, *56*, 3888.

[63] O. Gayer, Z. Sacks, E. Galun, A. Arie, *Appl. Phys. B* **2008**, *91*, 343.

[64] T. Toury, S. Brasselet, J. Zyss, *Opt. Lett.* **2006**, *31*, 1468.

[65] R. Castagna, A. Milner, J. Zyss, Y. Prior, *Adv. Mater.* **2013**, *25*, 2234.

[66] “NLOQO 43.1-4, p. 97–131 – Old City Publishing,” can be found under. https://www.oldcitypublishing.com/journals/nloqo-home/nloqo-issue-contents/nloqo-volume-43-number-1-4-2012/nloqo-43-1-4-p-97-131/, **n.d.**

[67] M. Jarema, A. C. Mituś, J. Zyss, *Sci. Rep.* **2021**, *11*, 2234.

[68] F. Michelotti, E. Toussaere, R. Levenson, J. Liang, J. Zyss, *Appl. Phys. Lett.* **1995**, *67*, 2765.

[69] W. Margulis, U. Österberg, *J. Opt. Soc. Am. B* **1988**, *5*, 312.

[70] N. B. Baranova, A. N. Chudinov, B. Ya, *Opt. Commun.* **1990**, *79*, 116.

[71] R. J. Glauber, *Phys. Rev.* **1963**, *130*, 2529.

[72] S. Bidault, S. Brasselet, J. Zyss, *Opt. Lett.* **2004**, *29*, 1242.

[73] M. Dumont, *J. Opt. Soc. Am. B: Opt. Phys.* **2009**, *26*, 1057.

[74] R. Piron, S. Brasselet, D. Josse, J. Zyss, G. Viscardi, C. Barolo, *J. Opt. Soc. Am. B* **2005**, *22*, 1276.

[75] S. Brasselet, J. Zyss, *Opt. Lett.* **1997**, *22*, 1464.

[76] J. A. Armstrong, N. Bloembergen, J. Ducuing, P. S. Pershan, *Phys. Rev.* **1962**, *127*, 1918.

[77] F. C. Zumsteg, J. D. Bierlein, T. E. Gier, *J. Appl. Phys.* **1976**, *47*, 4980.

[78] B. Boulanger, J. P. Fève, G. Marnier, B. Ménaert, X. Cabirol, P. Villeval, C. Bonnin, *J. Opt. Soc. Am. B: Opt. Phys.* **1994**, *11*, 750.

[79] V. Petrov, V. V. Badikov, D. V. Badikov, K. Kato, G. S. Shevyrdyaeva, K. Miyata, M. Mero, L. Wang, Z. Heiner, V. L. Panyutin, *J. Opt. Soc. Am. B* **2021**, *38*, B46.

[80] F. Guo, E. Boursier, P. Segonds, A. Peña, J. Debray, V. Badikov, V. Panyutin, D. Badikov, V. Petrov, B. Boulanger, *Opt. Lett.* **2022**, *47*, 842.

[81] E. Boursier, P. Segonds, J. Debray, P. L. Inácio, V. Panyutin, V. Badikov, D. Badikov, V. Petrov, B. Boulanger, *Opt. Lett.* **2015**, *40*, 4591.

[82] R. C. Miller, *Appl. Phys. Lett.* **2004**, *5*, 17.

[83] Rank and order for a given tensor are different as per conventions: n stands for the rank (i.e. number of Cartesian indices) and n - 1 stands for the order (i.e. degree of the monomial product of input components). In this section, the tensors are labelled with the rank n, instead of the order n - 1.

[84] L. Schwartz, *Théorie Des Distributions*, Hermann, Paris, **1997**.

[85] A. Karnieli, A. Arie, *Opt. Express* **2018**, *26*, 4920.

[86] W. A. Benalcazar, B. A. Bernevig, T. L. Hughes, *Science* **2017**, *357*, 61.

[87] W. A. Benalcazar, B. A. Bernevig, T. L. Hughes, *Phys. Rev. B* **2017**, *96*, 245115.

[88] F. Hache, M.-C. Schanne-Klein, H. Mesnil, M. Alexandre, G. Lemercier, C. Andraud, C. R. *Phys.* **2002**, *3*, 429.

[89] A.-M. Pena, T. Boulesteix, T. Dartigalongue, M.-C. Schanne-Klein, *J. Am. Chem. Soc.* **2005**, *127*, 10314.

[90] P. Woźniak, I. D. Leon, K. Höflich, G. Leuchs, P. Banzer, *Optica* **2019**, *6*, 961.

[91] We chose to use here the term kernel to distinguish the integral response, such as for the linear expression $\int \tilde{R}^1(t; \tau)E(\tau)d\tau$ where the time shift runs from $-\infty$ to $+\infty$, that is over past as well as future with no causality constraints at this stage, from the convolution response integral under $\int \tilde{R}^1(t - \tau)E(\tau)d\tau$ the assumption that the response is invariant by translation in time, to be further explained.

[92] We could have chosen other sampling functions such as featuring a triangular shape provided the convergence to a Dirac function is ensured.

[93] $\Delta\tau \to 0$ is equivalent to $N \to \infty$ as $\Delta\tau = 2A/N$.

[94] By a well behaved (w.b.) function, we refer here to one that is at least continuous and if needed, differentiable to any order needed to validate the demonstration.

[95] Assuming here «strong» convergence of the $f_i$'s towards $f$, that is $\sup[f] \text{ over } [a,b]$ tends to zero. “Weaker” conditions would refer to a point by point convergence whereby at all times $t_i f_i(t) \to f(t)$, indeed a weaker condition than the former one where global convergence is ensured. A more rigorous demonstration of the above would be needed to decide on the convergence condition needed to ensure the extended form of the Riemann integral definition. As this would not fall within the scope of this review, and we decided to opt for the stronger of the two conditions.

[96] J.-M. Nunzi, F. Charra, C. Fiorini, J. Zyss, *Chem. Phys. Lett.* **1994**, *219*, 349.

<figure><img src="image_12.png" /></figure>

Julien Duboisset is assistant professor at the Aix-Marseille University, France. He has obtained his Ph.D. in 2009 at the University Claude Bernard Lyon 1 France, on the study of the second harmonic generation of biomolecules. He then spent a year as a postdoctoral fellow at NTNU (Norway), working on the molecular structure of amyloids using fluorescence microscopy. Since 2010, he conducts research at the Fresnel Institute, where he develops polarization-resolved nonlinear optical methods for imaging and spectroscopy. His research interests include the molecular-scale organization of water and electrolytes, as well as biomineralization.
<figure><img src="image_13.png" /></figure>

**Benoit Boulanger** received the Ph.D. degree in material sciences from the University of Nancy in 1989. He was researcher at the National Center for Scientific Research (CNRS) from 1989 to 2000 at the University of Nancy, Stanford University, and the University of Bourgogne. In 2000, he joined Grenoble-Alpes University as professor. He does his research at Institut Néel CNRS. He got a position of distinguished professor at Tianjin University of Technology (2018–2021). His work is at the frontiers between nonlinear crystal optics, material engineering and quantum optics. He has authored over 280 papers in refereed journals and conference proceedings.

<figure><img src="image_14.png" /></figure>

**Sophie Brasselet** is an optical physicist. She has obtained her Ph.D. in 1997 at the University of Paris Sud France, on the study of multipolar molecules applied to nonlinear optics, and spent two years as a postdoctoral fellow at UCSD and Stanford University (USA), working on the detection of single molecules in cells. She is now a CNRS research director at Fresnel Institute in Marseille, France. She has developed methods in non-linear optical microscopy and fluorescence imaging, based on the control of light polarization, with pioneering contributions for structural nanoscale imaging in biological samples.

<figure><img src="image_15.png" /></figure>

**Patricia Segonds** received her Ph.D. degree in physics from the University Grenoble Alpes followed by a post-doctoral fellowship at Ohio State University. In 1989 she was appointed associate professor at the University of Bordeaux. Her research at the Centre de Physique Moléculaire et d'Optique Hertzienne (CPMOH) focused on third-order nonlinearities and ultrafast all-optical switching in glasses. She joined the University Grenoble Alpes in 1998 and is currently professor. She extended her expertise to frequency conversion from second-order nonlinearities in new crystals and the implementation of optical parametric oscillators, at the Laboratoire Interdisciplinaire de Physique (LIPHY) and since 2007 at NEEL Institute.

<figure><img src="image_16.png" /></figure>

**Joseph Zyss** is an emeritus professor in Physics at Ecole Normale Supérieure Paris-Saclay since 2015. He graduated from Ecole Polytechnique (1969–72) and earned a Ph.D. in physics at Pierre et Marie Curie University in Paris (1982). He has played a pioneering role in molecular, crystalline and polymer based nonlinear optics with related photonics applications, at nano-, micro-, and nanoscales, based on an original mix of algebraic, chemical and physical concepts gearing toward new directions in organic photonics and related fields. He has been at the Centre National d'Etudes des Télécommunications from 1975 to 1996 and then joined Ecole Normale Supérieure Paris Saclay, where he founded and directed the Laboratory for Molecular Quantum Photonics, the multidisciplinary D'Alembert Institute and a joint laboratory between CNRS and the Weizmann (2008–2015). He has been a recipient of the Meierof fellowship from the Weizmann Institute, the IBM and Yves Rocard prizes from the French Physical Society, the French German Humboldt-Gay Lussac award and is a fellow of the Optical Society of America.