1. INTRODUCTION

In this last part, we give voice to the computer by presenting a case of strong coupling between the different frequencies and by showing a case where the nonlinearity is such that we are no longer in the comfort zone of the pump wave nondepletion hypothesis.

2. NUMERICAL METHODS

A. Generalities

In this section, it is obviously out of the question to deal with the numerical solution of linear or nonlinear partial differential equations (PDEs): a book would not be enough. We implicitly suppose that we have at our disposal a finite element software capable of correctly solving a PDE (or a system of PDEs) in a domain of finite extension of 2D or 3D space, with suitable conditions at its boundary. This poses a double problem of infinity. On the one hand, fields have an unfortunate tendency to decrease very slowly: it is out of the question to mesh an area such that at its boundary the field is considered negligible. On the other hand, the sources are sometimes very far (in terms of wavelength) from the diffracting object, even to infinity if we consider that the incident wave is a plane wave. It will be necessary to find tricks to get around these problems.

B. Perfectly Matched Layers

The first problem that must be solved is to avoid meshing an area that is too large. This problem has been one of the key points and has contributed to the fact that for a long time the finite element method in the field of electromagnetic diffraction has remained ineffective. The introduction of perflectly matched layers (PMLs) in the 1990s (see, for instance, the seminal paper [1]) was decisive. This subsequently made it possible to solve electromagnetic diffraction problems using methods involving only finite elements, which, after discretization, made it possible to handle only sparse matrices. This technique has been applied in many contexts, for instance, for the study of gratings (e.g., [2]).

 

Fig. 1. Principle of the virtual antenna. (a) The scattering problem is described by a scatterer in a domain ${\Omega _d}$ and a source $S$ that may be far from the target. The scatterer is within the meshed area ${\Omega _m}$. (b) From an observer located within $\Gamma$, i.e., in ${\Omega _i}$, the field radiated by the virtual antenna located on $\Gamma$ is the same as those radiated by the remote source $S$. The source $S$ is still present as a “reminiscence” of the problem described in panel (a). In the text, the interior (resp. the exterior) of $\Gamma$ is denoted ${\Omega _i}$ (resp. ${\Omega _e}$).

Download Full Size | PDF

 

Fig. 2. Real parts of the electric field and the Poynting vector. In all sub-figures, the incident field to be simulated is a plane wave coming from the left. It can be seen that the simulated field is only correct inside the virtual antenna and that, outside, this antenna radiates a field that has no other interest than to meet outgoing wave conditions. (a) The virtual antenna is circular. The incident field has a wavelength that is equal to the radius. (b) The incident field has a wavelength that is smaller to the radius. (c) The incident field has a wavelength that is much smaller than the radius. (d) The virtual antenna is rectangular.

Download Full Size | PDF

C. Method of the Virtual Antenna

1. Generalities

The problem of sources has yet to be solved. To do this, all you have to do is change the unknown: switch from the total field to the diffracted field. This has a double advantage. On the one hand the sources are located in the diffracting object; on the other hand the diffracted field contrary to the total field has the good taste to check outgoing wave conditions, which is fundamental for PMLs to play their role. But in the case of nonlinear schemes, the game is not so easy. Because in the linear case, the diffracted field is a solution of the same PDEs as the total field: only the sources change (cf. [3], for example). One quickly becomes convinced that the diffracted field, contrary to the linear case, is not the solution of the same PDEs as the total field, which is very annoying. A subterfuge is then used: the virtual antenna, the principle of which will be now explained.

The idea is summarized in Fig. 1: we wish to simulate the incident field ${{\textbf{E}}^i}$ radiated by a source $S$ (or without any sources when dealing with plane waves, for instance) from a suitably chosen source inside ${\Omega _m}$. To fix the ideas this set ${\Omega _m}$ may be considered as the set to be meshed in the context of the finite element method. As for the source, it is a current located on a closed curve $\Gamma$ (a Jordan’s lace) in 2D (or a closed surface in 3D), itself contained in ${\Omega _m}$ and containing ${\Omega _d}$, which encompasses all the scatterers: the source ${{\textbf{j}}_{|\Gamma}} = {\textbf{j}} {\delta _\Gamma}$ is then represented by a singular distribution whose support is located on $\Gamma$. This curve thus defines an inside ${\Omega _i}$ and an outside ${\Omega _e}$. We therefore need to find the current $j$, which, when chosen correctly, generates the same field ${{\textbf{E}}^i}$ in ${\Omega _i}$ as the field ${{\textbf{E}}^i}$ created by $S$. It remains now to determine ${{\textbf{j}}_{|\Gamma}}$ from both $\Gamma$ and the incident field. Contrary to what one might think at first glance, it is not that difficult. Indeed, it is sufficient to first solve the problem of diffraction by an infinitely conducting obstacle placed in ${\Omega _i}$ illuminated by the incident field ${{\textbf{E}}^i}$. We therefore obtain the total field ${\textbf{E}}$, which allows us to calculate the magnetic field ${\textbf{H}}$. We then calculate the tangential trace of the magnetic field on $\Gamma$, ${\textbf{n}} \times {{\textbf{H}}_{|\Gamma}}$, which gives within a sign the current we are looking for. Indeed, the trick is that inside the infinitely conducting metal the total field vanishes, which means that the diffracted field is equal to the opposite of the incident field. It is worth noting that the details of the implementation can be found in [4].

 

