Abstract

In presence of weakly nonlinear media, it is tempting to neglect pump wave depletion when calculating the intensity of the various generated nonlinear components. However, even in the case of very weak pump wave depletion conditions, an optical structure that allows multipass optical paths, such as high finesse multilayer microcavities, could lead to significant cumulative pump wave depletion. In such conditions, neglecting pump depletion might lead to large computational errors. A matrix formalism devoted to such pump depletion in planar layered nonlinear structures without resort to the “bound” and “free” waves concept is described. A general approach that makes use of “intrinsic” and “extrinsic” waves concept is given, through a slight modification of the canonical propagation matrix. The theoretical results show that even in the case of very weak pump depletion conditions, the cumulative effect due to confinement actually leads to very significant effects. It turns out that taking into account the pump depletion is mandatory for numerous experimental conditions. This matrix formalism applies to absorbing media, and is extensible to include the case of anisotropic layers and cascading effects.

1. Introduction

From both classical electromagnetism and quantum electrodynamics point of view, there is a constantly growing interest with respect to nonlinear optics in confined structures [110]. In the last two decades, a large amount of theoretical and experimental efforts have been devoted to nonlinear processes in optical microcavities, including layered structures [1122]. Notable attention was given to harmonic generation. Besides a better understanding of basic phenomena in highly confined media, interest is also linked to the possibility of conceiving new technological devices. Indeed, it is now well established that in highly confined structures many optical processes are strongly modified [2329]. From the point of view of quantum electrodynamics, the modification of vacuum field fluctuations and the concomitant modulation of the density of states are responsible for such effects [30]. The resultant enhancement and inhibition for specific optical processes open doors toward new possibilities devoted to the manipulation of electromagnetic radiation. An important applications field is integrated optics and significant innovations are anticipated.

Because high-confinement structures imply many roundtrips of optical waves, even small pump depletion could cause a significant cumulative reduction in the pump wave amplitude. It is an important goal of the present paper to address this issue. Among other possible formal approaches, including purely numerical ones, here we focus on the transfer matrix method and show how a generalized version of the so-called “propagation matrix” leads to a convenient matrix formalism. Indeed, the transfer matrix method is widely used to calculate the optical properties of layered structures. This formalism automatically takes into account all multiple reflections for all waves components (pump and nonlinear waves) without explicit enumeration and cumbersome summations. Moreover, this formalism is well suited to computer calculations. This formalism could be especially useful in the context of multilayer microcavities. In brief, this paper describes an alternative approach and we believe that the picture and formalism linked to this formalism are both simple and allow further generalization.

2. System Description

Figure 1 shows a general planar multilayered structure made of N stratums (including the substrate and the superstrate). The layers are of infinite extension in the “x-y” plane. The incidence medium is labeled 0. The layer in concern is medium l (which should go over 0 to N in order to account for all nonlinear layers). The transmission medium is 𝑁+1. It is assumed that both media 0 and 𝑁+1 are semi-infinite and transparent. All layers are nonmagnetic (𝜇𝜇𝑜), can be nonlinear and anisotropic. However, they are homogeneous. For every layer, including the substrate and the superstrate, the coherence length of the incident radiation is assumed much larger than its optical thickness. The substrate and superstrate can thus be accounted by the matrices described below. If the substrate or the superstrate have an optical thickness much larger than the incident radiation coherence length and are made of linear material, one takes account of their presence through the reflectance of their respective interface with surrounding media. The intermediate context of layers thickness similar to the incident radiation coherence length needs specific approach. It is omitted here. As for any layer, when the substrate or the superstrate are absorbing, one includes the exponential decrease absorption law (Beer-Lambert-Bouguer).

The incident radiation is a plane wave. Section 4 describes the formal details. In a few words, each layer l has a thickness 𝑠𝑙 and a complex refractive index 𝜂𝑙𝑛𝑙𝑖𝜅𝑙 (The layer is absorbing when 𝜅𝑙 is positive). This choice presumes that a forward traveling wave is described by an exponential phase factor of the form (𝜔𝑡𝑘𝑟). Care must be practiced with various possible conventions [31]. The electric field of the light wave components is 𝐸. The input field, 𝐸+in, is associated to the incident (external) pump wave while the output fields, 𝐸out and 𝐸+out, are linked to the nonlinear outgoing waves emerging from the structure. The outgoing components 𝐸out and 𝐸+out are propagating in the backward (incidence side) and forward (transmission side) direction respectively.

3. Formalism Overview

