Analytical model of surface second-harmonic generation

The process of second-harmonic generation (SHG) in a finite one-dimensional nonlinear medium is analyzed in parallel by the Green-function technique and the Fourier-transform method. Considering the fundamental pump field propagating along a given direction and eliminating back-reflections at the boundaries the terms giving the surface second-harmonic fields in the particular solution of the wave equation are uniquely identified. Using these terms the flow of energy corresponding to the surface second-harmonic fields is analyzed in the vicinity of the boundaries. The formula giving the depth of the nonlinear medium contributing to the surface SHG is obtained. Both approaches for describing the SHG are compared considering complexity and quantization of the interacting fields. In addition, a theoretical model of surface SHG in centrosymmetric media is proposed. The model is built upon assumption that the second-order nonlinearity decays exponentially with distance from the boundary. As an important example, the generation of surface SHG from a thin dielectric nonlinear layer placed on a silicon substrate is analyzed by the proposed model.

Despite the fact that this phenomenological approach to surface SHG does not directly incorporate local changes of the electronic structure of atoms at the interface and thus neglects certain part of the surface SHG it has been found very useful. The reason is that the underlying formulas are written only for the optical-field amplitudes and so they can be relatively easily applied to more complex nonlinear optical systems. This is especially important under two different situations. First, when complex nonlinear photonic structures are analyzed. In this case, this simplified approach allows to keep certain control above the structure of the nonlinear contributions without sheer relying on the results of complex numerical calculations that involve the dynamics of induced electric dipoles. The second case occurs when quantum properties of the interacting fields are important. Here, the consideration of quantum fields poses the question how to consistently quantize an optical field in the simplest way compatible with the description of the nonlinear process.
In the framework of this phenomenological approach, two different methods have been developed: The first one is based upon the rigorous (and apparently mathematically consistent) approach that uses the Green-function technique. In the second method, the Fourier transform is applied to the interacting fields in the infinite volume to find a general ('generic') solution and the boundary conditions are incorporated subsequently using the field amplitudes' continuity relations at the boundaries.
The second method has been used much more frequently in the past as it is easier for application due to its clearly defined steps that allow for safe implementation even in more complex nonlinear photonic structures (e.g., using the transfer-matrix approach). Moreover, the interacting optical fields in this method may be conveniently quantized in an infinite volume after their decomposition into orthonormal monochromatic plane waves 17 . In fact, this formulation is close to the original description by Blombergen and Pershan who wrote a particular solution of the wave equation for the second-harmonic field in the form of harmonic plane waves and, alternatively, applied the Fourier transform method without considering the boundaries 7,14,15,18 . They arrived at the surface SHG just after considering the continuity requirements for the electric-and magnetic-field amplitudes at the boundaries. The behavior of optical fields at the boundaries is consistent with the optical-field quantization procedure as long as we invoke the quantization of photon-field flux (not energy, as it is usually done in quantum optics) 19 . For this type of quantization, only 'continuous' propagation of photons through both the involved linear and nonlinear media is required.
On the other hand, the first phenomenological method gives us the mathematically rigorous formula for second-harmonic field in terms of the Green functions of the nonlinear wave equation considered together with the boundary conditions. The Green function gives the solution with an impulse source function placed at a varying position inside the nonlinear medium 20,21 . So it allows for a physically-sound interpretation based on identifying individual material nonlinear dipoles at the positions of the impulse source functions, where they serve as physical sources of the emitted second-harmonic field. However, the structure of such solution is much less transparent and also the relation to the quantization scheme based on monochromatic plane waves is not straightforward.
Moreover, the analysis of surface SHG performed up to now assumed different indices of refraction for the media at different sides of a boundary. Thus, the nonlinear effects around the boundary were influenced by linear properties of the interacting fields. In the limiting, but common, case of surface SHG generated in photonic structures with more than one boundary, the effects of surface SHG emission and zig-zag propagation of the SHG field are mixed such that they cannot be separated. Both effects can only be safely separated when the effects of linear propagation are eliminated by considering the special case in which both linear and nonlinear media inside a photonic structure share the same index of refraction. Only in this very specific case, that is analyzed in the paper, the genuine properties of the SHG field emitted around the boundary are revealed. In the paper, we analyze this case from the point of view of both methods and show how the surface SHG is incorporated in both of them. Moreover, by comparing side by side all terms contained in two methods, we demonstrate in detail their close structural similarity that results in identical formulas for the second-harmonic field. In this comparison, we clearly identify both photons emitted in the volume and photons emerging at the surfaces of different media. To address all important physical aspects of the SHG process, it is sufficient to analyze the simplest case of collinear (1D) propagation of both fields. This brings us, among others, detailed understanding of the flow of the Poynting vector of the surface second-harmonic field in the vicinity of the interface. This extends the conclusions drawn both in the original papers 14,15 and recent literature 7,16 . We show that surface SHG propagates in the vicinity of the boundary always in the direction against the pump fundamental field, even inside the nonlinear medium. The overall second-harmonic field, consisting both from the volume and surface contributions, then propagates also against the pump-field direction in the close vicinity of the boundary and begins to co-propagate with the pump field typically several tens of nm beyond the boundary.
In addition, we address in the paper the problem of description of surface SHG in centrosymmetric nonlinear crystals with zero volume χ (2) susceptibility. We show that the considered phenomenological model allows to treat such SHG fields by defining a suitable spatial profile of χ (2) susceptibility around the boundary. If this profile decays exponentially from the boundary into the volume of the centrosymmetric medium, simple analytical formulas are obtained. Applying this model, surface SHG from a thin layer placed on a silicon substrate is analyzed as an important practical example.
The paper is structured as follows. We formulate and solve the model of SHG in section Second-harmonic generation treated by the Green-function technique first. Then, the model is analyzed in section Second-harmonic generation treated by the Fourier-transform method. The subsequent section brings the Poynting-vector analysis of second-harmonic generation near the boundaries. The model of surface SHG in centrosymmetric media and its application to a nonlinear layer on a substrate is presented in section Surface second-harmonic generation in centrosymmetric crystals. In the last section, Conclusions are drawn.