Fig. 3. Virtual current on a circular antenna together with the radiated field. (a) Real part of the virtual current together with the radiated field. (b) Imaginary part of the virtual current together with the radiated field.

Download Full Size | PDF

2. Some Examples

In this section four 2D examples of virtual antennas are given: five circular antennas [Figs. 2(a)–2(c), 3(a), and 3(b)], for which a closed formula can be obtained in the form of a series operating at increasingly higher frequencies, and a square antenna [Fig. 2(d)], which requires an additional numerical step, that of calculating the electromagnetic field radiated by an infinitely conducting square obstacle.

Table 1. With Wavelength in Nanometers, the Dimensionless Refractive Index and the Entrances We Will Be Using for Tensors of Susceptibility, Respectively Dimensionless, in ${\rm{m}}/{\rm{V}}$ and in ${{\rm{m}}^2}/{{\rm{V}}^2}$

View Table | View all tables in this article

D. Some Full Examples

1. Scattering by a Nonlinear Material

In this example, we study a material that is described by the set of nonlinear equations I.27, for a simple cylinder of circular cross section. This circular rod is supposed to be filled with barium borate (BBO) and surrounded with a vacuum. For this material, the parameters given in Table 1 have been used.

In Fig. 4, we see several concentric circles. These represent, with increasing radii, the section of the diffracting object, the circle on which we carry out energy balances (to illustrate the energy transfer developed in part II—they are not detailed here), and the virtual antenna. For simulations, we still need two more circles, marking the inner boundary (we see some of them on the four corners of the images) and the outer end of the PML.

 

Fig. 4. Scattering by an infinite cylinder of BBO, as modeled with system 1. The incident plane waves oscillates at 1064 nm, very close to the cylinder diameter, 1000 nm. The three columns give the three harmonics at a common scale. The amplitude increases from line to line. We successively observe an almost linear regime (${A^i}= 10^6\;{\rm{V}}/{\rm{m}}$), second harmonic generation (${A^i}= 10^9 \;{\rm{V}}/{\rm{m}}$), third harmonic generation (${A^i} = 2 {10^9}\; {\rm{V}}/{\rm{m}}$), and perturbation of the first harmonic with the optical Kerr effect and depletion of the pump beam (${A^i}\geqslant 310^9 \;{\rm{V}}/{\rm{m}}$).

Download Full Size | PDF

2. Electromagnetic Scattering by a Circular Cylinder, in the Order and at Degree Three: Second Generation and Third Harmonics and Kerr-Optics Effect

We are in the situation described by the system I.27, where the frequencies ${\omega _I}$, $2 {\omega _I}$, and $3 {\omega _I}$ appear, namely,

On the other hand, in order not to consider a case that is too complicated and too time and memory consuming, we wished to remain in a 2D scalar case, where the incident field is linearly polarized along the invariance axis of the supposedly infinitely long diffracting object along the $O{x_3}$ axis. To be more precise, the incident field is a plane wave: ${{\textbf{E}}^i} = {A^i}\exp (i\frac{{{\omega _I}}}{c}{x_1}){{\textbf{e}}_3}$. But the game was not over at that point; the material must have characteristics such that the diffracted field itself remains linearly polarized so that the problem is reduced to a scalar problem. It turns out that in order to have a material with the required characteristics, it is sufficient that the $\chi _{(n)}^{{i_1},{i_2}, \cdots ,{i_n}}$ tensors do not have a component that contains 3 unless they only contain 3! For example, $\chi _{(1)}^{i,3}$ vanishes unless $i = 3$. As for the pseudo-tensor ${\chi _{(2)}}$, all the terms $\chi _{(2)}^{i,j,3}$ vanish unless $i = j = 3$, and so forth. In that case, if ${{\textbf{J}}_1} = {J_1} {{\textbf{e}}_3}$, the aforementioned set becomes by writing ${u_p} = {{\textbf{E}}_p} \cdot {{\textbf{e}}_3}$

