## Abstract

We study theoretically and numerically the second harmonic generation in a nonlinear crystal with random distribution of ferroelectric domains. We show that the specific features of disordered domain structure greatly affect the emission pattern of the generated harmonics. This phenomena can be used to characterize the degree of disorder in nonlinear photonic structures.

©2010 Optical Society of America

## Introduction

Quasi-phase matching is one possible method of choice to compensate for the lack of exact phase synchronism among interacting waves in quadratic nonlinear media ensuring a high efficiency of frequency mixing [1,2]. It employs spatial periodic changes in the sign of nonlinearity. This is achieved by periodic poling of ferroelectric nonlinear crystals such as lithium niobate, lithium tantalate, etc [3]. However, the exact periodic poling ensures the fulfilment of the phase matching condition only for a particular choice of wavelengths of the interacting waves; thus restricting its practical applicability. To remove this constrain, quasi-periodic poled structures or sequences of periodic regions with different periods have been realized [4, 5]. Naturally grown ferroelectric crystals such as Strontium Barium Niobate (SBN) often exhibit an irregular multi-domain structure with rod-like domains [6–9]. Because of the random size and spatial distribution of these domains, the unpoled SBN crystal belongs to the class of disordered nonlinear quadratic crystals. Such media have attracted a considerable attention in recent years due to their potential benefits for nonlinear parametric processes [10]. Despite being usually considered as disadvantageous in typical nonlinear optical phenomena, the randomness of the nonlinear media offers potential applications in broadband frequency generation [11–18] and ultra-short pulse monitoring [19]. These applications utilize the fact that the crystal with disordered anti-parallel domains forms, in fact, a two-dimensional photonic nonlinear structure with an infinite number of *χ*
^{(2)} reciprocal vectors corresponding to various local nonlinear gratings. Consequently, such structure enables phase matching of various parametric processes over an ultra-broad spectral range and in different spatial directions, including a direction transverse to the beam propagation [15, 20, 21]. Figure 1(a) represents a schematic illustration of such type of random crystal placed in a typical experimental set-up. The inverted ferroelectric domains are represented by bright and dark colors, respectively.

However, recent experiments with random domain SBN crystals also demonstrated that the emission pattern of the generated waves may differ drastically from sample to sample [see Fig. 1(b), 1(c)], being either spatially homogeneous or exhibiting distinct intensity peaks [16, 22]. This behavior points towards the difference in the actual distributions of ferroelectric domains in different samples. In fact, observations of random domains ranging from tens of nanometers to few microns have been reported in the literature [6, 8, 9, 13]. Moreover, it has been even suggested that the domain-walls induced light scattering of the fundamental wave plays significant role in the SH generation [17].

In order to obtain a better understanding of the nonlinear interaction in random media, in this work we analyze theoretically and numerically the role of the domain distribution on the efficiency and transverse characteristics of the SH emission in a quadratic random crystal. The numerical simulations obtained are supported by our experimental evidence of the SH emission (see for example Fig. 1(b), 1(c) and Ref. [20]).

## Modeling random ferroelectric domain pattern

