Third-harmonic generation microscopy with Bessel beams : a numerical study

We study theoretically and numerically third-harmonic generation (THG) from model geometries (interfaces, slabs, periodic media) illuminated by Bessel beams produced by focusing an annular intensity profile. Bessel beams exhibit a phase and intensity distribution near focus different from Gaussian beams, resulting in distinct THG phase matching properties and coherent scattering directions. Excitation wave vectors are controlled by adjusting the bounding aperture angles of the Bessel beam. In addition to extended depth-of-field imaging, this opens interesting perspectives for coherent nonlinear microscopy, such as extracting sample spatial frequencies in the λ /8 λ range in the case of organized media. © 2012 Optical Society of America OCIS codes: (180.4315) Nonlinear microscopy; (190.4160) Multiharmonic generation; (170.3880) Medical and biomedical imaging. References and links 1. W. Mohler, A. C. Millard, and P. J. Campagnola, “Second harmonic imaging of endogenous structural proteins,” Methods 29, 97–109 (2003). 2. J.-X. Cheng, and X. S. Xie, “Coherent anti-Stokes Raman scattering microscopy: instrumentation, theory, and applications,” J. Phys. Chem. B 108, 827–840 (2004). 3. A. Volkmer, “Vibrational imaging and microspectroscopies based on coherent anti-Stokes Raman scattering microscopy,” J. Phys. D: Appl. Phys. 38, R59–R81 (2005). 4. D. Yelin, and Y. Silberberg, “Laser scanning third-harmonic generation microscopy in biology,” Opt. Express 5, 169–175 (1999). 5. J. Mertz, and L. Moreaux, “Second-harmonic generation by focused excitation of inhomogeneously distributed scatterers,” Opt. Commun. 196, 325–330 (2001). 6. J.-X. Cheng and X. S. Xie, “Green’s function formulation for third harmonic generation microscopy,” J. Opt. Soc. Am. B 19, 1604–1610 (2002). 7. D. Débarre, W. Supatto, and E. Beaurepaire, “Structure sensitivity in third-harmonic generation microscopy,” Opt. Lett. 30, 2134–2136 (2005). 8. C. Maurer, A. Jesacher, S. Furhapter, S. Bernet, and M. Ritsch-Marte, “Tailoring of arbitrary optical vector beams,” New J. Phys. 9, 78 (2007). 9. M. R. Foreman, S. S. Sherif, P. R. T. Munro, and P. Török, “Inversion of the Debye-Wolf diffraction integral using an eigenfunction representation of the electric fields in the focal region,” Opt. Express 16, 4901–4917 (2008). 10. C. Lutz, T. S Otis, V. DeSars, S. Charpak, D. A. DiGregorio, and V. Emiliani, “Holographic photolysis of caged neurotransmitters” Nat. Methods 5, 821–827 (2008). 11. N. Olivier and E. Beaurepaire, “Third-harmonic generation microscopy with focus-engineered beams: a numerical study,” Opt. Express 16, 14703–14715 (2008). #173834 $15.00 USD Received 3 Aug 2012; revised 3 Oct 2012; accepted 3 Oct 2012; published 16 Oct 2012 (C) 2012 OSA 22 October 2012 / Vol. 20, No. 22/ OPTICS EXPRESS 24886 12. K. Yoshiki, K. Ryosuke, M. Hashimoto, T. Araki, and N. Hashimoto, “Second-harmonic-generation microscope using eight-segment polarization-mode converter to observe three-dimensional molecular orientation,” Opt. Lett. 32, 1680–1682 (2007). 13. V. V. Krishnamachari and E. O. Potma, “Detecting lateral interfaces with focus-engineered coherent anti-stokes raman scattering microscopy,” J. Raman Spectr. 39, 593–598 (2008). 14. Y. Barad, H. Eisenberg, M. Horowitz, and Y. Silberberg, “Nonlinear scanning laser microscopy by third harmonic generation,” Appl. Phys. Lett. 70, 922–924 (1997). 15. M. Müller, J. Squier, K. R. Wilson, and G. J. Brakenhoff, “3D-microscopy of transparent objects using thirdharmonic generation,” J. Microsc. 191, 266–274 (1998). 16. D. Oron, D. Yelin, E. Tal, S. Raz, R. Fachima, and Y. Silberberg, “Depth-resolved structural imaging by thirdharmonic generation microscopy,” J. Struct. Biol. 147, 3–11 (2004). 17. C. K. Sun, S.-W. Chu, S.-Y. Chen, T.-H. Tsai, T.-M. Liu, C.-Y. Lin, and H.-J. Tsai, “Higher harmonic generation microscopy for developmental biology,” J. Struct. Biol. 147, 19–30 (2004). 18. D. Débarre, W. Supatto, A.-M. Pena, A. Fabre, T. Tordjmann, L. Combettes, M.-C. Schanne-Klein, and E. Beaurepaire, “Imaging lipid bodies in cells and tissues using third-harmonic generation microscopy,” Nat. Methods 3, 47–53 (2006). 19. N. Olivier, M. Luengo-Oroz, L. Duloquin, E. Faure, T. Savy, I. Veilleux, X. Solinas, D. Débarre, P. Bourgine, A. Santos, N. Peyriéras, and E. Beaurepaire, “Cell lineage reconstruction of early zebrafish embryos using labelfree nonlinear microscopy,” Science 339, 967–971 (2010). 20. J. Durnin, J. J. Miceli, and J. H. Eberly, “ Diffraction-free beams,” Phys. Rev. Lett. 58, 1499–1501 (1987). 21. J. Durnin, “Exact solutions for nondiffracting beams. i. the scalar theory,” J. Opt. Soc. Am. A 4, 651–654 (1987). 22. J. E. Durnin, J. J. Miceli, Jr., and J. H. Eberly, “Comparison of Bessel and Gaussian beams,” Opt. Lett. 13, 79–80 (1988). 23. D. McGloin and K. Dholakia, “Bessel beams: diffraction in a new light,” Contemp. Phys. 46, 15–28 (2005). 24. A. Boivin, “On the theory of diffraction by concentric arrays of ring-shaped apertures,” J. Opt. Soc. Am. 42, 60-64 (1952). 25. W. T. Welford, “ Use of annular apertures to increase focal depth,” J. Opt. Soc. Am. 50, 749–752 (1960). 26. C. J. R. Sheppard, “The use of lenses with annular aperture in scanning optical microscopy,” Optik 48, 329–334 (1977). 27. S. W. Hell, P. E. Henninen, A. Kuusisto, M. Schrader, and E. Soini, “Annular aperture two-photon excitation microscopy,” Opt. Commun. 117, 20–24 (1995). 28. E. J. Botcherby,R. Juškaitis, and T. Wilson, “Scanning two photon fluorescence microscopy with extended depth of field,” Opt. Commun. 268, 253–260 (2006). 29. N. Olivier, A. Mermillod-Blondin, C. B. Arnold, and E. Beaurepaire, “Two-photon microscopy with simultaneous standard and extended depth of field using a tunable acoustic gradient-index lens,” Opt. Lett. 34, 1684–1686 (2009). 30. T. A. Planchon, L. Gao, D. E. Milkie, M. W. Davidson, J. A. Galbraith, C. G. Galbraith and E. Betzig, “Rapid three-dimensional isotropic imaging of living cells using Bessel beam plane illumination,” Nat. Methods 8, 417– 423 (2011). 31. Z. Bouchal, J. Wagner, and M. Chlup, “Self-reconstruction of a distorted nondiffracting beam,” Opt. Commun. 151, 207–211, (2008). 32. F.O. Fahrbach, P. Simon, and A. Rohrbach, “Microscopy with self-reconstructing beams,” Nat. Photonics 4, 780–785, (2010). 33. K. Shinozaki, X. Chang-Qing, H. Sasaki, and T. Kamijoh, “A comparison of optical second-harmonic generation efficiency using Bessel and Gaussian beams in bulk crystals,” Opt. Commun. 133, 300–304 (1997). 34. J. Arlt, K. Dholakia, L. Allen, and M. J. Padgett, “Efficiency of second-harmonic generation with Bessel beams,” Phys. Rev. A 60, 2438–2441 (1999). 35. C. F. R. Caron and R. M. Potvliege, “Optimum conical angle of a bessel–gauss beam for low-order harmonic generation in gases,” J. Opt. Soc. Am. B 16, 1377–1384 (1999). 36. S. P. Tewari, H. Huang, and R. W. Boyd, “Theory of third-harmonic generation using bessel beams, and selfphase-matching,” Phys. Rev. A 54, 2314–2325 (1996). 37. V. E. Peet and R. V. Tsubin, “Third-harmonic generation and multiphoton ionization in bessel beams,” Phys. Rev. A 56, 1613–1620 (1997). 38. V. E. Peet and S. V. Shchemeljov, “Spectral and spatial characteristics of third-harmonic generation in conical light beams,” Phys. Rev. A 67, 013801 (2003). 39. S. Yang and Q. Zhan, “Third-harmonic generation microscopy with tightly focused radial polarization,” J. Opt. A 10, 125103 (2008). 40. B. Richards, and E. Wolf, “Electromagnetic diffraction in optical systems II. Structure of the image field in an aplanetic system,” Proc. Royal Soc. A 253, 358–379 (1959). 41. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge Univ. Press, 2006). 42. D. Débarre, N. Olivier, and E. Beaurepaire, “Signal epidetection in third-harmonic generation microscopy of turbid media,” Opt. Express 15, 8913–8924 (2007). #173834 $15.00 USD Received 3 Aug 2012; revised 3 Oct 2012; accepted 3 Oct 2012; published 16 Oct 2012 (C) 2012 OSA 22 October 2012 / Vol. 20, No. 22/ OPTICS EXPRESS 24887 43. Q. Zhan “Second-order tilted wave interpretation of the Gouy phase shift under high numerical aperture uniform illumination,” Optics Commun. 242, 351–360 (2004). 44. A. Piskarskas, V. Smilgevičius, A. Stabinis, V. Jarutis, V. Pašiškevičius, S. Wang, J. Tellefsen, and F. Laurell, “Noncollinear second-harmonic generation in periodically poled KTiOPO(4) excited by the bessel beam,” Opt. Lett. 24, 1053–1055 (1999). 45. V. V. Krishnamachari and E. O. Potma, “Multi-dimensional differential imaging with fe-cars microscopy,” Vib. Spectrosc. 50, 10–14 (2009). 46. D. Débarre and E. Beaurepaire, “Quantitative characterization of biological liquids for third-harmonic generation microscopy,” Biophys. J. 92, 603–612 (2007). 47. S. Carrasco, B. E. A. Saleh, M. C. Teich, and J. T. Fourkas, “Secondand third-harmonic generation with vector Gaussian beams,” J. Opt. Soc. Am. B 23, 2134–2141 (2006). 48. N. Olivier , F. Aptel, K. Plamann, M-C. Schanne-Klein and E. Beaurepaire, “Harmonic microscopy of isotropic and anisotropic microstructure of the human cornea,” Opt. Express 5, 5028–5040 (2010). 49. O. Masihzadeh, P. Schlup, and R. A. Bartels, “Enhanced spatial resolution in third-harmonic microscopy through polarization switching,” Opt. Lett. 34, 1240–1242 (2009). 50. F. Lu and W. Zheng and Z. Huang, “Coherent anti-Stokes Raman scattering microscopy using tightly focused radially polarized light,” Opt. Lett. 34, 1870–1872 (2009).