Our approach is based on a matrix formalism [3235]. In brief, column vectors describe the fields everywhere along the “𝑧” axis while square “transmission matrices” describes the system. Any field vectors pairs are linked through a square transmission matrix. In presence of isotropic media, the field vectors have two rows, and are linked through 2×2 matrices. In presence of anisotropic media, the field distribution within each layer can be expressed as a sum of four partial waves (eigenwaves) [3446]. Then, one deals with four-row field vectors and 4×4 matrices. Fortunately, in numerous practical applications we deal with anisotropic transparent layer having one axis normal to the layers' plane. Then, working with plane polarized pump waves, orienting properly its polarization plane, and restraining the incidence to the normal angle in case of biaxial material, the isotropic matrix formalism is suitable. Working with effective nonlinear coefficient is another possibility when dealing with anisotropic nonlinear media [47].

The point of view underlying the current formal approach can be summarized simply with the concept of “intrinsic” and “extrinsic” waves. This concept was already put forward with success [15] (the concept of “intrinsic” and “extrinsic” waves, equivalent of the one of “self" and “add" waves found in [15], were coined without knowledge of [15], during a study on microcavities in France Telecom's laboratories, in 1995). In brief, a contrario to the point of view based on the concept of “bound” and “free” waves [4850], here there is no need for a source vector and, in addition, we use a modified “propagation matrix” (details specified in Section 5).

For the purpose of the present discussion, we are dealing with two semi-infinite media where the intrinsic wave does not suffer any reflection. From the intrinsic and extrinsic waves point of view, a nonlinear wave generated within a bulk medium is depicted as a unique wave propagating along with the pump wave and having amplitude that varies in space. Because this wave does not originate from outside the bulk medium nor depends on external factors such as the presence of interfaces, this wave is thus said to be “intrinsic”. Consequently, in context of harmonic generation the nonlinear wave is not given by the superposition of constant amplitude “bound” and “free” waves from which interference governs the spatial amplitude modulation. As shown in Figure 2, an intrinsic nonlinear wave propagating in a medium, 𝑙, originates from within that medium. A contrario, an “extrinsic” wave propagating in a medium, l, originates from an interface or an external stratum. The most important feature of an extrinsic wave is that it travels in space with constant amplitude (except when linear absorption occurs).

Figure 2 also depicts the spatial evolution of the pump and nonlinear waves as they go deeper into a nonlinear medium. Both waves have variable amplitude. The (internal) pump field is the one inside the nonlinear medium, not the one incident on the medium 𝑙.

Our approach involves each nonlinear layer separately. For each such layer, the specific goal is to calculate the field amplitude of the corresponding intrinsic (nonlinear) and pump waves versus depth inside the medium, at every wavelength. From now on, this pair of solutions is the so-called “non-linear-layer solution.” Once that solution is known for the layer 𝑙, the matrix formalism allows to easily relates the nonlinear fields inside the structure to its outside. Of course, the use of the Manley-Rowe relation is invaluable to compute pump depletion.

Wave mixing is also important in nonlinear optics. Here we wish to make distinction between mixing of two, or more, incident pump waves and mixing of generated nonlinear waves with pump wave or other nonlinear waves (cascading). We refer to the former as “primary mixing” while the later is designated as “secondary mixing”. In the context of primary wave mixing, the overall pump field in the layer 𝑙 is given by the superposition of the incident pump waves. On the other hand, the intrinsic wave is alternatively one of the nonlinear components making the global “nonlinear layer solution”. Thus, in case of two pump waves, 𝜔1 and 𝜔2, in a second-order nonlinear medium, one computes every component 𝜔1+𝜔2, 𝜔1𝜔2, 2𝜔1, and 2𝜔2 separately. Even if only one of these nonlinear components is under concern, one has to compute all components if accurate pump depletion is required.

In presence of many nonlinear layers, when secondary mixing is neglected, the overall output is obtained by running 𝑙 through all layers and summing all the outgoing nonlinear waves (including the phase) propagating in the transmission side. From now on, this is the so-called “process.” Provided that secondary mixing entails negligible depletion, the matrix method presented below can account for cascading. Indeed the secondary mixing waves are simply seen as usual pump waves.

In presence of weak nonlinearity, it would be tempting to neglect pump depletion. However, in case of large number of nonlinear layers, the cumulative effect of even weak pump depletion in each layer could lead to large computational error. Thus, in order to take account of pump depletion that might occur in nonlinear layers, the pump field distribution within the structure is treated as a whole (self-consistency).

The linear and nonlinear properties of every layer are supposed known. This adds to material's axis orientations when anisotropy can't be neglected. It will soon become obvious to the reader that the extension of our approach to anisotropic layers is straightforward, though more algebraically cumbersome. The difficulties are only computational, not conceptual. However, it adds very little to the comprehension of the current approach to keep track of every field components and advance with a complete formal treatment of N anisotropic layers. Because the isotropic and anisotropic formalisms are isomorphic, the generalization to the case of anisotropic layers is left as a forthcoming paper. Specific formalisms are better adapted to this extension than others [4146]. Then, for the sake of simplicity the formalism put forward in this paper is devoted to isotropic layers. It provides the basic description of our formalism. Subsidiary parameters such as, inhomogeneity, interface roughness, coherence length, beam size, anisotropy and so forth are omitted. This enables us to focus on the basic ideas without the added algebraic complications.

In presence of isotropic layers only, an arbitrary pump wave has to be separated into two orthogonal eigenmodes (“s” and “p” polarization). The resultant, overall outgoing waves are calculated separately and combined at the end of the “process”. Because of their more general structure, no such double treatment is required when anisotropic media matrices are in use.

4. Definitions, Conventions, and Notation System

For sake of simplicity, and in order to keep track of all phase shifts, the complex amplitude (phasor) notation is used to fully describe the electric field. Because it doesn’t add any supplementary information, the use of complex conjugate when dealing with real quantities is redundant and not systematically shown in equations.

In order to avoid confusion, the formalism makes use of a subscripts-superscripts notation system, for both coordinates and fields. The superscripts are sometime bracketed to avoid mixup with mathematical exponents. As shown in Figure 1, the incident waves propagate along the “z” axis, forward from the right to the left. The superscripts “+” and “−” are used to indicate the propagation direction when applied to E. These superscripts are also applied to z with a slightly different meaning. A self-explained example is shown in Figure 1. It designates the amplitude () of the intrinsic (I), nonlinear (N), and electric field (E) of a wave located on the backward side (−) of the coordinate 𝑧𝑙 (the interface between layers m and l) and propagating into layer l in the forward (+) direction.

For the purpose of easy reading of some important equations, we sometimes write the abbreviated form of selected specific fields. According to the notation system, the abbreviated forms for the incident pump wave field (𝐸+in) and the generated outgoing waves (𝐸out,𝐸+out) are defined by 𝐸+in𝐸(+,P)0𝑧0,𝐸(1)out𝐸(,N)0𝑧0,𝐸(2)+out𝐸(+,N)𝑁+1𝑧+𝑁.(3) In (1), we keep the superscript + in order to distinguish from an ongoing pump wave eventually incident on the structure from the transmission side. Indeed, the picture described below remains valid. Nevertheless, here we set the negative input component to 0 (𝐸in=0). We define the abbreviated forms for intrinsic (𝐸+int and 𝐸int) and pump (𝐸+pump and 𝐸pump) fields, shown in Figure 3, as follows 𝐸+int𝐸(+,N,I)𝑙𝑧𝑙,𝐸(4)int=𝐸(,N,I)𝑙𝑧+𝑘,𝐸(5)+pump=𝐸(+,P)𝑙𝑧+𝑘,𝐸(6)pump=𝐸(,P)𝑙𝑧𝑙.(7) Finally, our convention for isotropic media is that the forward traveling waves amplitude figures in the first row of the column field vectors.

5. Basic Matrices

In regard to the current formalism, four basic system matrices are required (“coupling,” “continuity,” “propagation,” and “characteristic” matrix). Care must be practiced with various possible conventions [31] when writing matrices.

5.1. Coupling Matrix

On both sides of an interface, say interface l/k, the field vectors are related through a “coupling matrix," 𝐶𝑙𝑘, such that (incidence medium is labeled by the rightmost index) 𝐸+𝑙𝑧+𝑘𝐸𝑙𝑧+𝑘=𝐶𝑙𝑘𝐸+𝑘𝑧𝑘𝐸𝑘𝑧𝑘.(8) The form of the coupling matrix does not depend on the polarization type and is given by 𝐶𝑙𝑘1𝜏𝑘𝑙1𝜌𝑘𝑙𝜌𝑘𝑙1.(9) However, the Fresnel coefficients 𝜌 and 𝜏 (reflection and transmission, resp.) depend on the polarization type. Here, we adhere to the normal (Fresnel) convention [31]. We refer the reader to literature for their explicit expression.

5.2. Continuity Matrix

On both sides of an interface, say interface 𝑙/𝑘, the field vectors are also related through “continuity matrices” (This matrix is also known as “admittance matrix" and “dynamical matrix"), 𝐷𝑙 and 𝐷𝑘, such that 𝐷𝑙𝐸+𝑙𝑧+𝑘𝐸𝑙𝑧+𝑘=𝐷𝑘𝐸+𝑘𝑧𝑘𝐸𝑘𝑧𝑘.(10) These matrices originate from the tangential continuity condition for electric and magnetic fields at an interface (in absence of interfacial current density). They can be used as an alternative definition to the coupling matrix. Indeed, from (8) and (10), we conclude that 𝐶𝑙𝑘=𝐷𝑙1𝐷𝑘.(11) The continuity matrix depends on polarization. It is given by (Y is the admittance and 𝜃 the propagation angle) 𝐷𝑙=11𝑌𝑙𝜃cos𝑙𝑌𝑙𝜃cos𝑙(12) for s type polarization, and by 𝐷𝑙=𝜃cos𝑙𝜃cos𝑙𝑌𝑙𝑌𝑙(13) for p type polarization.

5.3. Propagation Matrix

On both intrinsic sides of a single layer, say layer l, the field vectors are related through a “propagation matrix" (we prefer “propagation matrix" to the denomination “phase matrix" because that matrix contains information on amplitude in addition to translational phase shift), 𝑃𝑙, such that 𝐸+𝑙𝑧𝑙𝐸𝑙𝑧𝑙=𝑃𝑙𝐸+𝑙𝑧+𝑘𝐸𝑙𝑧+𝑘.(14) The propagation matrix does not depend on polarization and is given by 𝑃𝑙=𝐴+𝑙exp𝑖𝜑𝑙00𝐴𝑙exp+𝑖𝜑𝑙.(15) This matrix also agrees with the convention that a forward traveling wave is described by an exponential phase factor of the form (𝜔𝑡𝑘𝑟). The reader would notice that, the propagation matrix defined here has a slightly different form than usual. Indeed, we introduced the “nonlinear amplitude coefficients” 𝐴+𝑙 and 𝐴𝑙. This allows accounting for the amplitude variation other than linear absorption (Amplification in case of lasing medium). Indeed, linear absorption is fully described by the complex phase 𝜑𝑙. In context of linear optics the amplitude coefficients both reduce to 1. Formally, we define the amplitude coefficients for both internal (not the incident or external) pump waves as, using (6) as follows: 𝐴+𝑙𝐸+𝑙𝑧𝑙𝐸+𝑙𝑧+𝑘=𝐸+𝑙𝑧𝑙𝐸+pump=𝐴+𝑙𝐸+pump(16) and, using (7) as 𝐴𝑙𝐸𝑙𝑧𝑙𝐸𝑙𝑧+𝑘=𝐸pump𝐸𝑙𝑧+𝑘=𝐴𝑙𝐸pump.(17) In presence of depletion, the pump waves are such that |𝐴+𝑙|<1 and |𝐴𝑙|>1. Most of the formal difficulties come from the nonlinear behavior of the optical process. Indeed, nonlinearity imposes that 𝐴+𝑙 and 𝐴𝑙 depend on the field amplitudes 𝐸+pump and 𝐸pump respectively. We discuss the formal implications below. For the moment, we assume that the amplitude coefficients are known. As stated in Section 3, the most important feature of the extrinsic waves is that they are described by propagation matrices having amplitude coefficients 𝐴+𝑙 and 𝐴𝑙 set to 1.

5.4. Charateristic Matrix

With the help of continuity and propagation matrices, one defines the “characteristic matrix" (“Characteristic matrices” are also known as “Abelès matrices") as 𝑀𝑙𝐷𝑙𝑃𝑙𝐷𝑙1.(18) The characteristic matrix contains all the necessary information about the layer l.