In order to simulate the 2D ferroelectric domain distribution for our study we will consider a simplified model which includes the most significant aspects of the experimentally reported distributions in SBN crystals [9, 17]. To mimic the process of domain formation we adopt the following strategy. Firstly, we assume that the individual domains (our building blocks) have form of rods with a circular transverse profile. While the assumption about circular shape follows from the experimental observations [23] it is, in fact, not that critical because, as we see below, the macroscopic pattern of nonlinearity distribution consists of regions of opposite sign of the nonlinearity which may have very complicated (not necessarily circular) shape. As we consider light propagation and emission in the plane perpendicular to the longer dimension of the domains only their transverse structure is taken into account here. A normal (Gaussian) distribution of the domain sizes is chosen with a certain mean diameter (*ρ*
_{0}) and variance (σ). Using this distribution a number of domains is randomly generated and randomly placed in a rectangular area representing the size of the sample. In all the cases discussed here at least two thousand domains were used to create a random domain pattern even if the domains number that take part in the process may depend on actual geometry of interaction. As a rule of thumb this number should be large compared to the actual size of the input optical beam [24]. In the next step each domain is randomly assigned a particular sign of the nonlinearity (positive/negative). Some particular realizations of the 2D domains pattern for different choice of statistical parameters are shown in the top row of Fig. 2(a)–2(c). The black and white colors denote the antiparallel oriented domains with opposite signs of the quadratic nonlinearity. These kind of generated patterns are directly used in the numerical simulations. A completely random assignment leads sometimes to large agglomeration of domains with the same sign of the *χ*
^{(2)} nonlinearity [Fig.2(c) top row]. To remedy this, and to create structures with very fine domain pattern the sign of each particular domain was assigned taking into account the signs of its closest neighborhood, such as in Figs. 2(a) and 2(b). The algorithm implemented here integrates the area close to the domain under study; if in this area there is a prevalence of positive *χ*
^{(2)} sign, a negative sign will be chosen and vice versa. The non occupied zones among the domains are considered to form the domain walls and they amount to around 10% of the whole area of the sample which is in agreement with the experimental observations in similar poled structures [25]. In the region of the domain walls the nonlinearity is assigned a null value. However, we allowed in this region for a small linear refractive index change of the order of 10^{-3} to account for the strong discontinuity of the vector of spontaneous polarization [25, 26]. This weak linear index variation will produce a linear scattering of light. We will talk about this effect later in the paper.

The insets in Fig. 2(a) and 2(b) depict corresponding histograms of normalized number of domains versus domain size. Such histograms were previously used to characterize actual samples of as grown SBN crystals with disordered domain pattern [17] but do not contain all the needed information on the distribution. For example, it is not possible to retrieve the distribution of the sign of the nonlinearity. Moreover, the histograms can be ambiguous since it is possible to choose either the diameter or the area of the domains as discriminating parameter. Notice how, starting from the same statistical parameters, it is possible to draw two completely different random patterns by either avoiding [Fig. 2(b)] or allowing the agglomeration [Fig. 2(c)].

## Second harmonic emission by disordered domain structure

In order to fully characterize the domain pattern as a quasi-phase matching nonlinear photonic structure we use the spatial Fourier spectrum of the nonlinearity distribution. Such spectrum represents the domain of reciprocal vectors *G⃗* which will be employed to fulfil the phase matching condition *k⃗*
_{2} − 2*k⃗*
_{1} − *G⃗* = 0 which determines the direction and efficiency of the quadratic process. Here *k⃗*
_{2}, *k⃗*
_{1} are the wave-vectors of the second harmonic and fundamental waves, respectively, and *G⃗* can be expressed as *G⃗* = (*π*/Λ)*n⃗* where Λ represents periodicity along the direction *n⃗*. In the bottom row in Fig. 2 we depict the modulus of the Fourier spectrum corresponding to the domain patterns shown in the top row. In these graphs the strength of the spectral components is directly represented by their brightness. It is clear that different domain distributions lead to drastically different structure of the reciprocal space which, as we will show later, will subsequently affect the second harmonic emission. In fact, this property has been utilized to control the shape of the generated SH beam in specially designed structures [27, 28]. As it is clear from the bottom row of the Fig. 2(a), 2(b), a real spatially random domain pattern leads to a circularly shaped Fourier spectrum. The radius of the circle will be directly related to the mean domain size and its width from the radius will be related with the standard deviation. In particular, domain distributions with narrow (broad) standard deviation lead to the correspondingly narrow (broad) width of the reciprocal distribution of the *G⃗* vectors. And domain distributions with small (big) mean domain size lead to the correspondingly large (small) radius of the circular distribution. As comparison, the bottom row of Fig. 2(c) depicts the spatial Fourier spectrum of the distribution with large rate of domains agglomeration. Here the structure consists effectively of large area domains of the size of the agglomerated areas, and consequently its spatial spectrum shrinks toward very small components. In order to see how the particular nonlinear structures affect emission of the second harmonic we need to consider the phase matching condition. In Fig. 3 we plot the three particular scenarios associated with the random domain structures. The green circle represents possible directions of the wave-vector of the second harmonic while the yellow ring represents the more considerable region of the reciprocal vectors *G⃗* provided by the randomness of the domain pattern. As the pattern is isotropic so is the distribution of the *G⃗* vectors. The most efficient second harmonic emission is expected to occur in the directions determined by the intersection of the green circle and yellow rings where the phase matching conditions is fulfilled with the stronger coefficients of the reciprocal spectrum. It is clear that because of availability of isotropically oriented *G* vectors the generation of the second harmonic is expected to be non-collinear with a broad spatial distribution of the intensity. The two cases shown in Fig. 3(a), 3(b) correspond to identical domain pattern but they differ in the wavelength of the fundamental wave. Hence, in the first case (a) the SH is expected to be emitted predominantly in two angular directions. On the other hand the same structure used with shorter fundamental wave [as in Fig. 3(b)] should lead to mostly forward emission centered around the direction of the fundamental wave. Finally, the domain pattern with much smaller average size is expected to lead to strongest emission oriented at large angles as in Fig. 3(c). As the geometric regions of the strongest phase matching shown in these three examples differ significantly one may expect this property to be reflected in the SH emission pattern. However, the actual spatial intensity distribution of the generated SH wave will be determined not only by the phase matching condition but also by the strength of the nonlinearity for the particular emission direction which is strongly affected by the randomness of the system.

