Comparative analysis of ferroelectric domain statistics via nonlinear diffraction in random nonlinear materials

We present an indirect, non-destructive optical method for domain statistic characterization in disordered nonlinear crystals having homogeneous refractive index and spatially random distribution of ferroelectric domains. This method relies on the analysis of the wave-dependent spatial distribution of the second harmonic, in the plane perpendicular to the optical axis in combination with numerical simulations. We apply this technique to the characterization of two different media, Calcium Barium Niobate and Strontium Barium Niobate, with drastically different statistical distributions of ferroelectric domains. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (000.2700) General science; (190.0190) Nonlinear optics; (190.4400) Nonlinear optics, materials. References and links 1. I. Freund, “Nonlinear Diffraction,” Phys. Rev. Lett. 21(19), 1404–1406 (1968). 2. G. Dolino, J. Lajzerowicz, and M. Vallade, “Second-harmonic Light Scattering by Domains in Ferroelectric Triglycine Sulfate,” Phys. Rev. B 2(6), 2194–2200 (1976). 3. J. A. Armstrong, N. Bloembergen, J. Ducuing, and P. S. Pershan, “Interactions between light waves in nonlinear dielectrics,” Phys. Rev. A 127(6), 1918–1939 (1962). 4. M. M. Fejer, G. A. Magel, D. H. Jundt and R. L. Byer, “Quasi-Phase-matched Second harmonic Generation: Tuning and Tolerances,” IEEE Journal Quant. Electron. 28(11), 2631–2654 (1992). 5. V. Berger, “Nonlinear Photonic Crystals,” Phys. Rev. Lett. 81(19), 4136–4139 (1998). 6. M. H. Chou, K. R. Parameswaran, M. M. Fejer, and I. Brener, “Multiple-channel wavelength conversion by use of engineered quasi-phase-matching structures in LiNbO3 waveguides,” Opt. Lett. 24(16), 1157–1159 (1999). 7. N. G. R. Broderick, G. W. Ross, H. L. Offerhaus, D. J. Richardson, and D. C. Hanna, “Hexagonally Poled Lithium Niobate: A Two-Dimensional Nonlinear Photonic Crystal,” Phys. Rev. Lett. 84(19), 4345–4348 (2000). 8. Y. Sheng, J. Dou, B. Ma, B. Cheng, and D. Zhang, “Broadband efficient second harmonic generation in media with short-range order,” Appl. Phys. Lett. 91, 011104 (2007). 9. R. C. Miller, “Optical harmonic generation in single crystal BaTiO3,” Phys. Rev. A 134(5), A313–A1319 (1998). 10. 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(3), 505–512 (2009). 11. S. Kawai, T. Ogawa, H. S. Lee, R. C. DeMattei, and R. S. Feigelson, “Second-harmonic generation from needlelike ferroelectric domains in Sr0.6Ba0.4Nd2O6 single crystals,” Appl. Phys. Lett. 73(6), 768–770 (1998). 12. A.R. Tunyagi, M. Ulex and K. Betzler, “Noncollinear Optical frequency Doubling in Strontium Barium Niobate,” Phys. Rev. Lett. 90(24), 243901 (2003). 13. J. Trull, I. Sola, B. Wang, A. Parra, W. Krolikowski, Y. Sheng, R. Vilaseca, and C. Cojocaru, “Ultrashort pulse chirp measurement via transverse second-harmonic generation in strontium barium niobate crystal,” Appl. Phys. Lett. 106(22), 221108 (2015). Vol. 26, No. 2 | 22 Jan 2018 | OPTICS EXPRESS 1083 #309558 https://doi.org/10.1364/OE.26.001083 Journal © 2018 Received 25 Oct 2017; revised 28 Dec 2017; accepted 2 Jan 2018; published 11 Jan 2018 Corrected: 06 February 2018 14. A. S. Aleksandrovsky, A. M. Vyunishev, A. I. Zaitsev, A. A. Ikonnikov, and G. I. Pospelov, “Ultrashort pulses characterization by nonlinear diffraction from virtual beam,” Appl. Phys. Lett. 98(6), 061104 (2011). 15. B. Wang, C. Cojocaru, W. Krolikowski, Y. Sheng, and J. Trull, “Transverse single-shot cross-correlation scheme for laser pulse temporal measurement via planar second harmonic generation,” Opt. Express 24(19), 22210– 22218 (2016). 16. 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 Y. Kivshar, “Third-harmonic generation via broadband cascading in disordered quadratic nonlinear media,” Opt. Express 17(22), 20117–20123 (2009). 17. Y. Sheng, A. Best, H. J. Butt, W. Krolikowski, A. Arie, and K. Koynov, “Three-dimensional ferroelectric domain visualization by Cerenkov-type second harmonic generation,” Opt. Express 18(16), 16539–16545 (2010). 18. Y. Le Grand, D. Rouede, C. Odin, R. Aubry, and S. Mattauch, “Second-harmonic scattering by domains in RbH2PO4 ferroelectric,” Opt. Commun. 200(1-6), 249–260 (2001). 19. V. Roppo, W. Wang, K. Kalinowski, Y. Kong, C. Cojocaru, J. Trull, R. Vilaseca, M. Scalora, W. Krolikowski, and Y. Kivshar, “The role of ferroelectric domain structure in second harmonic generation in random quadratic media,” Opt. Express 18(5), 4012–4022 (2010). 20. M. Ayoub, J. Imbrock, and C. Denz, “Ferroelectric domain diagnostics near the phase transition by Čerenkov second-harmonic generation,” Opt. Mater. Express 7(9), 3448 (2017). 21. M. Ayoub, J. Imbrock, and C. Denz, “Second harmonic generation in multi-domain χ(2) media: from disorder to order,” Opt. Express 19(12), 11340 (2011). 22. Th. Woike, T. Granzow, U. Dörfler, Ch. Poetsch, M. Wöhlecke, and R. Pankrath, “Refractive indices in congruently melting Sr0.61Ba0.39Nb2O6,” Phys. Status Solidi 186, R13–R15 (2001). 23. M. Eßer, M. Burianek, P. Held, J. Stade, S. Bulut, C. Wickleder, and M. Mühlberg, “Optical characterization and crystal structure of the novel bronze type CaxBa1-xNb2O6 (x = 0.28; CBN-28),” Cryst. Res. Technol. 38(6), 457– 464 (2003). 24. J. Trull, C. Cojocaru, R. Fischer, S. M. Saltiel, K. Staliunas, R. Herrero, R. Vilaseca, D. N Neshev, W. Krolikowski, and Y. S. Kivshar, “Second-harmonic parametric scattering in ferroelectric crystals with disordered nonlinear domain structures,” Opt. Express 15(24), 15868–15877 (2007). 25. E. Soergel, “Visualization of ferroelectric domains in bulk single crystals,” Appl. Phys. B 81(6), 729–751 (2005). 26. M. Scalora and M. E. Crenshaw, “A beam propagation method that handles reflections,” Opt. Commun. 108(46), 191–196 (1994).