second-Harmonic Generation treated By the Green-Function technique
At the phenomenological level, the process of SHG is described by the second-order term in the Taylor expansion of the nonlinear polarization vector P nl (r, t) considered as a function of the fundamental (pump) electric-field amplitude vector E(r, t) 7,17,22 In Eq. (1), nonlinear properties of the medium responsible for SHG are described by the second-order nonlinear tensor χ (2) and ε 0 means the vacuum permittivity. According to Eq. (1), if the medium is irradiated by a strong fundamental field at frequency ω, the induced nonlinear polarization P nl generates a second-harmonic field E at doubled frequency 2ω.
To present the analysis in its simplest form, that, however, keeps all the substantial features of the nonlinear process, we concentrate on 1D propagation of the fundamental and second-harmonic fields (along the z axis). We put the nonlinear medium in between the boundaries positioned at z 1 and z 2 (see Fig. 1). We assume all the media being homogeneous and isotropic. Also we consider the fundamental field to be undepleted by the nonlinear process.
Under these conditions, the cartesian components E j , j ∈ {x, y, z}, of the second-harmonic field at frequency 2ω obey the following form of the wave equation: using the fundamental electric-field amplitude E(r, ω). The solution of the wave equation (2) can be written as a sum of a homogeneous solution E H (z, 2ω) and a particular solution E P (z, 2ω). In our model, the electric-field amplitude E H (z, 2ω) of the homogeneous solution has to be perpendicular to the z axis and so it can be decomposed into the plane waves polarized along the x and y axes propagating in both forward (index F) and backward (B) directions. So, we can express the homogeneous solution E H in the form The complex amplitudes A j a belong to plane waves with wave vectors k a propagating in either forward (a = F, k F ≡ k) or backward direction (a = B, k B ≡ −k) and polarized along either the x (e x ) or y (e y ) axis. We note that k(ω) = n(ω)ω/c. The complex amplitudes A j a are commonly used to describe the fields that enter into the analyzed nonlinear medium through outside at the boundaries at z 1 and z 2 or are scattered by the boundaries 21 . Figure 1. Scheme of the investigated nonlinear structure positioned between z = z 1 and z = z 2 . In our 1D model, the structure is infinite in the transverse directions x and y. All media have the same linear indices of refraction n. The pump beam (with angular frequency ω) responsible for SHG impinges on the structure from the left and it propagates only along the +z axis (as indicated by the red arrows). Contrary to this, the second-harmonic field (with angular frequency 2ω) generated by the nonlinear process propagates both along the −z and +z axis depending on the position (see the purple arrows). The fact that the second-harmonic field propagates along the −z axis not only in the area z < z 1 but also for certain z > z 1 shows that the surface SHG is not a 'point process' , but it occurs in certain surroundings of the boundary.
www.nature.com/scientificreports www.nature.com/scientificreports/ On the other hand, the particular solution E P describes the fields emitted by the sources inside the nonlinear medium, in the process of SHG. The determination of the second-harmonic field thus means to find an appropriate particular solution E P (in parallel with the homogeneous one E H ). In this section, we reveal the looked-for particular solution E P by the Green-function technique. This technique is based on the determination of the Green function G j (z, z′, 2ω) that satisfies the following equation: symbol δ stands for the Dirac function. We note that the Green function G is in general a second-order tensor for vectorial tasks. However, in our configuration the only nonzero terms are the diagonal ones. Applying the Green function G j (z, z′, 2ω), we may express a j-th electric-field component E j P of the particular solution along the formula In Eq. (6) the Green-function components G j determine the second-harmonic electric field vector E P (z, 2ω) at coordinate z emitted by nonlinear polarization vector P nl (z′, 2ω) at position z′. Physically, the Green function relates the second-harmonic field at position z with its source represented by a single oscillating nonlinear dipole at position z′. Since there is a continuum of dipoles in between the coordinates z 1 and z 2 , the overall electric field E P (z, 2ω) is given by an integral that adds contributions from all dipoles oscillating in the nonlinear medium. Explicitly, the z component G z (z, z′, 2ω) of the Green function is equal to It follows from Eqs (6) and (7) that the z component (longitudinal) E z P of the second-harmonic electric-field amplitude is nonzero only in the nonlinear medium z ∈ (z 1 , z 2 ) where it attains the simple form: On the other hand, the x and y components G x (z, z′, 2ω) and G y (z, z′, 2ω) of the Green function are derived in the form of fields propagating to position z along the +z or −z direction from a dipole at position z′: θ(z) stands for the Heaviside step function, i.e. θ(z) = 1 for z > 0, θ(z = 0) = 1/2 and is zero otherwise.
To obtain the general solution described by the Green functions (7) and (9) explicitly, we assume the fundamental pump field to be composed of monochromatic plane waves with wave vectors k b (ω) and complex ampli- Then, the particular solution of second-harmonic electric-field amplitude E P (z, 2ω) provided by the Green-function technique is written as:   1 2 In Eqs (11) and (12), (2) and symbol: shorthands the susceptibility tensor χ (2) with respect to its last two indices. Function rect [a,b] (z) equals 1 for z ∈ 〈a, b〉 and is zero otherwise.
According to Eq. (12), the z component ω E z ( , 2 ) z P of the second-harmonic electric-field amplitude is nonzero only inside the nonlinear medium where it does not depend explicitly on the positions z 1 and z 2 of the boundaries. Similarly, this property is inherited to the first terms on the third and fourth lines of Eq. (11) that give the x and y y P , respectively. Local values of the electric-field amplitudes described by these terms thus do not depend on the positions of the boundaries and so we call these terms as generic, i.e. being general to any geometry of the nonlinear medium. Moreover, these terms contain only the doubled wave vectors k(ω) of the fundamental field.
Contrary to this, the remaining terms in Eq. (11) involve spatial harmonic variations with wave vectors ±k(2ω) of the second-harmonic field and explicit dependence on the positions z 1 and z 2 of the boundaries. These terms thus explicitly depend on the geometry of the nonlinear medium. In detail, the terms on the second [fifth] line of Eq. (11)  This division of all terms into the generic and specific ones reveals the structure of the particular solution for the second-harmonic field provided by the Green-function technique. The knowledge of this structure qualitatively helps us in the analysis of more complex nonlinear structures as well as structures emitting weak, i.e. quantum, fields. It allows to reach the appropriate particular solution in two subsequent steps. In the first step, the generic terms are obtained depending on the type of nonlinear interaction and form of the pump field. In the second step, specific terms are added to the generic terms such that the second ones accord with the boundary conditions. This means that the continuity of optical fields at the boundaries has to be guaranteed. This is typically accomplished by applying the transfer-matrix method. The discussed division is also important for quantum fields. It allows to justify the quantization scheme developed for infinite volume 23 . In this scheme, the quantization is based on the generic terms. If we invoke the photon flux quantization 19 , we immediately get the quantization that is consistent with the propagation of light at the boundary.