In order to explore this further we employ the approach of Dolino [29] and Le Grand [24]. In this approach the intensity of the SH wave generated in the medium consisting of randomly distributed antiparallel domains can be expressed as

where *I*
_{(ω)} is the intensity of the fundamental wave and *d _{eff}* is the effective

*χ*

^{(2)}nonlinearity which depends on material parameters and geometry of the interaction. The function

*f*(

*q*,

*ρ*

_{0}, σ) [Eq. (2)] represents the effect of disorder in the domain distribution.

where *L* is the propagation distance and *q* = |Δ*k*| is the absolute value of the phase mismatch. The linear dependence of *I*
_{(2ω)} on propagation distance is the signature of disorder which makes the process of SH build-up to be incoherent [10, 30–32]. As the phase mismatch parameter *q* is in our geometry of interaction uniquely linked to the angle between propagation direction of the SH and fundamental waves, the relation (1) provides, in fact, angular distribution of the SH intensity. The contour plots in Fig. 4 illustrate the relation Eq. (1) as a function of emission angle and/or statistical parameters of the domain distribution. From Fig. 4(a) it is evident that for low domain dispersion the structure is almost perfectly periodic leading to very directed transverse emission determined by the phase matching condition which in this case is satisfied in two particular directions of the SH. As the domain dispersion increases the more G vectors contribute towards the phase matching which results in a forward emission of the SH. Plot in Fig. 4(b) illustrates the effect of the mean value of the domain diameter. For fine domain pattern (at given dispersion) many G vectors participate in the SH generation leading to very broad emission. As the average domain diameter increases the SH emission becomes less spread approaching the forward emission for large *ρ*
_{0}. Finally, in Fig. 4(c) we demonstrate the effect of varying the wavelength of the fundamental wave on the SH generation. In this particular example *ρ*
_{0} = σ = 1*μ*m. The dashed line in all graphs in Fig. 4 represents emission of the SH determined solely by the phase matching condition taking into account only one reciprocal vector |*G⃗mean*| = 2*π*/(2*ρ*
_{0}) corresponding to the mean diameter of the domains. It is clear that while this simplified approach is sufficiently accurate for weak disorder it leads to erroneous results for high degree of randomness of the periodic structure. In addition it cannot be used to determine the spatial distribution of the emitted SH.