5.5. Transmission Matrix

Writing an ordered product of coupling and propagation matrices, we can establish a connection anywhere along the 𝑧-axis between any field vectors pairs, say the fields within the layers 𝑣 and 𝑢 as follows: 𝐸+𝑣𝑧+𝑣1𝐸𝑣𝑧+𝑣1=𝐶𝑣(𝑣1)𝑃𝑣1𝐶(𝑣1)(𝑣2)𝑃𝑣2𝑃𝑢+2𝐶(𝑢+2)(𝑢+1)𝑃(𝑢+1)𝐶(𝑢+1)𝑢𝐸+𝑢𝑧𝑢𝐸𝑢𝑧𝑢.(19) Even in presence of many nonlinear layers causing pump depletion, the formalism remains well founded if one uses the propagation matrices through (15). Introducing the characteristic matrices, we reduce the number of required matrices by a factor of two. Making use of (11), (18) and (19) we easily get 𝐸+𝑣𝑧+𝑣1𝐸𝑣𝑧+𝑣1=𝐷𝑣1𝑀𝑣1𝑀𝑣2𝑀𝑈+2𝑀𝑢+1𝐷𝑢𝐸+𝑢𝑧𝑢𝐸𝑢𝑧𝑢.(20) This allows defining the system transmission matrix as 𝑇𝑣𝑢=𝐷𝑣1𝑀𝑣1𝑀𝑣2𝑀𝑢2𝑀𝑢+1𝐷𝑢.(21) Equation (20) then takes the abbreviated form 𝐸+𝑣𝑧+𝑣1𝐸𝑣𝑧+𝑣1=𝑇𝑣𝑢𝐸+𝑢𝑧𝑢𝐸𝑢𝑧𝑢.(22)