3. Results and Interpretation

It is sometimes asserted in the literature ([5], for example) that the third harmonic is mainly generated by a cascade effect. When the amplitude of the incident field is of the order of ${A^i} \approx {10^9} \;{\rm{V}}/{\rm{m}}$, we have roughly $||{{\textbf{E}}_1}|| \approx {10^9} \;{\rm{V}}/{\rm{m}}$, $||{{\textbf{E}}_2}|| \approx {10^8} \;{\rm{V}}/{\rm{m}}$, and $||{{\textbf{E}}_3}|| \approx {10^6}\; {\rm{V}}/{\rm{m}}$. Also, taking the approximations for the BBO ${\mu _r} \approx 1$, ${\chi _{(1)}} \approx 1$, ${\chi _{(2)}} \approx \,{10^{- 12}} \;{\rm{m}}/{\rm{V}}$, and ${\chi _{(3)}} \approx {10^{- 23}} \;{{\rm{m}}^2}/{{\rm{V}}^2},$ the orders of magnitude of each term in the equation are obtained for ${{\textbf{E}}_p}$, which are reported in Table 2. In Fig. 4, the real parts of the fields ${{\textbf{E}}_1}$, ${{\textbf{E}}_2}$, and ${{\textbf{E}}_3}$ corresponding to the frequencies ${\omega _I}$, $2{\omega _I}$, and $3{\omega _I}$ are shown. These frequencies correspond to the wavelengths in vacuum at ${\lambda _1} = 1064 \;{\rm{nm}}$ (infrared), ${\lambda _2} = 532\; {\rm{nm}}$ (green visible light), and ${\lambda _3} = 355 \;{\rm{nm}}$ (near ultraviolet), respectively. We have therefore represented ${u_1}$, ${u_2}$, and ${u_3}$ as solutions of the system 1 with a current ${J_1}$ that simulates a plane wave. Figure 4 corresponds to different and increasing incident field amplitudes from ${A^i}= 10^6 \;{\rm{V}}{{\rm{m}}^{- 1}}$ (first line of sub-figures) to ${A^i} = 5 {10^9} \;{\rm{V}}{{\rm{m}}^{- 1}}$ (last line), on the same energy scale. When ${A^i}= 10^6 \;{\rm{V}}{{\rm{m}}^{- 1}}$, the oscillating fields at $2{\omega _i}$ and $3{\omega _i}$ are imperceptible. When ${A^i}= 10^9 \;{\rm{V}}{{\rm{m}}^{- 1}}$, a field begins to appear at $2{\omega _i}$ (i.e., in the green). It is not until the amplitude ${A^i} = 2 {10^9} \;{\rm{V}}{{\rm{m}}^{- 1}}$ that third harmonic generation appears. We can see that the field at ${\omega _I}$ starts to change slightly, but the hypothesis of the nondepletion of the pump wave is still credible. When the amplitude of the incident field further increases, the assumption of nondepletion of the pump wave must be seriously questioned. Oscillations that are typical of double and triple frequencies are clearly visible. From the amplitude ${A^i} = 3 {10^9} \;{\rm{V}}{{\rm{m}}^{- 1}}$, we see increasingly strong effects of the nonlinearity, and the wave at the frequency ${\omega _I}$ (thus oscillating at the pump wave) is increasingly modified. From system (1), this is caused by the two nonlinear terms, $2\langle \langle {u_{- 1}},{u_2}\rangle \rangle$ and $3\langle \langle {u_{- 1}},{u_1},{u_1}\rangle \rangle$. However, according to Table 2, the former is one order of magnitude greater that the latter.

Table 2. Orders of Magnitude, in ${\rm{V}}/{\rm{m}}$, of the Different Contributions ${{\textbf{E}}_p}$ When ${A_i} \approx {10^9}\;{\rm{V}}/{\rm{m}}$

View Table | View all tables in this article

3. MATHEMATICAL SURVIVAL GUIDE

A. Levi-Civita’s Symbol

In this paper, and especially in Section 2.D, we have used Levi-Civita’s symbol. We recall briefly the definition. Let ${{\textbf{e}}_i}$ be the $i$th vector of the standard basis defined as follows: ${{\textbf{e}}_i} = ({\delta _{1,i}}, \cdots ,{\delta _{j,i}}, \cdots ,{\delta _{n,i}}{)^T}$. In other words, every component of ${{\textbf{e}}_i}$ is null except the $i$th, which is equal to 1. Let $I$ be a set of $n$ integer $I: = ({i_1}, \cdots ,{i_n})$; Levi-Civita’s symbol [6] is the application from $I$ towards the set $\{{0, - 1,1} \}$ defined by