## Numerical simulations

To test the above discussed analytical predictions we resorted to numerical simulations of second harmonic generation in random domain structures. We numerically solved the Maxwell’s equations assuming that the random quadratic medium (Fig. 2) is illuminated with a Gaussian fundamental pulse. We used the material parameters which correspond to those of as-grown SBN crystals [26]. It should be stressed here that Vidal and Martorell [33] were perhaps the first ones to investigate numerically second harmonic emission in multi-domain crystal. They considered only 1D structure and used the matrix transfer approach. Our simulations of full Maxwell equations were conducted using optical pulses and finite beams in two dimensional medium. Therefore we could describe SH emission at arbitrary angle and did not have to perform statistical averaging in order to obtain relevant physical quantities such as the energy of the generated waves. To simulate the light propagation through SBN crystal we use a numerical model similar to that discussed before in [34] with the addition of a transverse coordinate to include diffraction and propagation in a two-dimensional plane. We use a Lorentz oscillator model to describe material dispersion and, for simplicity, we assume a TM-polarized incident pump field. The electric and magnetic fields may be written as a superposition of the harmonics as follows:

where *H ^{lω}_{y}*(

*z,y,t*),

*E*(

^{lω}_{y}*z,y,t*) and

*E*(

^{lω}_{z}*z,y,t*) are generic, spatially- and temporally-dependent, complex envelope functions;

*k*and

*ω*are carrier wave vector and frequency, respectively, and

*l*is an integer. Equations (3), (4) are a convenient representation of the fields, and no

*a priori*assumptions are made about the field envelopes. The second order nonlinear polarization is

**P**

_{NL}=

*χ*

^{2}

**E**

^{2}. Expanding the field into its components yields nonlinear polarization terms at the fundamental and harmonics frequencies (taking only the leading terms):

*P*(

_{ω}*z,t*) = 2

*χ*

^{2}

_{ω}

*E*

^{*}

_{ω}

*E*

_{2ω}and

*P*2

*ω*(

*z,t*) = 2

*χ*

^{2}

_{2ω}

*E*

^{2}

_{2ω}. The linear response of the medium is described by a Lorentz oscillator model:

*ε*(

*ω*) = 1 −

*ω*

^{2}

_{p}/(

*ω*)

^{2}+

*iγω*−

*ω*

^{2}

_{r}) and

*μ*(

*ω*) = 1, where

*γ, ω*, and

_{p}*ω*are the damping coefficient, the plasma and resonance frequencies, respectively. Assuming that polarization and currents may be decomposed as in Eqs. (3), (4), we obtain the following Maxwell-Lorentz system of equations for the

_{r}*l*field components, in the scaled two-dimensional space (

^{th}*ỹ*, ξ) plus time (τ) coordinate system:

In Eqs. (5)–(10), the functions *J,P,P _{NL}* refer to linear electric currents, polarization, and nonlinear polarization, respectively. The coordinates are scaled so that

*ξ*=

*z*/

*λ*

_{0},

*ỹ*=

*y*/

*λ*

_{0},

*τ*=

*ct*/

*λ*

_{0},

*ω*

_{0}= 2

*πc*/

*λ*

_{0}, where

*λ*

_{0}= 1

*μ*m is just reference wavelength;

*γ, β*= 2

_{lω}*πlω*/

*ω*

_{0},

*β*= 2

_{rω}*πω*/

_{r}*ω*

_{0},

*ω*are the scaled damping coefficient, wave-vector, resonance and electric plasma frequencies for the

_{p}*l*harmonic, respectively.