Introduction
Coherent nonlinear microscopy techniques such as second-harmonic generation (SHG), thirdharmonic generation (THG) or coherent anti-Stokes Raman scattering (CARS) microscopies are receiving considerable attention. These imaging modalities are all compatible with twophoton excited fluorescence (2PEF) microscopy and provide complementary information on sample intrinsic nonlinear optical properties with micrometer three-dimensional (3D) resolution [1][2][3][4]. One specificity of coherent nonlinear imaging techniques is that they are sensitive to both the sub-micrometer structure of the sample and the spatial distribution of the excitation laser field near its focus [5][6][7]. Indeed, the detected signal results from the coherent superposition in the detection plane of waves originating from different locations near focus, and the visibility of a particular distribution of scatterers is determined by interference phenomena. Engineering the spatial distribution of the intensity, phase and polarization of the focused driving field is possible by controlling the field at the pupil of the objective [8][9][10]. In coherent nonlinear imaging, focus shaping generally results in a modulation of phase matching and far-field scattering patterns [11][12][13], which potentially provides information about the sample microstructure.
In this article, we will focus on third-harmonic generation (THG) microscopy [14,15], an imaging modality which has proven useful for biological imaging [4,[16][17][18][19] and which is particularly sensitive to the interplay between sample and field structure. The most remarkable characteristic of THG microscopy with a Gaussian excitation beam is that no signal is obtained from a homogeneous positively dispersive sample [14]. This is due to the large phase mismatch accumulated between the fundamental and the harmonic beam across the focal region, owing to the Gouy phase shift of focused Gaussian beams. In contrast, a non-zero THG signal is obtained around χ (3) (−3ω; ω, ω, ω) inhomogeneities, with an efficiency depending on the relative axial sizes of the inhomogeneity and of the focused Gaussian beam [7].
Regarding focus shaping, we will concentrate on Bessel beams [20][21][22][23] which, as we will see, provide a convenient means to control the excitation wave vectors for coherent nonlinear processes such as THG. We will specifically consider Bessel beams as obtained by focusing an annular intensity distribution using a high numerical aperture (NA) objective lens. The use of an annular aperture in the back focal plane of a lens has been investigated even before Bessel beams were defined as such, and this arrangement was shown to increase the depth of field of an imaging system [24]. Since then, Bessel beams have received significant attention for microscopy applications. Their main characteristic for imaging applications is that they provide an extended depth of field without compromising the lateral resolution of the system, which has been explored in linear microscopy [25,26], as well as in 2PEF microscopy [27][28][29][30]. The fact that they exhibit "self-reconstructing" properties [31] also makes them attractive for imaging within complex media [32]. Finally, their phase and intensity distribution results in specific phase-matching mechanisms in nonlinear optics, which attracted some attention for SHG [33,34] and THG [35][36][37][38][39].
The use of Bessel beams for THG microscopy is therefore interesting for two reasons. Firstly they can provide an extended depth of field; secondly, they are expected to modify the phasematching mechanisms and therefore the contrast of THG microscopy images, and possibly provide additional structural information. Our aim in this article is to shed light on these processes in a simple context: we only consider isotropic media, and an incoming linear polarization. In section 2 we outline the geometry considered, we discuss basic properties of Bessel beams, and we investigate THG phase matching in the context of Gaussian and Bessel beams. In section 3 we subsequently use a vectorial numerical model adapted to the high NA focusing conditions used in microscopy [6,11] to study THG from axial slabs, interfaces, and axially periodic samples excited using tightly focused Gaussian and Bessel beams.

Model
The method used to calculate THG from arbitrary excitation fields and sample geometries is detailed in [11]. Briefly, the focal field distribution is first calculated using a Debye-Wolf diffraction integral [40], as further discussed in the next subsection. Because the signal emitted from an homogeneous medium is zero (see subsection 2.3), the non linear susceptibility is assumed to be zero outside the sample. The induced nonlinear polarization in the sample volume is then calculated for a given spatial distribution of the sample nonlinear susceptibility χ (3) . The resulting nonlinear field is finally propagated into the far field using Green's functions [6,41], and signal levels and scattering patterns are analyzed. The notations used in the article are summarized in Fig. 1. Apart from the case of imaging axially periodic media, we will consider THG in the forward direction only (|θ T HG | < π/2) since backward THG is in most cases orders of magnitude weaker than forward THG because of larger phase mismatch [42].

Non-paraxial focusing of Bessel beams
For reference, we remind the reader that a perfect (theoretical) Bessel beam would be obtained by focusing an infinitely narrow transverse intensity distribution and can be described as: Such an ideal Bessel beam is non-diffractive (i.e. the intensity does not depend on z) and is characterized by a rapid on-axis phase variation.
We consider now experimentally more realistic Bessel beams produced by propagating an incident x-polarized Gaussian beam through an annular aperture at the back aperture of a high NA objective, and we neglect diffraction at the aperture edges. The phase and amplitude at the back aperture are supposed to be constant over the annulus. Such beams can be described using two parameters: NA B and ΔNA with NA B = n.sin(θ B ) and ΔNA ≈ 2n.cos(θ B ).Δθ , where θ B depends on the radius of the annular aperture, and ΔNA represents the width of this aperture. The focal field E f can be expressed as a function of θ B and Δθ as: with: Here J 0 , J 1 and J 2 are the zero, first and second order Bessel functions. We note that a Gaussian beam with numerical aperture NA is obtained when setting Δθ = θ B = 1 The most important phenomenon for scanning microscopy is the increase in axial extension, which makes Bessel beams a prime target for extended depth of field imaging [27,28]. This is particularly clear when comparing the Gaussian beam with the fourth Bessel beam in Fig. 2A (NA B = 0.95 ± 0.05), for which the axial extension is increased by a factor ≈ 10 for a similar maximum NA. The secondary rings are generally an issue in linear microscopy because they generate considerable background. However their influence on THG contrast is expected to be more limited due to the nonlinear dependence of the signal generation on the excitation intensity. In addition, since the lateral extension of the third-order point spread function (PSF) is mostly determined by the principal excitation focus, its reduced lateral extension is a side benefit corresponding to a gain in the lateral resolution of the image. This improvement comes at the cost of energy being redistributed in the secondary rings, implying that higher excitation power may be required experimentally. Figure 2B summarizes these considerations by showing the intensity full width at half maximum (FWHM) in z and in x as a function of NA B and ΔNA. Whereas the lateral resolution is principally determined by NA B , the axial extension strongly depends on ΔNA, as can be seen on the two graphs showing the lateral and axial FWHM as a function of ΔNA for NA B = 0.95. For a given lateral resolution, it is therefore possible to increase the axial extension by decreasing ΔNA. Additionally, a crucial parameter for a coherent nonlinear effect such as THG is the spatial distribution of the phase. In line with previous reports [43], Fig. 2D shows that the axial phase gradient exhibited by the beam near its focus increases as NA B is increased, and increases weakly as ΔNA is decreased. This phenomenon results in different phase-matching properties for Bessel beams as compared to Gaussian beams.

Phase-matching mechanisms
It has been previously noted that third-harmonic generation by a Bessel beam can be described as a Bessel-Bessel interaction [36]. In the case of a perfect Bessel beam, the phase matching condition reads: where k ω (resp. k 3ω ) is the wave vector at the fundamental (resp. harmonic) frequency, i.e. |k ω | = 2πn ω /λ . The phase matching condition is thus satisfied in the case of a non-dispersive medium for which n ω = n 3ω along a THG scattering cone with angle θ T HG = θ B . As a result, THG is obtained from an homogeneous medium. This is in contrast with the case of Gaussian excitation, which in the case of a homogeneous medium and a Gaussian-Gaussian interaction yields: where k g,ω (resp. k g,3ω ) is a vector oriented along the z axis opposite to the beam propagation, which describes the Gouy phase shift around the focus of the excitation (resp. scattered) beam plotted on Fig. 2D [6]. In the case of a Gaussian beam with confocal parameter b, |k g,ω | ≈ π/4b. Because of this additional wave vector, phase matching is only achieved in media with negative dispersion, and no THG is obtained in the case of non dispersive or positively dispersive media. Instead, when considering samples with infinite transverse (x − y) extension, THG is obtained from interfaces and objects whose longitudinal size is of the order of the coherent construction length (or effective coherence length) l f w defined as: When considering THG from a non-perfect Bessel beam produced by focusing an annular intensity distribution of finite angular spread ΔNA = 0, a similar behaviour is obtained. In this situation, THG can be described as a non-perfect Bessel -Bessel interaction and the phase matching condition follows Eq. (4), with (see derivation in Appendix): where ΔNA T HG is the THG angular spread. From these equations, it can be seen that no bulk emission is obtained in the case of positive dispersion, or of zero dispersion unless k g,ω = k g,3ω = 0, which corresponds to the perfect Bessel -Bessel interaction for which ΔNA = ΔNA T HG = 0. Instead, coherent emission is obtained from structures with longitudinal size of l f w or less.
To further understand the mechanisms that govern THG emission in media with positive dispersion, two limits should be considered. First, if dispersion is zero or negligible (i.e. (k 3ω − k ω ) << k g ), emission is dominated by the phase mismatch due to tripling of the excitation Gouy phase shift. The longitudinal effective coherence length is then determined by the numerical aperture and angular spread of the excitation beam. Conversely in the case of large positive dispersion or Bessel beams with very long axial extension (i.e. (k 3ω − k ω ) >> k g ), emission is dominated by the material index dispersion between the fundamental and harmonic frequencies (n 3ω − n ω ), and depends little on the excitation beam parameters (numerical aperture, angular dispersion).
Additionally, following [36], in the case of a Bessel excitation and a laterally infinite medium, phase matching should be investigated both in the lateral (x and y) and in the longitudinal (z) direction. Lateral phase matching appears through the integral: This integral is zero for k 3ω sin(θ T HG ) > 3k ω sin(θ B ), so that emission can only be obtained for sin(θ T HG ) ≤ sin(θ B )n ω /n 3ω , thereby preventing longitudinal phase matching for n 3ω > n ω . This function also reaches a maximum for k 3ω sin(θ T HG ) = k ω sin(θ B ), which corresponds to sin(θ T HG ) = sin(θ B )n ω /(3n 3ω ). This situation has been termed self-phase matching [36] since the excitation and emission Bessel beams have identical lateral profile, thus maximizing coherent coupling. The self-phase matching angle θ T HG does not however simultaneously satisfy the longitudinal phase matching condition for non dispersive or positively dispersive media. Coherent construction of the signal for the self-phase matching angle is thus limited to a finite length, so that the magnitude of the signal created at this angle oscillates as a function of sample thickness.
In summary, the THG produced by a perfect or imperfect Bessel beam is mostly scattered along two directions : θ T HG ≈ θ B , which satisfies the longitudinal phase-matching condition; and θ T HG ≈ θ B /3, which satisfies the lateral self-phase matching condition.  The first sample geometry that we consider in this study is that of a slab of variable width (e) with an infinite extension perpendicularly to the beam propagation direction (z). This geometry has been extensively studied in the case of Gaussian beams. The dependence of THG on slab width e gives straightforward information about the phase-matching efficiency along the zaxis [6,7,11].

THG from slabs
For this set of calculations, we neglect dispersion effects, i.e. we assume that n ω = n 3ω = n. We obtain therefore k 3ω = 3k ω , and the effective coherence length along the optical axis for a Gaussian-Gaussian interaction (or along θ T HG = θ B for an imperfect Bessel-Bessel interaction) can be expressed as: This length is a good approximation of the slab width giving maximum THG. In the case of a Gaussian excitation, as the thickness e increases from zero, we first expect a coherent increase of the signal (THG scales as the number of emitters squared, i.e. as e 2 for a slab sample) until the effective coherence length is reached, and then a decrease as a function of e because of destructive interferences until the THG signal reaches zero for a homogeneous sample, as illustrated in Fig. 3A (black curve).
The same holds for a Bessel excitation, for which the decrease in THG intensity occurs for thicker samples due to the axial spread of the Gouy phase shift as ΔNA decreases. In addition, THG exhibits oscillations as a function of the slab thickness, which correspond to the oscillations of the self-phase matched and on-axis emission with periods ≈ k(cos(θ B /3) − cos(θ B )) and ≈ k(1 − cos(θ B )), respectively. The oscillations of the on-axis scattered component are illustrated on Fig. 3B. We normalized the simulations by considering a constant power at the back aperture of the objective. While the expected THG signal is significantly decreased (by a factor ≈15 for NA B = 0.95 ± 0.05) when using Bessel beams, due to energy being redistributed in the secondary rings, it is by no means negligible.
The scattering patterns for different slab thicknesses are investigated in more detail in Fig. 4. Figure 4A illustrates THG scattering from a slab illuminated by a Bessel beam (NA B = 0.95 ± 0.05) as a function of slab width e and scattering angle. Three regimes are highlighted: on-axis scattering (red), intermediate-angle scattering (orange) and off-axis scattering (green) corresponding to three detection NA ranges, namely 0-0.2, 0.2-0.4 and 0.8-1. The integrated signal over these detection ranges is shown in Fig. 4B, along with the signal detected over the whole NA range (0-1.33). These different ranges can be linked to the three phase matching situations as discussed in the previous section: indeed in the absence of dispersion, THG scattering mostly occurs along the longitudinal (axial phase matching, θ T HG = θ B , off-axis scattering, green) and lateral (self-phase matching, θ T HG ≈ θ B /3, intermediate-angle scattering, orange) phase matching directions. Along the self-phase matching direction, we observe as expected dampened oscillations with a larger amplitude and period than along the on-axis direction. In contrast, at large angles we obtain emission from large e, indicating that the effective coherence length is larger than the slab widths investigated here. The diminishing total THG signal observed in Fig. 4A (black) on the considered range of slab thicknesses is therefore mostly due to the poorly phase-matched directions (red, orange). The signal created around θ T HG = θ B (green curve) is however also expected to decrease as the slab width is further increased, on a length scale given by the imperfect Bessel-Bessel coherence length (see Eqs. (5)-(7)), consistently with previous predictions [36].

Influence of dispersion
We now consider the influence of linear dispersion, i.e. we assume that n ω = n 3ω . Indeed, while we focus mostly on non-dispersive media in this study to simplify the discussion of phase-matching mechanisms, most biological samples are positively dispersive (typically Δn ≈ 0.02 − 0.03) and dispersion cannot be neglected in practice. In the case of THG with a tightly focused Gaussian beam, phase-matching is dominated by the Gouy phase mismatch, meaning that dispersion has only a minor effect. However, we have seen that in the case of Bessel-Bessel coupling, the geometrical phase matching conditions are more favorable and the interaction length is larger, so the influence of dispersion may become significant when the medium coherence length defined as l c = π/Δn is shorter than the dispersion-free imperfect Bessel-Bessel interaction length. Figure 5 shows the THG signal as a function of the width of an axial slab for several values of n 3ω − n ω for a Gaussian beam (Fig. 5A, NA = 1.2) and a Bessel beam (Fig. 5B, NA B = 0.95 ± 0.05). As expected, THG is obtained in the case of a homogeneous negatively dispersive medium, both for the Gaussian and for the Bessel excitation, but this bulk emission disappears for zero or positively dispersive media. Nevertheless, the influence of dispersion on the curves of THG as a function of the slab thickness differs significantly for the two excitation geometries. In the case of the Gaussian excitation, apart from the bulk emission efficiency, the THG dependence on the sample size does not change much with dispersion. In particular, the slab thickness for which the maximum THG emission is obtained is little affected. On the contrary, the location of the maximum of the signal curve obtained for a Bessel excitation depends strongly on Δn. This is a straightforward consequence of the change in the non-dispersive effective coherence length of the beams compared to the dispersion coherence length l c = π/Δn. In particular, increasing the positive dispersion of the nonlinear medium reduces the slab width for which the maximum THG signal is obtained in the case of a Bessel excitation. This can be further understood by considering the scattering diagram for zero and positive dispersion (Fig. 5C). As described in section 2.3, THG propagation angles cannot exceed θ B for a slab perpendicular to the propagation axis, and dispersion results in reduced phase-matching for waves scattered at this angle. This component therefore exhibits a bell-shaped variation as e increases, corresponding to the first period of an oscillating behavior with period l c = π/Δn. Intermediate-angle scattering retains its oscillating behavior in the presence of dispersion, although the frequency and amplitude of the oscillations are shifted.

Quasi phase-matching Bessel beams
In light of the calculation results obtained for the case of homogeneous slabs, we now consider THG from longitudinally periodic media illuminated by Bessel beams. We have seen that phase matching can be simply described in the case of Bessel excitation beams thanks to the limited variation of both phase gradient and intensity over long distances. Therefore Bessel beams are good candidates for quasi-phase matching (QPM) using axially periodic materials having a spatial period equal to twice the effective coherence length, which can compensate for the phase mismatch and produce efficient coupling between the excitation and harmonic beams over large propagation distances. For example, quasi-phase-matched SHG of Bessel beams has been experimentally demonstrated [44].
We consider in this section an axially periodic medium where χ (3) oscillates along z with a period e. For simplicity, we neglect the changes in linear indices. Figure 6 illustrates the calculated THG signal as a function of the spatial period for the same Bessel and Gaussian beams considered in Fig. 3. We point out several phenomena: 1. In each case, there is a resonant spatial frequency providing QPM, which is narrower for Bessel beams than for Gaussian beams, and gets narrower as the Bessel NA spread (ΔNA) gets smaller.
2. The enhancement of the total integrated signal for a given axial period (Fig. 6B) mostly stems from the resonance at the self-phase matching angle (θ T HG ≈ θ B /3), for which (quasi)-phase matching is then achieved both in the longitudinal and axial directions.  4. In the case of Bessel beams, the signal for large axial periods goes to zero only for large sample thicknesses, consistent with the calculations shown in Fig. 4A.
5. Since QPM can be achieved at the self-phase matching angle for which the secondary emission ring contribute efficiently to THG emission, the increased QPM efficiency yields a signal at resonance comparable to the Gaussian case.
It should also be noted that although the enhancement of the THG signal is larger for the self-phase matching angle, additional resonances occur for all emission angles for other axial periods. On-axis emission, in particular, is significant for a slightly smaller axial period.

Probing organized media using Bessel beams
This narrow near-axis quasi phase-matched scattering obtained with Bessel beams gives us one interesting possibility: probing organized media using these spatial resonances. We consider here the same geometry as in the previous paragraph, but assuming NA d,Max = n sin(θ d,Max ) = 0.2 to retain only the near-axis propagating THG (which would be simple to implement experimentally by using a low NA collection lens or adding a diaphragm after the collection optics), and we investigate the influence of excitation NA. Figure 7B illustrates the THG signal obtained from a axially periodic structure for several Bessel beams characterized by different numerical apertures and similar ΔNA. As NA B increases, the on-axis coherence length decreases (see section 3.1) so the spatial resonance period also decreases. Therefore, a very interesting property of Bessel beams is that they provide a way to probe the axial organization of a sample, since the spatial resonance period can be tuned from ≈ 500nm to ≈ 3μm by changing the excitation NA. We note that the resonance broadening observed for smaller NAs comes from the fact that we considered a constant ΔNA, and that it could be avoided by setting ΔNA ∝ NA B . In order to understand this dependence, we can express the on-axis wavevector mismatch mentioned in section 3.1 as: From this expression, we can predict the QPM period: This equation is plotted in black in Fig. 7E along with the QPM period calculated from the numerical simulations, and provides a good approximation of the resonance period. This confirms that THG efficiency for simple geometries (slab, periodic structure) can be described using simple phase-matching considerations, owing to the relative homogeneity of the intensity and phase gradient of the Bessel beams in the focal region. Similarly, the QPM period for the self-phase matching angle can be obtained as:  (12) which would in the same conditions correspond to approximately the same range of axial periods (blue line in Fig. 7E).
If the detection NA d,Max is increased in order to also detect the self-phase matched contribution, the signal is increased by a factor ≈ 4, and the resonance range is slightly shifted and broadened. This is illustrated by the the dashed curves in Fig. 7B and by the difference between the two theoretical resonant spatial periods in Fig. 7E , shown as black and blue curves). Finally, the use of higher detection NA would result in the detection of additional off-axis components and in a more complex signal (as in Fig. 6). We now consider backward QPM conditions, corresponding to the exact same sample geometry, but with a THG signal detected in the backward direction. We can see from Fig. 8D that we again observe spatial resonances, and that this resonant period increases with the NA considered.
From the geometrical phase-matching conditions illustrated in Figs. 8B and 8C, we can estimate the on-axis backward QPM resonance period as: which is in good agreement with the maximum resonant period that can be determined from Fig. 8D.
To estimate the range of axial periods producing backward THG scattering, we can consider the following two limit cases. First for a null NA B , the Bessel beam turns into a plane wave for which the coherence length is l c = λ /12n, and hence p = λ /6n. In the other limiting case of NA B ≈ n, the axial component of the excitation wave vector becomes negligible and we have: l c = λ /6n, and hence p = λ /3n.
If we now consider both forward and backward resonances, we can see from Eqs. (11) and (13) that we can probe periodicities ranging from p = λ /6n (low NA Bessel, backward scattering) to almost infinitely large periods. What limits this range is the geometrical focusing requirement: even at very high NA (e.g. 1.27 in water), the angle is still limited to 72.7 degrees which means that cos(arcsin(NA B /n)) can only vary from ≈ 0.3 to 1, excluding from the attainable range the region [λ /3n(1+0. 3) λ /3n(1−0.3)] ≈ [λ /4n λ /2n]. One possible workaround for this limit could be to use two different excitation wavelengths with non-overlapping exclusion ranges. While measuring large periods using QPM seems of limited interest, the possibility to measure periodicity in the range [λ /6n λ /4n] gives access to sub-wavelength information in ordered media.