We then have, for instance, ${\epsilon _{1,2,3}} = 1$ and ${\epsilon _{0,1,2,3}} = 1$. Finally, for convenience, we need to manipulate Levi-Civita’s symbol with the form $\epsilon _I^J$, where $I$ and $J$ are two sets of integers for applying Einstein’s rule. In that case, it is just a notation

For instance, we have $\epsilon _i^{j,k} = {\epsilon _{i,j,k}}$. Note that with the above convention, we have $\epsilon _I^J = \pm \epsilon _J^I$, depending on the length of $I$ and $J$. Besides, the reader can check the following useful results:

B. Minimum on Orthogonal Matrices

1. Generalities

In the following, the matrices are assumed to be real valued. In that case, a matrix $r$ is said to be orthogonal if its transpose ${r^T}$ is equal to its inverse, namely,

These matrices are of utmost importance because they preserve the scalar product in the sense that they are associated with linear transformations $\rho$ that preserve the scalar product in question. Their importance is such that mathematicians have given a name to the set of all these matrices, namely, $O(n;{\mathbb{R}})$. For all vectors of ${{\mathbb{R}}^n}$ (in practice, $n = 3$), ${\textbf{x}}$ and ${{\textbf{x}}^\prime}$

A matrix can be associated with each linear transformation in a particular ${\cal B}$ base. Let ${R_{|{\cal B}}}$ be the matrix associated with $\rho$ in the base ${\cal B}$, and let ${\underline {\textbf{x}} _{|{\cal B}}}$ (resp. ${\underline {{{\textbf{x}}^\prime}} _{|{\cal B}}}$) be the column vector associated with ${\textbf{x}}$ (resp. ${{\textbf{x}}^\prime}$) in the base ${\cal B}$; the equality written above becomes

If the basis ${\cal B}$ is orthogonal, the scalar product is simply written as a matrix product (raw matrix by a column matrix)

So we have

This equality is well verified since the matrix ${R_{|{\cal B}}}$ is orthogonal. Henceforth the reference to the base will be implicit, and unless otherwise stated we will write $R$ rather than ${R_{|{\cal B}}}$. We will now see that the matrices in $O(n;{\mathbb{R}})$ fall into two categories according to the sign of the determinant. We have indeed

Orthogonal matrices therefore preserve the algebraic volumes, which is not a surprise since they preserve the scalar products and therefore the lengths and the angles between arbitrary vectors. We therefore have this subdivision: matrices with a determinant equal to 1 and those with a determinant equal to ${-}1$. The matrices having a determinant equal to 1 belong to a group contrary to those having a determinant equal to ${-}1$, and the space they are living in has been christened $SO(n;{\mathbb{R}})$ by the mathematicians. In ${{\mathbb{R}}^2}$, it is easy to show that these matrices are written in the following manner:

In ${{\mathbb{R}}^3}$, it can be shown that a matrix $r$ in $SO(3;{\mathbb{R}})$ associated with a linear transformation $\rho$ leaves invariant all vectors colinear to a fixed non-null vector ${\textbf{u}}$ and is reduced to a rotation in ${{\mathbb{R}}^2}$ for all vectors orthogonal to ${\textbf{u}}$ [7]. The vector ${\textbf{u}}$ is then a vector characterizing the axis of rotation. Moreover, if ${\textbf{u}}$ is a unit norm vector ${\textbf{u}} = {({{u_1},{u_2},{u_3}})^T}$ expressed in an appropriate basis ${\cal B}$ and $R({\textbf{u}},\theta)$ the rotation matrix of axis $\Delta$, characterized by the vector ${\textbf{u}}$, and the angle $\theta$, the so-called Rodriguès’ formula can be shown:

We give, for example, the three rotation matrices ${R_1}(\theta)$, ${R_2}(\theta)$, and ${R_3}(\theta)$ about the three axes ${x_1}$, ${x_2}$, and ${x_3}$:

2. Matrix Writing versus Index Writing

Let us consider two linear applications ${\cal A}$ and ${\cal B}$. With the tensorial notation, the associated matrix $a$ (resp. $b$) is characterized by its components denoted $a_j^i$ (resp. $b_j^i$) instead of ${a_{i,j}}$ (resp. ${b_{i,j}}$). The matrix $c$ defined as the product $c = a b$ is such that

With this notation, how can we rewrite the relation II.31, namely,