^{th}*θ*is the angle of incidence of the pump field. The equations are solved using a split-step, fast Fourier transform-based pulse propagation algorithm that advances the fields in time [35]. Contrary to the most common nonlinear numerical codes, the one implemented here is able to solve the complete 2D vectorial Maxwell’s equations in the time domain without slowly varying envelope approximation. For this reason all angles of propagation and/or diffraction are allowed. Furthermore, no undepleted pump approximation is used. A number of numerical simulations have been performed with different working conditions. The fundamental pulse length in the plane of the nonlinear domain was in the range 6–30 optical cycles. We varied the fundamental wavelength from 790nm up to 1550nm. The data for the index of refraction were obtained fitting the Sellemier’s equation with the coefficient as in [26] to a general Lorentz curve. Figure 5 shows the results of some of these simulations obtained assuming the

_{i}*λ*= 800,1064nm. In this figure the plots in the left column depict the spatial intensity distribution of the fundamental and generated second harmonics for three different random samples (from the top to bottom). Middle column shows the transverse structure of the fundamental and second harmonics in the spatial Fourier space (or equivalently, in the far field) and the plots in the right column depicts the energy of the fundamental and second harmonics as a function of the time of propagation. In the Fig. 5(a) is depicted the typical propagation of the fundamental pulse. Initially located in the air, it impinges normally in the sample. Due to the index mismatch, part of the pulse is backward reflected and the remaining is forward propagating in the material. The k-spectrum of the fundamental field shows a very narrow peak with coordinates

_{F}*k*= 2

_{x}*π*/

*λ*and

*k*= 0 during the first part of the propagation in the air and

_{y}*k*= 2

_{x}*πn*(

*λ*)/

_{F}*λ*and

_{F}*k*= 0 after entering in the material, where

_{y}*n*(

*λ*) is the index of refraction at the fundamental frequency. Since the efficiency of the nonlinear process is low, the energy of the fundamental is practically constant during the propagation in the SBN. The graphs in Fig. 5 demonstrate dramatic variations in the far field of the emitted SH depending on the characteristics of the domain pattern. Note the differences in the emission angles in cases (b), (c) and (d) which correspond to the domain distributions from large to small average domain sizes, namely 3,0.9 and 0.3

_{F}*μ*m. The angular distribution of the intensity of the generated SH fields reflects the spatial distributions of the domains. As expected, the energy of the generated SH clearly displays a linear-like growth. This is indication of the incoherent character of the build-up of the SH and is in full agreement with earlier experiments and theoretical prediction [10, 30, 31]. The difference in the maximum reached values of the SH fields and energies is due to the different input intensities of the fundamental as initial condition. For example, in the simulation in Fig. 5(b) the input intensity is three orders of magnitude higher than the simulation in Fig. 5(c). Comparing the numerical simulations with the analytical predictions [Eq. (1)] we find them in very good agreement. Therefore the simple statistical approach appears to be a powerful tool in predicting spatial distribution of the generated harmonic in media with disordered domain structure. This suggests that experimental observation of second harmonic pattern can be employed to obtain information about the degree of the disorder of the ensuing domain structure. This aspect is currently being experimentally investigated. To give an idea of the power of this method, as it follows from the results of this study we can evaluate that to have second harmonic light propagating at the transverse and backward directions the sample needs to possess nonlinear domains at least of the order of 10 – 50nm. This consideration is very important as far as the choice of the domain visualization technique with appropriate resolution is concerned [36].

Finally, we would like to comment on the role of the light scattering on the SH emission. Due to the discontinuity in the direction of spontaneous polarization between domains a small jump of the linear index of refraction is present in the area of the domain walls. For instance, according to [26] this index modulation is Δ*n*
_{λ=852nm} = 1.0×10^{-4} and Δ*n*
_{λ=404nm} = 2.6×10^{-3}. Obviously the strength of light scattering induced by the domain walls will be at least one order of magnitude weaker for infrared than visible light. Thus, the SH light will always experience stronger scattering than the fundamental wave while propagating in the same sample. It has been suggested before that the wide angle emission of the second harmonic in experiments with multi-domain SBN crystals could be attributed to the emission of the SH signal by scattered fundamental beam [17]. To test this hypothesis we included the effect of index variation in the domain walls in our numerical simulations. In fact all graphs depicted in Fig. 5 have been obtained assuming the index jumps at fundamental frequency of Δ*n* = 0.01. No appreciable effect in both fundamental beam nor the second harmonic has been observed. This confirms that the transverse emission of the SH is caused solely by the spatial modulation of the nonlinearity resulting from the random domain pattern. However, it should be stressed that in the situation when the strong scattering of the second harmonic is present it may indeed affect its overall intensity distribution.