second-Harmonic Generation treated by the Fourier-transform Method
We show in this section, that the Fourier-transform method that first solves the wave equation for a nonlinear medium in the whole infinite space and then incorporates the boundary conditions exactly fits into the structure of the particular solution of the second-harmonic field as described in the last paragraph of the previous section.
In the Fourier-transform method the second-harmonic field E P (z, 2ω) is decomposed into harmonic plane waves whose complex amplitudes determine the spatial spectrum The solution of the wave equation (2) in its Fourier domain is then simply written as: Considering the fundamental field E(z, ω) in the form (10) and formula (3) for the nonlinear polarization P nl , the inverse Fourier transform applied to the spectral solution in Eq. (14) gives us the z component E z (z, 2ω) of the second-harmonic field in the form (12) reached by the Green-function method. Contrary to this, the x and y components of the second-harmonic field are obtained in the following form: The terms written in Eq. (15) coincide with those in Eq. (11) that were identified as the generic terms. Thus, the Fourier-transform method provided us all generic terms of the Green-function approach, which is needed to accomplish the first step in the determination of the particular solution as outlined in the previous section. We note that, from the mathematical point of view, the solution in Eq. (15) also represents a particular solution to the wave equation (2), which is however different from that found by the Green-function technique. We note that this solution was found for the first time by Blombergen and Pershan 14 .
We also remind the striking feature of the solution in Eq. (15): It does not contain the terms evolving with the second-harmonic wave vectors ±k(2ω). It contains only the terms with doubled spatial frequency k(ω) of the pump field. Thus the free-field propagation of the second-harmonic field from the emitting nonlinear dipoles is not explicitly visible in the solution. The reason is that, for an infinitely extended nonlinear medium, the (2019) 9:4679 | https://doi.org/10.1038/s41598-019-39260-9 www.nature.com/scientificreports www.nature.com/scientificreports/ second-harmonic waves propagating freely with wave vectors −k(2ω) from the right-hand side of an arbitrary point z completely destructively interfere with those coming from the left-hand side with wave vectors +k(2ω). As a consequence, the formula in Eq. (15) gives only the electric-field amplitude generated by the local nonlinear polarization P nl at position z.
The solution in Eq. (15) diverges in the phase-matching limit 14,15 . This divergence is then naturally removed when a finite-length nonlinear medium is analyzed, i.e. when the positions z 1 and z 2 of the boundaries are taken into account. Indeed, the relations when applied to the expressions in the third and fourth lines of the solution (11) reached by the Green-function technique guarantee regularity of the whole formula in Eq. (11). Now we perform the second step in which we include the boundary conditions at positions z 1 and z 2 . This means that we extend the particular solution (15) by a suitable homogeneous solution such that the boundary conditions for the electric and magnetic fields are fulfilled. The looked-for particular solution is such that there are no external second-harmonic fields impinging onto the nonlinear medium from outside. This means that the electric-field components E j of the looked-for particular solution have the following general form outside the nonlinear medium:  (18) below, parameterize the general solution in a jth area specified in the lower index (1 -in front of the nonlinear medium, 2 -inside the nonlinear medium, 3 -beyond the nonlinear medium). On the other hand, a general solution in the nonlinear medium is written in the form: The continuity requirements at the boundaries at z 1 and z 2 leaves us with the following set of four linear algebraic equations for amplitudes A j 1,  The obtained general solution described by the amplitudes in Eq. (20) coincides with the particular solution E j P obtained by the Green function method in Eq. (11). At the end, we would like to note that the amplitudes of the homogeneous solution E H [see Eq. (20)] that corrects the original particular solution E P depend linearly on the nonlinear coupling constant d. This means that, from the physical point of view, such homogeneous terms should be automatically considered as a part of the particular solution as they describe the second-harmonic field emerging in the nonlinear process.
When we return back to the particular solution for the second-harmonic field obtained by the Green-function technique and written in Eqs (11) and (12) we may formulate the following physical picture of SHG. The second-harmonic field E(z, 2ω) is given as a coherent superposition of local second-harmonic fields emitted by www.nature.com/scientificreports www.nature.com/scientificreports/ individual nonlinear dipoles. Since a radiating dipole emits its radiation equally in the forward and backward directions [see the Green function (9) symmetric with respect to the exchange z ↔ −z] the second-harmonic field is expected in general in both the forward and backward directions. Especially, in our geometry, the second-harmonic field propagating along the −z axis from the nonlinear medium (z < z 1 ) is expected. This field is built by the radiation emitted into this direction by the nonlinear dipoles situated near the boundary at z 1 . The reason is that such radiation cannot be compensated by the radiation propagating in the opposed direction. However, this argumentation can be applied only to a narrow layer surrounding the boundary. This means that the surface SHG is not a 'point' process found just in the plane of the boundary, but it occurs in a narrow layer around the boundary. We note that the original studies by Blombergen et al. 14,15 as well as the newer investigations by Mlejnek et al. 16 did not address this issue.