Introduction
The spatial modulation of the linear refractive index offers extensive control of light wave propagation properties, providing changes in the propagation direction (diffraction gratings), reflection and transmission (Bragg mirrors) or molding the dispersion relation by using photonic crystals. Apart from modulating the linear susceptibility, the spatial modulation of the nonlinear coefficient has proven to effectively extend the mechanisms of light-matter interaction, and to be especially important for the field of nonlinear optics. In this regard, scattering of the second harmonic radiation obtained when a fundamental beam propagates inside ferroelectric crystals has been a centerpiece of nonlinear optics since its inception [1,2]. The idea of periodically inverting the sign of the quadratic nonlinear susceptibility, while keeping constant the refractive index, to enhance the nonlinear interaction efficiency, was proposed by Bloembergen in one of his seminal papers [3]. The so-called quasi phase matching (QPM) method was successfully implemented by M. Fejer [4], and since then its applications have grown.
It is well known that a one dimensional (1D) structure with periodic inversion of the quadratic nonlinearity sign can compensate the phase mismatch (Δk) of a particular nonlinear process. To this end, the spatial periodicity (Λ) has to fulfill the condition Λ = 2π/Δk. This process can be understood by assuming that the spatial periodicity introduces a reciprocal lattice vector (RLV) G = 2π/Λ that compensates the phase mismatch, so the condition Δk + G = 0 holds. This idea can be extended to more sophisticated nonlinear structures having 2D periodic modulation of the nonlinearity [5]. 2D QPM structures increase the available reciprocal lattice vectors and can simultaneously satisfy phase matching conditions for noncollinear propagation or multiple wavelengths. Later on, with the proposal of chirped QPM structures or quasi-periodic distributions an even broader RLV space became available [6][7][8].
In this spirit, structures displaying homogeneous linear properties, but a 2D random distribution of domains with inverted sign of quadratic nonlinearity, further increase the RLV space by providing a continuous set of reciprocal vectors that span around a given central value determined by domain size. Some as-grown ferroelectric crystals, as for example SBN and BaTiO3 possess this kind of domain distribution and provide phase matching over a broadband range of wavelengths [9, 10]. These crystals have been proposed as potential candidates for broadband frequency converters [11]. As-grown ferroelectric random crystals show a needle-like structure of domains oriented along the optical axis with a set of RLV confined in the plane perpendicular to this axis. When a fundamental beam propagates in this plane, the crystal emits a SH signal over the whole plane, as shown schematically [12]. The angular spreading of this planar SH emission is strongly dependent on the RLV distribution, which in turn depends on the particular domain size and distribution within the crystal. Any change in the domain pattern, which may be accomplished by either partial or total crystal poling, will be reflected in the transverse profile of the SH intensity. Planar SH generation in these crystals has been successfully used for different applications such as characterization of ultrashort optical pulses [13][14][15], cascaded THG [16], and nonlinear domain visualization [17].
In this work we present a combined experimental and numerical study of the angular distribution of the planar SH emission in different types of random media. We demonstrate that the emitted SH pattern variation as a function of fundamental wavelength provides relevant information about the statistics of the ferroelectric domain pattern in the crystal. We use a numerical model based on a 2D split-step fast Fourier transform beam propagation method to verify and complement the experimental observations. Finally, we show that our results are consistent with the direct analysis of domain statistics determined by means of near-field nonlinear microscopy.