## Conclusions

In conclusion, we have studied the effect of random domain distribution on the emission of second harmonic generation in a nonlinear crystal with disordered ferroelectric domains. We demonstrated that there is direct connection between the statistical properties of the domain structure and the far field of the second harmonic. This effect can be employed to characterize the degree of randomness of the quasi-phase matching domain pattern. We are currently working on experimental verification of these results in SBN crystals with disordered domain structure.

## Acknowledgements

This work was supported by the Australian Research Council, Spanish Government (FIS2008-06024-C03-02), Catalan Government (2009 SGR 1168 and BE 2009) and the COST Action MP0702. V.R. thanks Cristian Ciracì for helpful discussions and suggestions with the images. The authors thank D. Neshev for collaboration.

## References and links

**1. **J. A. Armstrong, N. Bloembergen, J. Ducuing, and P. S. Pershan, “Interactions between Light Waves in a Nonlinear Dielectric,” Phys. Rev. **127**, 1918–1939 (1962). [CrossRef]

**2. **P. A. Franken and J. F. Ward, “Optical harmonics and nonlinear phenomena,” Rev. Mod. Phys. **35**23–39 (1963). [CrossRef]

**3. **M. M. Fejer, G. A. Magel, D. H. Jundt, and R. L. Byer, “Quasi-phase-matched 2nd harmonic-generation - tuning and tolerances,” IEEE J. Quantum Electron. **28**, 2631–2654 (1992). [CrossRef]

**4. **S. Zhu, Y. Y. Zhu, and N. B. Ming, “Quasi-phase-matched third-harmonic generation in a quasi-periodic optical superlattice,” Science **278**, 843–846 (1997). [CrossRef]

**5. **U. K. Sapaev and G. Assanto, “Efficient high-harmonic generation in engineered quasi-phase matching gratings,” Opt. Express **16**, 1–6 (2008). [CrossRef] [PubMed]

**6. **K. Terabe, S. Takekawa, M. Nakamura, K. Kitamura, S. Higuchi, Y. Gotoh, and A. Gruverman, “Imaging and engineering the nanoscale-domain structure of a Sr0.61Ba0.39Nb2O6 crystal using a scanning force microscope,” Appl. Phys. Lett. **81**, 2044–2046 (2002). [CrossRef]

**7. **M. O. Ramirez, D. Jaque, L. Ivleva, and L. E. Bausa, “Evaluation of ytterbium doped strontium barium niobate as a potential tunable laser crystal in the visible,” J. Appl. Phys. **95**, 6185–6191 (2004). [CrossRef]

**8. **V. V. Shvartsman, W. Kleemann, T. Lukasiewicz, and J. Dec, “Nanopolar structure in Sr_{x}Ba_{1-x}Nb_{2}O_{6} single crystals tuned by Sr/Ba ratio and investigated by piezoelectric force microscopy,” Phys. Rev. B **77**, 054105 (2008). [CrossRef]

**9. **R. V. Gainutdinov, T. R. Volk, O. A. Lysova, I. I. Razgonov, A. L. Tolstikhina, and L. I. Ivleva, “Recording of domains and regular domain patterns in strontium barium niobate crystals in the field of atomic force microscope,” Appl. Phys. B **95**, 505–512 (2009). [CrossRef]

**10. **M. Baudrier-Raybaut, R. Haidar, Ph. Kupecek, Ph. Lemasson, and E. Rosencher, “Random quasi-phase-matching in bulk polycrystalline isotropic nonlinear materials,” Nature (London) **432**, 374–376 (2004). [CrossRef] [PubMed]