6. Pump and Nonlinear Waves in Layered Structures

In the act of a starting point we assume that the whole structure contains only one (isotropic) nonlinear layer, 𝑙. From (22) we connect the field of the nonlinear wave emerging from the structure in the transmission side to the layer 𝑙 through 𝐸(+,N)𝑁+1𝑧+𝑁0=𝑇(N)(𝑁+1)𝑙𝐸(+,N,I)𝑙𝑧𝑙+𝐸(+,N,E)𝑙𝑧𝑙0+𝐸(,N,E)𝑙𝑧𝑙(23) or equivalently as 𝐸(+,N,E)𝑙𝑧𝑙𝐸(,N,E)𝑙𝑧𝑙=𝑇1(N)(𝑁+1)𝑙𝐸(+,N)𝑁+1𝑧+𝑁0𝐸(+,N,I)𝑙𝑧𝑙0.(24) On the other hand, we connect the field of the nonlinear wave emerging from the structure in the incidence side to the layer 𝑙 as 𝐸(+,N,E)𝑙𝑧+𝑘𝐸(,N,E)𝑙𝑧+𝑘=𝑇(N)𝑙00𝐸(,N)0𝑧00𝐸(,N,I)𝑙𝑧+𝑘.(25) Equations (24) and (25) can now be connected via the propagation matrix 𝑃(N,E)𝑙 such as 𝐸(+,N,E)𝑙𝑧𝑙𝐸(,N,E)𝑙𝑧𝑙=𝑃(N,E)𝑙𝐸(+,N,E)𝑙𝑧+𝑘𝐸(,N,E)𝑙𝑧+𝑘.(26) In (26) we connect extrinsic waves. As stated above, the extrinsic waves are described by propagation matrices having amplitude coefficients 𝐴+𝑙 and 𝐴𝑙 set to 1 (superscript E is used as a reminder). This connection leads to the following: 𝐸(+,N)𝑁+1𝑧+𝑁0𝑇(N)(𝑁+1)00𝐸0(,𝑁)𝑧0=𝑇(N)(𝑁+1)𝑙𝐸(+,N,I)𝑙𝑧𝑙0𝑇(N)(𝑁+1)𝑙𝑃(N,E)𝑙0𝐸(,N,I)𝑙𝑧+𝑘.(27) Making use of (2), (3), (4), and (5), we abbreviate this equation for easier reading as 𝐸+out0𝑇0𝐸out=𝑇𝑙𝐸+int0Π𝑙0𝐸int,(28) where 𝑇𝑇(N)(𝑁+1)0,𝑇𝑙𝑇(N)(𝑁+1)𝑙,Π𝑙𝑇(N)(𝑁+1)𝑙𝑃(N,E)𝑙.(29) Equation (28) set up an explicit link between the pump wave incident on the structure and the outgoing nonlinear components. Indeed, the field amplitudes 𝐸+int and 𝐸int are determined from the “nonlinear-layer solution,” and are known functions. From (28) we get (the matrix element indexes are bracketed to avoid confusion with layers' labels) 𝐸out𝑇=𝑙[]21𝐸+intΠ𝑙[]22𝐸int𝑇[]22,𝐸+out𝑇=[12]𝑇[22]𝑇𝑙[21]𝑇𝑙[11]𝐸+int+𝑇[]12𝑇[]22Π𝑙[]22Π𝑙[]12𝐸int.(30) As illustrated in Figure 3, the process of calculating 𝐸+int and 𝐸int requires the amplitude of the pump waves 𝐸+pump and 𝐸pump given by (6) and (7) respectively. They are determined from (𝜌 is the amplitude reflection coefficient of the whole structure at pump wavelength) 𝐸+pump𝐸(,P)𝑙𝑧+𝑘=𝑇(P)𝑙01𝜌𝐸+in𝐸,(31)(+,P)𝑙𝑧𝑙𝐸pump=𝑃(P)𝑙𝑇(P)𝑙01𝜌𝐸+in=Π(P)𝑙1𝜌𝐸+in.(32) Here, 𝜌 is given by 𝐸𝜌(,P)0𝐸(+,P)0𝑇=(P)[]21𝑇(P)[]22,(33) with 𝑇(P) now defined by 𝑇(P)𝑇(P)(𝑁+1)0.(34) In presence of pump depletion, we make use of the generalized propagation matrix introduced in (15). The amplitude coefficients 𝐴+𝑙 and 𝐴𝑙 are determined from the “nonlinear layer solution” for 𝑙 and are known. One understands that in presence of pump depletion, the amplitude coefficients 𝐴+𝑙 and 𝐴𝑙 cannot be set to 1. Unfortunately, from (6), (7), (15), (16), (17), (18), (31), and (32) one realizes that the fields 𝐸+pump and 𝐸pump are themselves required for their own determination. Indeed, the pump fields are required in the computation of 𝑃(P)𝑙, and therefore also for the computation of 𝜌. This self-consistency requirement leads to two complex implicit equations for 𝐸+pump and 𝐸pump. More specifically, writing (31) and (32) in the following compact forms 𝐸+pump𝐸(,P)𝑙𝑧+𝑘=𝑓1𝑓2𝐸+in,𝐸(+,P)𝑙𝑧𝑙𝐸pump=𝑔1𝑔2𝐸+in,(35) and the implicit equations are (here pump depletion occurs only in one layer) 𝐸+pump=𝑓1𝐸+pump,𝐸pump𝐸+in,𝐸pump=𝑔2𝐸+pump,𝐸pump𝐸+in.(36) In (36) we explicitly introduced the fields as variable for clarity. From these equations, we now introduce the “deviation functions” Δ+Δ+re+𝑖Δ+im and ΔΔre+𝑖Δim as Δ+𝑓1𝐸+pump,𝐸pump𝐸+in𝐸+pump=Δ+𝐸+pump,𝐸pump,Δ𝑔2𝐸+pump,𝐸pump𝐸+in𝐸pump=Δ𝐸+pump,𝐸pump.(37) In general, these functions do not have simple analytic representation. Finding analytic solution to this system could be difficult. However, (37) are easily computable once the amplitude coefficients 𝐴+𝑙 and 𝐴𝑙 are known. Because the pump fields are complex numbers, each deviation function depends on four parameters (in addition to 𝐸+in). From a geometrical point of view, the deviation functions describe two complex (Δ+ and Δ), or evenly four real (Δ+re, Δ+im, Δre and Δim), hypersurfaces in a five dimensional space. The solution for the pump fields 𝐸+pump and 𝐸pump is clearly given at the crossing point of the hypersurfaces in the hyperplane Δ=0. This procedure is equivalent to solving a system of four-variable nonlinear algebraic equations. It is also analogous to load-line analysis practiced in electronics when dealing with nonlinear electronic devices and searching for the operating (quiescent) point. Thus, one solves self-consistency approximately through numerical calculations. However, the accuracy level can be made arbitrary small. Developing algorithms that locate the crossing point is not a difficult task, especially if one puts forward the knowledge one has on the gradient of the hypersurfaces. The solution can be granted a posteriori by the application of conservation laws such as field continuity at all interfaces and/or energy balance accounting for linear and nonlinear waves (even in presence of absorption, if known). In case of many depleting layers, the approach described here still applies. One has to add four more dimensions to the hyperspace each time one includes a depleting layer. As an alternative, an iterative algorithm applied to the pump fields could also help. Finally, the above algorithms extend to the context of primary wave mixing. One adds four dimensions to the hyperspace for each pump wave. Nonetheless, it is not the purpose of this paper to further scrutinize the aspects of pump fields' self-consistency.

7. Testing the Formalism: Second Harmonic Generation

7.1. Structure Description

The cavity's thickness is defined by the thickness of the confined nonlinear layer. The microcavity involves two identical mirrors with one nonlinear layer imbedded in-between. Each mirror consists of 23×𝜆/4 dielectric layers. The low index 𝜆/4 layers are 205.4 nm thick SiO2 (𝑛=1,448) and the high-index layers are 151.7 nm thick ZrO2 (𝑛=1,961). The layers in contact with the nonlinear layer are both SiO2. All dielectric layers are assumed transparent, isotropic and achromatic. The cavity's finesse is close to 400. The incoming pump wave amplitude is 2,0×106 V/m. This corresponds to a fairly small pump intensity of 5,31 GW/m2 (531 kW/cm2, this intensity is approximately obtained with a 1 μJ, 50 ps duration pulse focalized on a 200 μm diameter spot.). The nonlinear medium is presumed to be 2-methyl-4-nitroaniline (MNA), a highly nonlinear material [51, 52], known as a “yellow material”. The only important nonlinear axis is along the cleavage plane. Fortunately, this axis is parallel to the layers’ plane [53]. Working with linearly polarized pump wave having their plane aligned along the nonlinear axis, an isotropic medium formalism is adequate. The incidence angle is set normal. The nonlinear layer is assumed transparent, although chromatic dispersion is included [54]. The second order susceptibility is 𝜒(2)=225 pm/V [55]. Finally, the intensity is obtained from 1𝐼=||𝐸||2𝑍2.(38)

7.2. Comparison of Conceptual Approaches

In order to test the current matrix formalism, we compare calculations, based on the intrinsic and extrinsic waves concept, with previous one [29], based on the concept of bound and free waves (Bethune's formalism [11]). Of course, because only the current formalism takes into account the pump depletion, for the purpose of comparison we neglected pump depletion. In [29], our aim was to calculate the outgoing second harmonic generation (SHG) from a layered microcavity embedding only one second-order nonlinear layer. The calculations are performed for a “thick” microcavity, that is a “30×𝜆/2 microcavity” (𝜆=1,19μm). Figure 4(a) shows calculation based on bound and free waves concept while Figure 4(b) shows calculation based on the current matrix formalism, founded on the concept of intrinsic and extrinsic nonlinear waves. The agreement appears very satisfying. Only minute differences (due to graphical resolution) can be seen. This comparison is an indication that both conceptual approaches are well founded.

Figure 4 does not show however that the intensity of modes is strongly wavelength dependent. Indeed, the intensity of the sharp peaks between 1,05 and 1,30 μm goes far beyond the graph limit of 4000 W/m2. This modulation and “pseudomodes” very close to the edges of the mirror's bandgap below 1,1 μm and above 1,3 μm are discussed elsewhere [29].

7.3. Nonlinear Layer Solutions

To determine the nonlinear-layer solutions, we solve the following one-dimensional nonlinear wave equation (here it is implicit that E applies to the nonlinear wave and v is its phase velocity): 𝜕2𝑧1𝐸(𝑧,𝑡)𝑣2𝜕2𝑡𝐸(𝑧,𝑡)=𝑆(𝑧,𝑡).(39) This equation is readily obtained from basic electromagnetic theory [47]. We assume that the source of nonlinear component is only the second-order nonlinear susceptibility. The nonlinear source function, 𝑆(𝑧,𝑡), is then given by (the subscript stands for “pump") 𝑆(𝑧,𝑡)=𝜀𝑜𝜇𝑜𝜒(2)𝐸2P(𝑧)𝜕2𝑡𝜔exp𝑖2P𝑡𝑘P𝑧.(40) As above, the permeability is set to 1. We search for solutions in the form of a modulated amplitude plane wave 𝐸𝐸[𝑖]=𝐸𝑖(𝑧,𝑡)=(𝑧)exp(𝜔𝑡𝑘𝑧)(𝑧)exp2𝜔P.𝑡𝑘𝑧(41) The occurrence of phase mismatch is formally described through Δ𝑘2𝑘p𝑘=4𝜋𝑛𝑛p𝜆p.(42) For convenience, we investigated four possibilities; (1) nonphase matching SHG without pump depletion (NMND), (2) nonphase matching SHG with pump depletion (NMWD), (3) phase matching SHG without pump depletion (PMND), (4) phase matching SHG with pump depletion (PMWD). The corresponding “nonlinear-layer solutions” are already known (can be found in any good book on nonlinear optics) and summarized in Table 1. They assume slowly varying amplitude approximation (the context of important amplitude variation is easier solved through numerical calculations.) and no absorption at pump and harmonic waves. Nonetheless, absorption at pump wavelength can be included in the calculation of pump fields 𝐸+pump and 𝐸pump given by (31) and (32), respectively.

In order to complete the pertaining “nonlinear-layer solution,” the pump wave amplitude is found through the use of Manley-Rowe relation expressed in the following form (by definition the phase of pump wave is set to 0) ||𝐸P||(𝑧)2=𝐸2P𝑛(0)𝑛P||||𝐸(𝑧)2.(43) We note that from (43), (* *), and (* * * *) in Table 1 one sees the implicit dependence of the amplitude coefficients on the pump fields.

From (43) and (* * * *) in Table 1, one determines that pump wave suffers an amplitude decline of about 15% when 12𝜋𝑛P𝜆P𝜒(2)𝐸P1(0)𝑧2.(44) From (38), this imposes a limit on pump wave intensity given by 𝐼lim1(15%)2𝑍𝑜||||1𝑛4𝜋P𝜆P𝑠𝑙1𝜒(2)||||2=132𝜋𝑍𝑜𝑛P𝜆P𝑠𝑙1𝜒(2)2.(45) For a 1 μm thick MNA pumped at 𝜆=1μm, this limit is about 1,64 PW/m2 (164 GW/cm2). Through focusing, such value is easily achievable with typical ps-lasers.

8. Results and Discussion

With the aim to characterize the significance, or not, of pump depletion in nonlinear multilayer microcavities, we investigate the effects of three important parameters; (1) intensity of the incident (external) pump wave (Figure 5), (2) thickness of the nonlinear layer (Figure 6), and (3) nonlinear susceptibility of the nonlinear layer (Figure 7). For all figures, a continuous (blue) curve means “without pump depletion” while a dashed (red) curve means “with pump depletion.”

In Figure 5, the calculations are performed in the cases of PMND and PMWD, for a “thin” microcavity, that is a “6×𝜆/2 microcavity” (𝜆=1,19μm). This thickness allows one mode in the center of the mirror's bandgap and two pseudomodes. In Figure 5 the intensity of the incident pump wave varies from 5,31 GW/m2 (as previously) to 531 GW/m2. The effect of pump depletion is clearly seen even at weak intensity and is dominant for the central mode, relatively to the pseudomodes. This is expected because the central mode suffers much more roundtrips than pseudomodes. Also anticipated and in agreement with the results from the undepleted case, an increase in incoming pump wave intensity by one order of magnitude increases the outgoing harmonic wave by two orders of magnitude. Also expected and shown in Figures 5(b) and 5(c), neglecting the pump depletion might lead to such a large computational error that the outgoing harmonic intensity can be calculated having a higher value than the incoming pump wave intensity, which is not the case when taking into account the pump depletion. Finally, it worth noting that even the “strong” incident pump wave intensity (531 GW/m2) is far from the strongest possible intensity.

In Figure 6, the calculations are performed in the cases of NMND and NMWD, at the “intermediate” incident pump wave intensity (53,1 GW/m2). In Figure 6 the thickness of the nonlinear layer varies from 10 μm to 11 μm. The 10 μm thickness corresponds very closely to an odd number, 11, of nonlinear coherence lengths (in this work, we adopt the definition of the nonlinear coherence length, 𝑙coh, as the distance between a minimum and a maximum in the spatial profile of SHG intensity. Formally, here, 𝑙coh=𝜋/|Δ𝑘|. Some authors define the nonlinear coherence length as twice this value), so that SHG is maximized. A contrario, The 11 μm thickness corresponds very closely to an even number, 12, of coherence lengths, so that SHG is minimized. As expected, the layer having an odd number of coherence lengths in thickness yield a greater SHG wave than for an even number of coherence lengths.

In Figure 7, the calculations are performed in the cases of PMND and PMWD, for a “thin” microcavity, that is a “6×𝜆/2 microcavity” (𝜆=1,19μm), at the “intermediate” incident pump wave intensity (53,1 GW/m2). In Figure 7 the nonlinear susceptibility varies from 𝜒(2)=100 pm/V to 𝜒(2)=10 pm/V. As expected from the undepleted case, a decrease in nonlinearity by one order of magnitude decreases the resulting outgoing harmonic wave by two orders of magnitude. From the numerical values found in Figure 7, neglecting the pump depletion for the current structure (a typical multilayered microcavity embedding a nonlinear medium the order of 100 pm/V) when submitted to a moderate incident pump wave intensity (~50 GW/m2), causes an inaccuracy as large as 350% in the intensity of the outgoing harmonic wave. This error is reduced to a mere 2% in presence of a moderately nonlinear medium (~10 pm/V). We expect a negligible error in case of a weakly nonlinear medium (~1 pm/V). Of course, this is not a general rule and various structures embedding arbitrary nonlinear medium and submitted to variable incident pump wave intensity should be covered with proper calculation.

9. Summary, Conclusions, and Perspectives

In summary, in presence of optically confining structures embedding nonlinear medium, even submitted to a low incident pump wave intensity, neglecting the pump depletion could lead to large errors in the computation of the intensity of the outgoing harmonic waves. This could have important consequences when light control and management in stratified nanostructures is involved. The matrix formalism described in this paper appears very suitable to this task and could be helpful in numerous interesting situations where the simulation of the optical properties on nonlinear-layered nanostructures is required. The complete formalism presented here is easily implemented on a computer in few hours.

The approach given in this paper is a step forward; a very general formalism devoted to the theory of nonlinear processes in layered structure. Especially when multiple layer (cumulative) pump depletion might occur. It would also be of interest to study different structures incorporating some nonlinear layers. The effect of various parameters such as their thickness and/or separation would strongly affect the emerging nonlinear waves. A forthcoming next step could be the extension of this approach to include subsidiary parameters neglected above such as, anisotropy, inhomogeneity, and strongly depleting secondary wave mixing.

Acknowledgment

One of the authors (SG) is very grateful to Professor Joseph Zyss for his assistance and support at the early stage of this work. The authors thank the Conseil de recherches en sciences naturelles et en génie du Canada, the Fondation de l’innovation du Nouveau-Brunswick and the Faculté des études supérieures et de la recherche de l’Université de Moncton for their financial support.