Models and experimental results
Second harmonic scattering from materials with a random distribution of the quadratic nonlinearity has been studied and reported by different authors. Among them, Dolino et al.
[2] showed that the intensity of the SH emission should be proportional to the square modulus of the Fourier transform of the nonlinear polarization distribution in the material, which corresponds to the RLV distribution provided by the structure. This indicates that the emission of SH in a particular direction within the crystal will take place if RLV compensate for the phase mismatch in that direction. Figure 1(b) shows schematically how this phase mismatch compensation takes place.
A perfectly homogeneous distribution of G vectors spanning all possible directions in the plane should give rise to a homogeneous planar SH signal as shown in Fig. 1(a). Any deviation from this ideal frame should be reflected in the angular dependence of the observed SH pattern. The RLV modulus value, G, needed to compensate the phase mismatch at a given fundamental wavelength and at a given emission angle may be obtained readily from the scheme shown in Fig. 1 where kω and k2ω are the wave vectors of the fundamental and SH, respectively, and θ corresponds to the external propagation angle of the SH signal with respect to the fundamental beam propagation direction.
Here, I(ω) is the intensity of the fundamental wave, deff the effective quadratic nonlinear coefficient which depends on material parameters and interaction geometry, and G = |Δk| is the absolute value of the phase mismatch. The last term represents the effect of disorder in the domain distribution [18]. This model constitutes the basis of many different studies of SH scattering by ferroelectric crystals with disordered domains [19][20][21]. The typical emission pattern corresponding to Eq.
(2) is shown in Fig. 2(a). A structure with well-defined average domain diameter ρo and small standard deviation gives rise to an emission pattern having two well-resolved maxima which correspond to the emission angles where the phase mismatch is compensated by the particular RLV provided by G  . For the typical domain size dispersion found in common ferroelectric crystals, as for example SBN, one can show that a change in wavelength results in a change of the angular emission because the G vector distribution remains unaltered (fixed by the crystal). However, the phase matching condition changes due to material dispersion. Figure 2(b) shows an increased emission angle when the wavelength increases.
In this work we study two different nonlinear crystals with disordered domain structure: a sample of Strontium Barium Niobate (Sr0.61Ba0.39Nb2O6, SBN61) (5x5x10mm) and a sample of Calcium Barium Niobate (Ca0.28Ba0.72Nb2O6, CBN28) (5x5x5 mm). The CBN sample is a typical as-grown crystal with neither temperature nor electric field treatment. The SBN sample was exposed to a variable electric field to change its domain distribution in a nonprescribed way. Figure 3 shows a plot of the G vectors needed to compensate the phase mismatch for SBN (according to Eq (1)) at different angles of emission, when the fundamental wavelength between 800 nm and 1600 nm. For these calculations, we used the available data for SBN [22] and similar results are obtained for our second sample when using the refractive index of CBN [23]. Both crystals are transparent over a very broad wavelength range (typically between 350 nm and 4 μm) and can provide SH radiation over this whole wavelength range. This kind of broadband SHG is not possible in the usual crystals used in nonlinear optics, which provide PM through optical birefringence. The different colors in Fig. 3 represent the G modulus value needed to compensate phase mismatch at this particular external angle and wavelength. From these results we can assert that for the standard domain distributions both crystals exhibit basically the same behavior as far as the nonlinearity modulation-mediated second harmonic emission is concerned. Furthermore, as Eq.
(2) shows, for a crystal with a domain structure consistent with Le Grand model (i.e. narrow distribution of G vectors distributed around a particular central value Gm), the angle corresponding to the SH maxima would increase with wavelengths. However, as we show below, it is the specific feature of the actual domain distribution that strongly affects the observed spatial distribution of the emitted SH. In a first experiment we measured the SH pattern emitted by the CBN28 and SBN61 samples for two particular fundamental wavelengths: 800 nm (from a Ti:Sapphire oscillator with 1W average power, emitting 150 fs pulses at 76 MHz repetition rate) and 1064 nm (from a Nd:YAG laser of 200 mW average power, emitting 8 ns pulses at 10 Hz repetition rate). The collimated fundamental laser beam impinged on the crystal at normal incidence. The crystal optical axis (Z) was placed vertical such that the SH pattern was emitted in a horizontal plane. A rotating arm, centered on the crystal position allowed us to measure the SH as a function of the angle with respect to the direction of propagation of the fundamental beam. A photodiode with a 200-micron pinhole diaphragm was used to measure the SH signal in each particular direction. Ordinary or extraordinary polarization of the fundamental beam was selected using a half-wave plate placed in front of the crystal. The SH polarization is always extraordinary since the nonlinear susceptibility matrix of these crystals allows for both ee-e and oo-e interaction [24].
The experimental results are shown in Fig. 4, where we plot the angular distribution of the SH emitted at 400 nm and at 532 nm for an as-grown CBN28 crystal (sample 1) and an artificially poled SBN61 crystal (sample 2). These results clearly show a drastic difference in the angular SH intensity patterns for both crystals. In the CBN28 crystal the SH is emitted in a broad angular range up to 60°; the emission patterns for both fundamental wavelengths are similar. In contrast, the SBN sample leads to a much narrower SH angular emission with two clear maxima whose separation is sensitive to the fundamental wavelength. Moreover, an increasing wavelength causes these maxima to move in a direction opposite to that predicted by the Le Grand model illustrated in Figs. 2 and 3: longer wavelengths correspond to smaller emission angle. The SH angular distribution from the CBN (sample 1) represents the emission pattern from a medium with statistical domain and RVL distribution with large standard deviation (σ). As σ increases, the initially two well resolved SH maxima shown in Fig. 2 start to broaden and subsequently merge forming a single-peaked, broad intensity distribution similar to that of sample 1. A more detailed study of this effect can be found in [19]. However, the behaviour of the SBN61 crystal (sample 2) cannot be reproduced using the Le Grand model, thus highlighting its limits. It thus became evident that additional experiments with continuously varied input wavelengths were necessary to shed light on actual domain statistics in our SBN61 sample. To this end, we conducted an extensive study of the SH emission pattern in both crystals using a continuously tuned light source covering the spectral window between 800 nm and 1600 nm. Figure 5 shows a schematic representation of our experimental setup. We used two different tuneable laser sources: a femtosecond Ti:Sapphire oscillator (Coherent, Chameleon Ultra II) and an Optical Parametric Oscillator (APE GmbH, Chameleon Compact OPO) providing femtosecond pulses at wavelengths ranging from 800 nm to 1000 nm and from 1000 nm to 1600 nm, respectively. The laser beam propagates at normal incidence onto the sample surface (Sample1: as-grown CBN28, Sample 2: artificially poled SBN61). The polarization state and the average power of the incident fundamental beam were controlled by the combination of a half-wave plate and a polarizer. The output average power for all wavelengths and for both polarizations (ordinary and extraordinary) is fixed at the same value: 300 mW for CBN28 and 250 mW for SBN61. The polarization state of the input beam is selected either perpendicular for oo-e, or parallel for ee-e interactions, respectively. The incident laser beam is slightly focused in the crystal by means of a planoconvex lens (beam waist located in the middle of the crystal). The beam width at 1/e 2 in intensity was 400 μm for both crystals. The temporal pulse duration at FWHM was around 180 fs (pulse width vary form 140 fs to 200 fs due to its source wavelength dependence). We projected the emitted SH beam pattern on a diffusing screen (Lee Filter) and imaged the pattern into a CCD camera. A high pass optical filter (Thorlabs FGB-37-A) was used to block the strong IR fundamental beam before the CCD camera (Imaging Source DMK31BU). We used a stripe of millimetre paper attached to the screen to calibrate spatial dimension of the SHG profiles. The distance between screen and output facet of the crystal was measured by digital calliper. From the aforementioned experimental geometry, we were able to calculate angular intensity profiles.   The data depicted in Fig. 6 (and later in Fig. 8) was artificially zeroed near normal incidence to eliminate the experimentally observed strong collinear signal. For some wavelengths our laser emits additional weak spectral components, which were difficult to filter out because its wavelengths overlapped with the range of emitted SHG signals. Moreover, the SHG intensity emitted from the crystal was considerably weak, and for the most of the measurements, the spurious collinear signal exceeded the intensity of SHG signal at least by factor of two. These results show that the characteristic emission pattern remains almost undisturbed over the entire wavelength range with a weak signature of SH maxima appearing for the longest wavelengths. These results are compatible with those shown in Fig.  2(b) for a distribution with large standard deviation and consequently broad SH peaks. Since the peaks move apart with increased wavelength, one also expects the appearance of resolved maxima for larger wavelengths as experimentally observed.
In order to determine the most significant RLV participating in this SH emission we overlap these experimental results with the map shown in Fig. 3, calculated for CBN. In Fig.  7 we represent the width of the SH emission (at 1/e 2 level; indicated by the bar) and the position of the SH peaks (indicated by the circles). From this plot we assert that the crystal must possess a broad RLV distribution with values up to roughly 12 μm −1 . The domain size contributed to the SH peaks is around 2 microns, similar to previously reported as-grown samples [10]. The experimental results corresponding to the SBN61 crystal (sample 2) are summarized in Figs. 8 and 9. We can immediately see a drastically different angular emission of the second harmonic in this sample compared to CBN28. Firstly, the emission takes place in a much narrower angular range. Furthermore, these plots confirm the atypical tendency of the angular separation of second harmonic peaks to decreases with the wavelength corroborating the initial guess that a completely different domain distribution should be responsible for this emission.  In order to provide a meaningful interpretation of the results of sample 2 we need more information on the domain statistics. Previously reported results dealt with the evolution of the SH emission pattern under well-controlled conditions of poling of the samples. In our case, the poling action was performed in an uncontrolled manner, and we do not have reliable information about the domain distribution for sample 2 (SBN). To shed light into the question of domain statistics, we need a method that can give us an image of the ferroelectric domains pattern. Unfortunately, traditional ferroelectric domain visualization methods employ chemical etching and hence are destructive and do not allow the visualization of domains with a size below 1 μm [17,25]. As a solution, we used the SHG imaging microscopy technique. This method is based on scanning a tightly focused ultrashort pulsed beam (Ti:Sapphire oscillator (Coherent)) onto the sample. As the beam is scanned, the generated SHG signal from the sample is recorded, forming an image. Since the generation of the SH signal depends on the sign of the nonlinearity, this method finally allowed the visualization of the domains. The images of sample 1 and sample 2 are shown in Figs. 10(a) and 10(b), with 122 nm and 204 nm resolution respectively. In these images, we can identify very different domain distributions. For sample 1 (CBN): the inverted nonlinear domains with bright and dark colors are compactly arranged (high filling factor) within the crystal; the typical domain size is smaller than 1 μm but the agglomeration of domains leads to a broad distribution of domain sizes. This is in agreement with the well-known domain distribution in as-grown ferroelectric crystals. On the other hand, the artificially poled sample 2 (SBN) shows completely different domain statistics with bigger size and a low fraction inverted domains on otherwise homogeneous nonlinear background. This is precisely the nontrivial domain distribution responsible for the unusual SH emission pattern from the SBN sample. Based on these observations one can propose a possible mechanism responsible for the SH emission pattern. First, since the number of domains has decreased, one should expect a reduction in their contribution to the nonlinear SH diffraction process. Furthermore, as now there are large regions of the crystal with fixed sign of nonlinearity, one should expect increased contribution to the SH originating from these regions, especially since the phase mismatch decreases with longer wavelengths.
Inspection of the phase mismatch curve for the homogeneous SBN crystal effectively confirms this hypothesis as it shows that we are approaching the phase-matching condition as we increase the wavelength, the condition of PM being very close to 1600 nm for ooe interaction. This would appear in the emission pattern approaching forward emission as for uniform crystal.
As the PM condition is approached, it is reasonable to expect an increase of the SH signal at zero angle, while the effect of SH scattering is reduced by the decreased number of domains. Under these conditions one should expect the emitted SH angular pattern to follow the phase-mismatch curve of the material. Since this behaviour is intrinsically opposite to that expected from a fully random distribution (Le Grand et al. [18]), we can claim that the observation of the angular SH pattern as a function of wavelength is an effective method to obtain real-time information on the domain statistics of the nonlinear crystals.