Comparison of the THG imaging properties of focused Gaussian and Bessel beams
We finally summarize and compare the imaging properties of a Gaussian beam of NA = 0.7 and a Bessel beam defined by NA B = 1.15 ± 0.05. These two beams have similar axial extensions and therefore have fairly similar properties for incoherent imaging, and this axial extension is small enough that it may be used for microscopy imaging.  Figure 9 illustrates the difference between these two beams for THG. Figure 9A shows the phase and intensity distribution, with the main differences being the lateral extension and the phase distribution of the beams. Figure 9B displays the difference in THG signal as a function of a slab width, with the Gaussian beam providing maximum signal for a slab of ≈ 2.5μm, corresponding to the effective coherence length in the forward direction, while the THG produced by a focused Bessel beam exhibits a much slower decay because of the favorable off-axis phase-matching conditions. Figure 9C illustrates the average emission angle obtained from the previous geometry, and shows that THG scattering from a Gaussian beam resembles a Gaussian beam while THG from a Bessel beam is close to a Bessel beam. The last panel (Fig. 9D) shows the axial (NA d,Max = 0.2) forward-QPM of the two beams with a period of twice the coherence length in both cases, and a much sharper spatial resonance when using Bessel beams due to the larger interaction distance.

Conclusion and outlook
We have shown that THG from Bessel beams produced by focusing an annular intensity distribution with a high NA lens can be described using simple phase-matching considerations, which can be used to predict the signal expected from model geometries. THG with Bessel beams is interesting from two perspectives: 1. Rapid probing of large volumes, as in the case of extended-depth (incoherent) two-photon imaging.
2. Structural χ (3) microscopy with modified phase matching conditions. Bessel beams have very interesting characteristics for THG, such as the ability to probe organized structures.
We also point out that the combination of focus engineering and spatially selective detection is an interesting perspective, as already demonstrated in CARS [45], as it allows to separate different contributions, each providing different structural information. The methods explored in this article are generally transposable to other coherent imaging modalities.
Two effects have been neglected in this study which may be of importance in practical situations: 1. We considered optimal beam shaping by normalizing the power at the back focal plane of the objective. Even if there are efficient strategies to produce annular illumination before the focusing objective (e.g. [28]), some excitation energy will be lost in the shaping process, which is an important issue for a third-order process. Moreover, except in the case of QPM, the energy redistribution in the secondary rings implies a reduced nonlinear signal.
2. We only considered homogeneous media and neglected differences in optical indices. This model obviously needs to be refined in the case of extended depth-of-field imaging of heterogeneous media.
The increased depth of field conjugated with the increased coherence length open interesting possibilities for THG imaging. One is the possibility to image very large volumes, and to obtain quantitative information about the distribution of χ (3) . For example, since lipid droplets are one important source of contrast in biological tissues [18,46], THG with Bessel beams could be used to rapidly probe lipid contents in thick samples. Another perspective would be to measure unknown χ (3) more easily than with interference techniques [46], as the expected signal would be almost independent of the structure. However, the depth of field extension would of course reduce the sectioning ability of the microscope.
One perspective motivating the study of focal field engineering for coherent nonlinear microscopy is that the scattering patterns reflect the interplay between the (unknown) sample structure and a known driving field distribution. The phase and polarization properties of shaped focused beams have started to be investigated in regard to their consequence for coherent nonlinear imaging [13,45,[47][48][49][50]. Bessel beams appear as good candidates for harmonic imaging of organized media, because they provide a flexible means to control the excitation wave vectors.
To evaluate the variations of k g with NA B and ΔNA, we calculate the electric field distribution on axis (ρ = 0) around the focal point (z << 1/|k g |). Because we are only interested in a scaling law, we use a number of approximations that are detailed below. We have: In good agreement with more accurate simulations, we find that the sinc term that governs the intensity distribution corresponds to an axial width that varies as ∝ 1/(sin(θ B )Δθ ) ∝ 1/(NA B × ΔNA). The phase term should be evaluated after subtraction of the propagation phase term, which scales as ≈ e ikz(cos(θ B −Δθ )) . It is indeed easily seen that this phase term correctly becomes ≈ e ikz for θ B = Δθ , which corresponds to the Gaussian case. The remaining on-axis phase therefore becomes : e ik g z ≈ e ikz(cos(θ B )−cos(θ B −Δθ )) ≈ e −ikz sin(θ B )Δθ , (14) and hence :