poynting-Vector Analysis of second-Harmonic Generation Near the Boundaries
In this section, we analyze the formula for the second-harmonic field in Eq. (11) in the vicinity of the boundaries. This allows us to reveal the behavior of the surface contribution to the second-harmonic field. In the investigation, we determine the flow of energy of the second-harmonic field as given by the z component S z of the Poynting vector. We first analyze all terms of the second-harmonic field given in Eq. (11) with respect to the field behavior around the boundary at z 1 .
The first term in line two and the second term in line four in Eq. (11) depend on z 2 and so they describe the second-harmonic field propagating from the boundary at position z 2 in the −z direction. The second terms in lines two and three contribute to the field around the boundary at position z 1 . Physically the most important term is the second one on line two that describes constant flow of energy from the nonlinear medium from the area around z 1 (see Fig. 1) back to the free space (−z direction), i.e. in the propagation direction opposed to that of the pump beam. This term quantifies the second-harmonic field arising from the region around the boundary at z 1 . Fixing again the indices jkl and bg, the corresponding constant energy flow component We note that the energy flow is a scalar quantity in our 1D model. As the energy flow component − S z ( ) z 1 is continuous at the boundary at z 1 , it has to flow in the −z direction also inside the nonlinear medium, in the close vicinity of z 1 . In this area, the energy flow + S z 1 is approximately described as: where the energy flows contra-intuitively in the −z direction. This energy belongs to the surface second-harmonic generation [14][15][16] .
To complete the analysis of surface SHG, let us continue to determine the partial Poynting vector − S z ( ) z 2 describing the field propagating back from the area around the boundary at position z 2 , as given by the first term in line two and second term in line four in Eq. (11). Its magnitude apparently equals that of the partial Poynting vector − S z ( ) z 1 given in Eq. (21) for the field emitted around the boundary at z 1 . Moreover, the surface fields generated around both boundaries at z 1 and z 2 may interfere depending on the length L ≡ z 2 − z 1 of the nonlinear medium. This interference then leads to the complete Poynting vector S at position z 1 encompassing both surface second-harmonic fields: ). In this case, the surface second-harmonic field generated around the boundary at + z 1 perfectly compensates that emitted around the boundary at − z 2 . Another specific case occurs when both surface contributions maximally constructively interfere ]. The dependence of three analyzed Poynting vectors, S z 1 , − S z 2 and S, on position z around the boundary at z 1 is plotted in Fig. 2. Assigning sign -to the energy flow in the −z direction, the Poynting vector + S z 1 only increases with z, whereas the Poynting vector − S z 2 is constant. The overall Poynting vector S attains its minimal value at the boundary at z 1 when the maximally constructive interference between both terms occurs. In this case, the Poynting vector S increases with the increasing z and it becomes positive at z = 32 nm, i.e. at > + z z z 1 and + z z 1 is defined in Eq. (23). When we change the crystal length L such that the constructive interference between two surface second-harmonic www.nature.com/scientificreports www.nature.com/scientificreports/ fields weakens, the minimum of the overall Poynting vector S shifts into the nonlinear medium to position ∈ 〈 〉 + z z z , z 1 1 , attains larger values and the area with the negative overall Poynting vector S broadens first (see the solid curve with ○ in Fig. 2). Then, this area starts to shrink and the minima of the Poynting vector S approach to the value S = 0. At certain point, the overall Poynting vector S attains only non-negative values. This is the case of the completely destructive interference of two surface second-harmonic fields. If we continue with the change of crystal length L further, the negative minimum of the Poynting vector S found at the position z 1 sinks down and the extension of the area of z with the negative Poynting vector S gradually broadens from zero until the case of the completely constructive interference is reached. Gradual evolution of length l quantifying the extension of the area with the negative Poynting vector S observed for the increasing length L of the nonlinear medium in depicted in Fig. 3. According to the curves in Fig. 3, the greatest value of the length l equals 42 nm. We note that the discussed interference of surface second-harmonic fields coming from two boundaries can be avoided by either considering non-collinear geometry 15 or pulsed fields 16 .