Numerical simulations and discussion
A numerical simulation of the beam propagation and SH in these kind of samples was performed in order to test these hypothesis. A fundamental beam propagates in the y direction and SH diffraction emission occurs in the x-y plane. Since the ferroelectric domains have needle-like shape along the z axis, and the domain structure lies in the x-y plane, we adopt a 2D solution to simulate disordered nonlinear ferroelectric domains [19] based on the wellknown split-step fast Fourier transform (FFT) algorithm modified to handle the propagation of pulses of arbitrary duration. The split-step algorithm involves the "splitting" of free-space and material Eq.s that are generally of the first order. The material is described as an ensemble of nonlinear Lorentz oscillators excited by an external field and exhibiting a generic nonlinear polarization. Suffice it to say here that the use of FFTs allows the calculation of field derivatives to high accuracy, and that material equations are integrated using a modified, second-order-accurate, predictor-corrector algorithm. The interested reader will find more details on the method of integration in [19,26]. Regarding the structures, we assumed individual domains in the form of rods with a circular transverse profile; after selecting a mean domain size and a certain distribution and filling fraction, a given configuration is randomly generated and randomly placed in a rectangular area representing the size of the sample; the number of domains that takes part in the nonlinear process should be large compared to the actual size of the input optical beam. In the 2D domain pattern simulation for sample 1, domain size was assumed to have a Gaussian distribution with the experimentally retrieved mean diameter (Dmax = 2.0 μm). This assumption is based on the analysis of real domain image shown in Fig. 10(a). We have considered different values of the Gaussian distribution variance (σ) until we have obtained the same angular distribution as in the experimental results. As the total number of domains involved in the real experiment is much larger than the number of domains involved in this simulation, the implementation of a particular domain size distribution with broad variance in a single numerical realization was in general not possible. We conducted a number of different simulations and averaged over these domains realizations to achieve the desired σ distribution. This averaging is naturally performed in our real experiment during propagation of the fundamental beam along the 10 mm long crystals. Figure 11(a) depicts a typical generated 2D domain structure in the x-y plane, simulating sample 1. The white and black regions represent domains with opposite sign of the nonlinearity. Figure 11(b) is the corresponding Fourier spectrum of the domain structure in reciprocal space. The majority of reciprocal vectors are gathered around a broad ring with radius Gmax = π/Dmax = 1.5 μm -1 . To simulate sample 2, we have considered a domain distribution similar to the one revealed by the image shown in Fig. 10(b). The domains with positive χ (2) have a very low filling factor (around 6%) and they are placed on a homogeneous negative χ (2) background. We have simulated different Gaussian distributions around the diameter mean value of 1 μm until we have obtained the same SH angular distribution as in the experiment. The simulated 2D domain structure in the x-y plane is shown in Fig. 11(c) where a Gaussian distribution of domains with a mean size of Dmax = 0.9 μm has been assumed. Figure 11(d) is the corresponding Fourier spectrum of the domain structure in reciprocal space. The majority of reciprocal vectors are gathered around a narrow ring with radius Gmax = π/Dmax = 3.5 μm −1 .
A typical result of the propagation code is shown in Fig. 12 for a particular realization close to Le Grand type distribution similar to that shown in Fig. 11(a). Figure 12(a) represents the far field distribution in air (equivalent to the spatial transverse or angular spectrum) as a function of the propagation distance inside the crystal. The far-field angular distribution profile at the maximum propagation distance is plotted in Fig. 12(b). The SH emission in the form of well resolved maxima is clearly observed, as expected. Fig. 13. Comparison between the experimental results for as-grown CBN crystal (top) and numerical simulations (bottom) at different wavelengths for the domain distribution shown in Fig. 11(a) Averaging over 50 different realizations we obtained the simulation of the experimental results for sample 1, confirming the mean value and desired variance. Figure 13 shows a comparison between our simulated far-field angular distribution and the experimental results for both oo-e interaction (blue line) and ee-e interaction (red line). Simulation and experiment are in good agreement, showing that the process of averaging over different realizations can be an efficient way to simulate the large variance distribution.
The results for sample 2 show the combined effects of domain distribution and nonlinear background on the SH intensity angular distribution. From the far-field angular distribution results we surmised that the nonlinear domain distribution dominates the parametric process between 800 and 1200nm. The nonlinear background dominates the parametric process at longer wavelengths and it completely dominates at 1600 nm as we approach the phasematching condition. This is evident for the oo-e interaction when we are close to PM and a strong peak at 0° appears. However, because the coherence length for ee-e, LCeee, is far less than LCooe, for ee-e interactions the nonlinear background cannot completely dominate the parametric process even at 1600 nm. All parametric processes for ee-e interaction are a result of the combined effects of domain distribution and nonlinear background. Figure 14 shows a comparison between our simulated far-field angular distribution and the experimental results for both oo-e interaction (blue line) and ee-e interaction (red line). The simulation results are in good agreement with the experimental results.

Conclusions
We demonstrate an indirect non-destructive optical method for the characterization of nonlinear domain statistics based on the analysis of the second harmonic generation intensity angular distribution. The method is based on a combination of experimental measurements and numerical simulations. We begin by experimentally measuring a SH diffraction pattern from real crystals. Then we design 2D domain structures and numerically simulated the propagation through them using a spectral method, and found good agreement between our simulations and the experimental data. From these results we deduce that our domain designs are realistic. The method was implemented experimentally and tested in crystals with different types of domain distributions. The analysis of the dependence of the angular distribution of SHG on the fundamental beam wavelength informs us about complicated domain structures. This method could be used for real time monitoring of unknown domain distributions created by poling or crystal-growing processes.