The term $\sum\nolimits_{i_1^\prime = 1}^n {r_{{i_1},i_1^\prime}} {\chi _{i_1^\prime ,i_2^\prime}}$ corresponds to the element ${i_1},i_2^\prime $ of the product $r \chi$, denoted by ${(r \chi)_{{i_1},i_2^\prime}}$. So we have

This time, the dummy index appears twice on the column index: the right-hand term is thus not ${(r \chi r)_{{i_1},{i_2}}}$. So, in order to introduce the matrix product, it is necessary beforehand to invert the indices of rows and columns: for that, we have the transposition operator $T$ in the hands

In other words, Eq. (5) reads

C. Properties of the Forms $\langle \langle \cdots \rangle \rangle$

1. On the Tensor Nature of the Forms $\langle \langle \cdots \rangle \rangle$

We recall the definition of this multilinear form

2. Algebraic Properties

3. Diagrammatic Properties

In this subsection, we exemplify the intrinsic permutation symmetry by a diagram of the form $\langle \langle {{\textbf{E}}_2},{{\textbf{E}}_1},{{\textbf{E}}_{- 1}}\rangle \rangle$, which involves the tensor ${\chi _{(3)}}(2 {\omega _I};2 {\omega _I},{\omega _I}, - {\omega _I})$. In Fig. 5, six different processes involving the tensor ${\chi _{(3)}}$ at the frequency $2{\omega _I}$ are drawn. These six different terms are indistinguishable by the virtue of the intrinsic permutation symmetry: they mimic three-photon quantum interactions involving creation and annihilation.

 

Fig. 5. Intrinsic permutation symmetry shows the equalities of the following six forms: $\langle \langle {{\textbf{E}}_2},{{\textbf{E}}_1},{{\textbf{E}}_{- 1}}\rangle \rangle$, $\langle \langle {{\textbf{E}}_2},{{\textbf{E}}_{- 1}},{{\textbf{E}}_1}\rangle \rangle$, $\langle \langle {{\textbf{E}}_1},{{\textbf{E}}_2},{{\textbf{E}}_{- 1}}\rangle \rangle$, $\langle \langle {{\textbf{E}}_1},{{\textbf{E}}_{- 1}},{{\textbf{E}}_2}\rangle \rangle$, $\langle \langle {{\textbf{E}}_{- 1}},{{\textbf{E}}_1},{{\textbf{E}}_2}\rangle \rangle$, and $\langle \langle {{\textbf{E}}_{- 1}},{{\textbf{E}}_2},{{\textbf{E}}_1}\rangle \rangle$.

Download Full Size | PDF

D. Transposition Operators for the Tensors ${\chi _{(n)}}$

1. With a (1, 0)-Tensor

For any couples of vectors ${{\textbf{v}}_{\rm{a}}}$ and ${{\textbf{v}}_{\rm{b}}}$,

2. With a (2, 0)-Tensor

For any set of vectors ${{\textbf{v}}_{\rm{a}}}$, ${{\textbf{v}}_{\rm{b}}}$, and ${{\textbf{v}}_{\rm{c}}}$,

4. SOME WORDS ABOUT THE PROGRAM MATHEMATICA

In this program [8] it is assumed that the material at stake is characterized by two symmetries given by two rotation matrices, ${r_1}({\theta _1},{{\textbf{v}}_1})$ and ${r_2}({\theta _2},{{\textbf{v}}_2})$. In the following example the two rotation matrices are RR1 and RR2. RR1 (resp. RR2) corresponds to a rotation about the $x$ axis (resp. $z$ axis) with an angle of $\frac{\pi}{2}$ (resp. $\pi$):

Equation II.30 is translated into the Mathematica language for the tensor ${\chi _{(3)}}$ and the rotation matrix ${r_1}$ (RR1).

Eq. II.30 is re-indexed thanks to the instruction Flatten.

The system amounts then to addressing an homogeneous linear system that is “solved” by the NSolve instruction.

The program finally outputs the null terms and the independent terms as well as the dependencies of the terms between them as done in part II, Section 5.D.

5. CONCLUDING REMARKS

The world of nonlinear electromagnetism is an extremely rich field both from a theoretical point of view and from the point of view of concrete applications. Much has been overlooked. For example, it is to be deplored that the part devoted to numerical calculation is reduced to a minimum. I am aware of this, and it is possible that in the future we could make a tutorial only on this fundamental problem. This will be the occasion to illustrate the theoretical parts of the two previous parts and to show the immense variety of nonlinear effects. It is perhaps regrettable that there is no reference to quantum physics or that there is no interest in propagation conditions with the hypotheses of slowly varying envelopes, or solitons. Choices had to be made with the revolution in nonlinear nano-optics in mind. I hope this goal has been achieved.