surface second-Harmonic Generation in Centrosymmetric Crystals
In centrosymmetric crystals the bulk χ (2) susceptibility is zero. As a consequence, the above formulas predict no surface SHG. However, the above presented model can be simply modified to account for the presence of surface SHG in these crystals. At microscopic level, electronic structure of atoms just in the close vicinity of a boundary is modified such that it forms nonlinear dipoles that emit the SHG field 3,4,6 . The more distant the atom from a boundary is, the smaller the modification of its structure. We can phenomenologically describe such behavior by introducing a suitable profile of χ (2) susceptibility, e.g. in its simplest exponential form: (2)  www.nature.com/scientificreports www.nature.com/scientificreports/ where l nl gives the extension of the area with the modified electronic structure of atoms.
Considering this profile in Eqs (3) and (6), we arrive at the modified formulas (11) and (12) for the electric-field amplitudes: and l nl,j quantifies the extension of χ (2)   www.nature.com/scientificreports www.nature.com/scientificreports/ The particular solution ω E z ( , 2 ) j,1 P is given by Eqs (26) and (27) in which we substitute l nl by l nl,1 . The particular solution ω E z ( , 2 ) j,2 P arising from the χ (2) (z) profile around position z 2 is derived in a similar form as that given in Eqs (26) and (27): nl,2 and the susceptibility χ (2) (z 2 ) at position z 2 determines the nonlinear coefficients d jkl bg ,2 . Up to now, we have strictly assumed identical indices of refraction of all media to eliminate the effects of fields' back scattering. However, as back-scattering on boundaries separating media with different indices of refraction is an inevitable part of fields' propagation in real structures, we show here how to generalize the above formulas for a centrosymmetric crystal to account for different indices of refraction in front of and beyond the crystal. We assume that the medium in front of (beyond) the nonlinear crystal [denoted as 1 (3)] has index of refraction n 1 (n 3 ) whereas the nonlinear crystal [denoted as 2] is characterized by index of refraction n 2 . First, we apply the transfer-matrix formalism 21,24 to the fundamental field to reveal its forward-and backward-propagating components in all three media. This allows us to determine the particular solution E j P for the second-harmonic field of the form of Eqs (26) and (27), i.e. without including back-scattering effects. This corresponds to the application of the Fourier-transform method with spatially dependent χ (2) nonlinearity. We note that both the forward-and backward-propagating fundamental fields contribute to the solution E j P . This solution then gives the source terms to the transfer-matrix formulation of the second-harmonic field propagation in the nonlinear crystal. If we assume the x polarized second-harmonic electric-field amplitude E x P , the source terms S 1 and S 2 are given as: In Eq. (32), the transition matrices  l for medium l, l = 1, 2, 3, and the propagation matrix P 2 for the nonlinear crystal are defined as and the wave vector k 2 (2ω) belongs to the second-harmonic field inside the crystal.
In the second step, the application of the transfer-matrix formalism provides us the system of two linear algebraic equations for the amplitudes of the reflected ω A (2 )  using the matrix elements  jl . Finally, the electric-field amplitude E x (z, 2ω) of the second-harmonic field at an arbitrary position z is easily derived from the amplitudes ω A (2 ) and ω A (2 ) x 3, F using the Fresnel relations at the boundaries and free-field propagation.
As an example, using the developed model we numerically analyze the SHG from a centrosymmetric nonlinear layer of thickness L = 64.5 nm fabricated on a silicon substrate and surrounded by air (parameters are given in the caption to Fig. 5). The generation of surface second-harmonic field around the boundary at z 1 = 0 is much more intense compared to that around the boundary at z 2 = L, as documented by the curves giving the appropriate Poynting vectors S z 1 and S z 2 in Fig. 5. This unbalance in surface SHG originating around the boundaries at z 1 and z 2 is caused by back-scattering of the second-harmonic field that considerably suppresses the electric-field amplitude E x in the vicinity of the boundary at z 2 (see the dashed curve in Fig. 5). The complete Poynting vector S derived from the second-harmonic field then reflects interference of both terms arising around the boundaries at z 1 and z 2 . Comparing the curves for different Poynting vectors drawn in Fig. 5 we reveal that this interference is constructive (destructive) in the area around the boundary at z 1 (z 2 ). Moreover, the second-harmonic field propagating in air is more than twice intense compared to that propagating in the substrate.