**11. **M. Horowitz, A. Bekker, and B. Fischer, “Broad-band 2nd-harmonic generation in Sr_{x}Ba_{1-x}Nb_{2}O_{6} by spread-spectrum phase-matching with controllable domain gratings,” Appl. Phys. Lett. **62**, 2619–2621 (1993). [CrossRef]

**12. **S. Kawai, T. Ogawa, H. S. Lee, R. C. DeMattei, and R. S. Feigelson, “Second-harmonic generation from needlelike ferroelectric domains in Sr_{0.6}Ba_{0.4}Nb_{2}O_{6} single crystals,” Appl. Phys. Lett. **73**, 768–770 (1998). [CrossRef]

**13. **J. J. Romero, D. Jaque, J. García Solé, and A. A. Kaminskii, “Simultaneous generation of coherent light in the three fundamental colors by quasicylindrical ferroelectric domains in Sr_{0.6}Ba_{0.4}(NbO_{3})_{2},” Appl. Phys. Lett **81**, 4106–4108 (2002). [CrossRef]

**14. **J. J. Romero, C. Aragó, J. A. Gonzalo, D. Jaque, and J. García Solé, “Spectral and thermal properties of quasiphase- matching second-harmonic generation in Nd3+:Sr_{0.6}Ba_{0.4}(NbO_{3})_{2} multi-self-frequency-converter nonlinear crystals,” J. Appl. Phys. **93**, 3111–3113 (2003). [CrossRef]

**15. **A. R. Tunyagi, M. Ulex, and K. Betzler, “Noncollinear optical frequency doubling in strontium barium niobate,” Phys. Rev. Lett. **90**, 243901 (2003). [CrossRef] [PubMed]

**16. **R. Fischer, D. N. Neshev, S. M. Saltiel, W. Krolikowski, and Yu. S. Kivshar, “Broadband femtosecond frequency doubling in random media,” Appl. Phys. Lett. **89**, 191105 (2006). [CrossRef]

**17. **P. Molina, M. O. Ramirez, and L. E. Bausa, “Strontium Barium Niobate as a Multifunctional Two-Dimensional Nonlinear Photonic Glass,” Adv. Funct. Mater. **18**, 709–715 (2008). [CrossRef]

**18. **W. Wang, V. Roppo, K. Kalinowski, Y. Kong, D. N. Neshev, C. Cojocaru, J. Trull, R. Vilaseca, K. Staliunas, W. Krolikowski, S. M. Saltiel, and Yu. Kivshar “Third-harmonic generation via broadband cascading in disordered quadratic nonlinear media,” Opt. Express **17**, 20117–20123 (2009). [CrossRef] [PubMed]

**19. **R. Fischer, D. N. Neshev, S. M. Saltiel, A. A. Sukhorukov, W. Krolikowski, and Yu. S. Kivshar, “Monitoring ultrashort pulses by the transverse frequency doubling of counterpropagating pulses in random media,” Appl. Phys. Lett. **91**, 031104 (2007). [CrossRef]

**20. **V. Roppo, D. Dumay, J. Trull, C. Cojocaru, S. M. Saltiel, K. Staliunas, R. Vilaseca, D. N. Neshev, W. Krolikowski, and Yu. S. Kivshar, “Planar second-harmonic generation with noncollinear pumps in disordered media,” Opt. Express **16**, 14192–14199 (2008). [CrossRef] [PubMed]

**21. **S. M. Saltiel, D. N. Neshev, R. Fischer, W. Krolikowski, and Yu. S. Kivshar, “Spatiotemporal toroidal waves from the transverse second-harmonic generation,” Opt. Lett. **33**, 527–529 (2008). [CrossRef] [PubMed]

**22. **J. Trull, C. Cojocaru, R. Fischer, S. M. Saltiel, K. Staliunas, R. Herrero, R. Vilaseca, D. N. Neshev, W. Krolikowski, and Yu. S. Kivshar, “Second-harmonic parametric scattering in ferroelectric crystals with disordered nonlinear domain structures,” Opt. Express **15**, 15868–15877 (2007). [CrossRef] [PubMed]

**23. **J. J. Romero, D. Jaque, and J. Garca Sol, “Diffuse multiself-frequency conversion processes in the blue and green by quasicylindrical ferroelectric domains in Nd3 + :Sr0.6Ba0.4(NbO3)2 laser crystal,” Appl. Phys. Lett. **78**, 1961–1963 (2001). [CrossRef]

**24. **Y. Le Grand, D. Rouede, C. Odiu, R. Aubry, and S. Mattauch, “Second-harmonic scattering by domains in RbH_{2}PO_{4} ferroelectric,” Opt. Commun. **200**, 249–260 (2001). [CrossRef]

**25. **A. L. Aleksandrovskii, O. A. Gliko, I. I. Naumova, and V. I. Pryalkin, “Linear and nonlinear diffraction gratings in lithium niobate single crystals with a periodic domain structure,” Quantum Electron. **26**, 641–643 (1996). [CrossRef]

**26. **Th. Woike, T. Granzow, U. Dörfler, Ch. Poetsch, M. Wohlecke, and R. Pankrath, ”Refractive Indices of Congruently Melting Sr0:61Ba0:39Nb2O6,” Phys. Stat. Solidi A **186**, R13–R15 (2001). [CrossRef]

**27. **Y. Q. Qin, C. Zhang, and Y. Y. Zhu, et al. “Wave-front engineering by Huygens-Fresnel principle for nonlinear optical interactions in domain engineered structures,” Phys. Rev. Lett. **100**, 063902 (2008). [CrossRef] [PubMed]

**28. **T. Ellenbogen, N. Voloch-Bloch, A. Ganany-Padowicz, and A. Arie “Nonlinear generation and manipulation of Airy beams,” Nature Photonics **3**, 395–398 (2009). [CrossRef]

**29. **G. Dolino, “Effects of domain shapes on second harmonic scattering in Triglycine Sulfate,” Phys. Rev. B **6**, 4025–4035 (1972). [CrossRef]

**30. **E. Yu. Morozov, A. A. Kaminskii, A. S. Chirkin, and D. B. Yusupov, “Second Optical Harmonic Generation in Nonlinear Crystals with a Disordered Domain Structure,” JETP Lett. **73**, 647–650 (2001). [CrossRef]

**31. **E. Yu. Morozov and A. S. Chirkin, “Stochastic quasi-phase matching in nonlinear-optical crystals with an irregular domain structure,” Sov. J. Quantum Electron. **34**, 227–232 (2004). [CrossRef]

**32. **I. V. Shutov, I. A. Ozheredov, A. V. Shumitski, and A. S. Chirkin, “Second harmonic generation by femtosecond laser pulses in the Laue scheme,” Opt. Spectroscopy **105**, 79–84 (2008). [CrossRef]

**33. **X. Vidal and J. Martorell, “Generation of Light in Media with a Random Distribution of Nonlinear Domains,” Phys. Rev. Lett. **97**, 013902 (2006). [CrossRef] [PubMed]

**34. **M. Centini, V. Roppo, E. Fazio, F. Pettazzi, C. Sibilia, J. W. Haus, J. V. Foreman, N. Akozbek, M. J. Bloemer, and M. Scalora, ”Inhibition of linear absorption in opaque materials using phase-locked harmonic generation,” Phys. Rev. Lett. **101**, 113905 (2008). [CrossRef] [PubMed]

**35. **M. Scalora and M. E. Crenshaw, “A beam propagation method that handles reflections,” Opt. Commun. **108**, 191–196 (1994). [CrossRef]

**36. **E. Soergel, “Visualization of ferroelectric domains in bulk single crystals,” Appl. Phys. B **81**, 729–752 (2005). [CrossRef]