A flexible, fast and benchmarked vectorial model for focused laser beams
Qingfeng Li, Maxime Chambonneau, Markus Blothe, Herbert Gross, Stefan Nolte
11 A flexible, fast and benchmarked vectorial model forfocused laser beams Q INGFENG L I , M AXIME C HAMBONNEAU , M ARKUS B LOTHE , H ERBERT G ROSS , AND S TEFAN N OLTE Institute of Applied Physics, Abbe Center of Photonics, Friedirich-Schiller-University Jena, Albert-Einstein-Str. 15, 07745 Jena, Germany Fraunhofer Institute for Applied Optics and Precision Engineering, Albert-Einstein-Str. 7, 07745 Jena, Germany * Corresponding author: [email protected] March 1, 2021
In-bulk processing of materials by laser radiation has largely evolved over the last decades and still opensup new scientific and industrial potentials. The development of any in-bulk processing application relieson the knowledge of laser propagation and especially the volumetric field distribution near the focus.Many commercial programs can simulate this, but, in order to adapt them, or to develop new methods,one usually needs to create a specific software. Besides, most of the time people also need to measurethe actual field distribution near the focus to evaluate their assumptions in the simulation. To easilyget access to this knowledge, we present our high-precision field distribution measuring method andrelease our in-house software InFocus [1], under the Creative Commons 4.0 License. Our measurementsprovide 300-nm longitudinal resolution and diffraction limited lateral resolution. The in-house softwareallows fast vectorial analysis of the focused volumetric field distribution in the bulk. The simulationsof light propagation under different conditions (focusing optics, wavelength, spatial shape, propagationmedium) are in excellent agreement with propagation imaging experiments. The aberrations provoked bythe refractive index mismatch as well as those induced by the focusing optics are both taken into account.The results indicate that our proposed model is suitable for the precise evaluation of energy deposition.
1. INTRODUCTION
In the last two decades, laser processing in the bulk of opticalmaterials has attracted intensive attention in a wide range ofacademic researches and industrial engineering. In-volume laserdirect writing enables precise three-dimensional structuring andhas allowed many innovative applications that include the fabri-cation of channels [2–5], waveguides [6–10], gratings [11], datastorage [12], and photonics quantum gates [13, 14]. The natureof the in-bulk processing also innovates new manufacturingprocedures such as bonding [15–19] and dicing [20] of brittle ma-terials. Among all potential applications, the precise descriptionand control of the laser focusing and the energy deposition arecrucial. Numerous experimental and theoretical investigationson the laser propagation and the energy absorption have beenintensively carried out [21–25].The propagation of the electromagnetic (EM) field can be rig-orously described by the finite difference time domain (FDTD)method [26], however in general it requires significant com-putational resources. Considerable effort has been devoted todesigning propagation equations that on one hand preserve theircomputational simplicity and on the other hand preserve thecorrect description of nonparaxial and vectorial effects. In the nonlinear propagation regime, only until recently, the unidirec-tional Hertz vector propagation equation (UHPE) [27, 28] wasderived to provide a seamless transition from Maxwell’s equa-tions to the various envelope-based models, which considerablyreduces the computational time. Simulations of the UHPE, re-quire starting from input conditions, i.e., from the Hertz vectorin a plane z = z . When the focusing elements have a high nu-merical aperture (NA), the input conditions are then determinedby a phase correction to the field that simulates the action ofthe focusing element. The input conditions for the UHPE wereconstructed by a detailed calculation of diffraction by vectorialdiffraction integrals (VDIs) [29, 30], i.e. the linear propagationmodel. However, even though the vectorial effects have beenconsidered by VDIs, the residual aberrations of the focusingelements are often ignored, which leads to a deviation from thecorrectness in the real laser processing conditions.The purpose of this work is to accurately describe these inputconditions by taking into account all the potential aberrationsthat may occur in the laser in-bulk focusing, and to allow for afast analysis with a proper transformation of the VDIs.One of the most widely used integral for analyzing the vecto-rial diffraction is the Debye-Wolf integral. As demonstrated by a r X i v : . [ phy s i c s . op ti c s ] F e b Leutenegger et al. [31], the 3D vectorial field distribution at thefocus can be computed plane by plane under a proper transfor-mation of the original Debye-Wolf integral. At a given axial posi-tion, the EM field in this plane is obtained by a two-dimensionalFourier transform. Lin et al. [32] have also demonstrated thatthis method is applicable to focusing through an interface be-tween two media of mismatched refractive index. In this paper,we further adapt this method to the real lens conditions and pro-vide a fast analysis tool for the evaluation of the actual EM fielddistribution at the focus of the lens whose residual aberrationscannot be neglected. Meanwhile, a non-destructive experimen-tal method is introduced to provide 300-nm longitudinal anddiffraction limited lateral resolution measurements of the in-bulkvolumetric intensity distribution. Our numerical methods arebenchmarked with experiments relying on propagation imagingunder various conditions (focusing optics, wavelength, spatialshape, propagation medium).
2. MODEL DESCRIPTIONS
A. Three-dimensional Fourier Transform (3D-FT) representa-tion of the field vectors near the focus
Using the form developed by Richards and Wolf [33], the time-dependent electric and magnetic fields ( E and H ) in the imageregime of a system can be expressed by Eq. (1). Here e and h are the time-independent electric and magnetic vectors, ω is theangular frequency. E ( x , y , z , t ) = (cid:60){ e ( x , y , z ) e − i ω t } , H ( x , y , z , t ) = (cid:60){ h ( x , y , z ) e − i ω t } . (1) At any point P ( x , y , z ) in the image space, the electric and mag-netic vectors e and h can be expressed in the form as a summa-tion of the plane waves that are leaving the aperture: e ( x , y , z ) = − ik π (cid:90)(cid:90) Ω a ( s x , s y ) s z e ik [ Φ ( s x , s y )+ s x x + s y y + s z z ] ds x ds y , h ( x , y , z ) = − ik π (cid:90)(cid:90) Ω b ( s x , s y ) s z e ik [ Φ ( s x , s y )+ s x x + s y y + s z z ] ds x ds y (2) where Φ ( s x , s y ) is the aberration function describing the opticalpath difference between the aberrated and the spherical wave-front along s . Here s is a unit vector pointing from a point in theexit aperture to the focus, a and b are the electric and magneticstrength vectors of the unperturbed electric and magnetic fieldsin the exit aperture, k is the wave number, and Ω is the solidangle formed by all the geometrical optical ray. The phase factorshown in Eq. (2) contains two parts, one is the scalar product ofvector s and vector r p , another is the vectorial aberration func-tion. In this section henceforth we only discuss the electric fieldsince, apart from the strength vector, the two equations in Eq. (2)are identical.Now let us consider a laser in-bulk focusing scenario. Asshown in Fig. 1, after the focusing element, this configurationconsists of materials 1 and 2 with refractive indices n and n ,respectively.In material 1 and at the interface ( z = − d ), the electric field isgiven by e ( x , y , − d ) = − ik π (cid:90)(cid:90) Ω C ( (cid:126) s ) × exp [ ik ( s x x + s y y − s z d )] ds x ds y (3) zy x φ θ O E (0) (r, θ) f φ d ’ c n n s r s a ( s ) /s W (e) ( S ) e u d Fig. 1.
Diagram showing laser focused by a lens into two me-dia separated by a planar interface.where C ( (cid:126) s ) = a ( s x , s y ) e ik Φ ( s x , s y ) s z (4) by combining the strength vector a ( s x , s y ) and the aberrationphase factor e ik Φ ( (cid:126) s ) into the new complex strength vector C ( (cid:126) s ) .The remaining phase factor in Eq. (3) thus contains only thescalar product of vector s and vector r p .Since there is no optical coating for all the cubical materialspresented in this paper, we assume that each plane wave com-ponents refraction at the interface obeys the Fresnel equations.To determine the transmitted field in the second material, wealso assume that the field in the second material is constructedby the superposition of refracted plane waves. As the complexstrength vector of the plane wave upon the interface is describedas C ( (cid:126) s ) , the strength vector of the transmitted plane wave can bedescribed as a linear function of C ( (cid:126) s ) , i.e, T · C ( (cid:126) s ) , where T is arefraction operator which is a function of angle of incidence and n , n . Therefore, the transmitted field in the second materialcan be written as e ( x , y , − d ) = − ik π (cid:90)(cid:90) Ω T · C ( (cid:126) S ) × exp [ ik ( s x x + s y y − s z d )] ds x ds y (5) On the other hand, as Török et al. [34] suggested, we can alsorepresent the field in the second material again as superpositionof plane waves, which is a solution of time-dependent waveequation and can be written as e ( r p ) = − ik π (cid:90)(cid:90) Ω F ( (cid:126) s ) exp ( ik (cid:126) s · r p ) ds x ds y . (6) One can notice that Eq. (5) is the boundary condition of Eq. (6).Let us now establish the relation between (cid:126) s and (cid:126) s .According to the law of refraction, k ( (cid:126) u × (cid:126) s ) = k ( (cid:126) u × (cid:126) s ) , (7) where (cid:126) u is the unit vector that is normal to the interface. Whena planar interface is presented (cid:126) u = (
0, 0, 1 ) , and we have k s x = k s x , k s y = k s y . (8) By taking the coordinate transformation, Eq. (6) yields e ( r p ) = − ik π (cid:90)(cid:90) Ω F ( (cid:126) s ) × exp ( ik (cid:126) s · (cid:126) r p ) J ( s x , s y ; s x , s y ) ds x ds y , (9) where J is the Jacobian of the coordinate transformation ob-tained from Eq. (8): J = (cid:18) k k (cid:19) , (10) As Eq. (9) must satisfy the boundary condition represented byEq. (5), we have F ( (cid:126) s , (cid:126) s ) = ( k k ) T · C ( (cid:126) s ) exp [ id ( k s z − k s z )] . (11) By substituting Eq. (11) into Eq. (9) we obtain the electric fieldin the second material: e ( x , y , z ) = − ik π k (cid:90)(cid:90) Ω T · C ( (cid:126) s ) × exp [ id ( k s z − k s z )] exp ( ik s z z ) × exp [ ik ( s x x + s y y )] ds x ds y . (12) The first phase factor exp [ id ( k s z − k s z )] stands for the aber-ration induced by the interface. The second phase factorexp ( ik s z z ) accounts for the phase accumulation when propa-gating along the z-axis, and the third term exp [ ik ( s x x + s y y )] represents the phase difference of the wave front at off-axispoints (x, y, z) with respect to the on-axis point (0, 0, z). Depend-ing on the chosen coordinates, the following forms of the wavevectors are equivalent: (cid:126) k = k x k y k z = k − s x − s y s z = k − sin φ cos θ − sin φ sin θ cos φ , (cid:126) k = k x k y k z = k − s x − s y s z = k − sin φ sin θ − sin φ cos θ cos φ . (13) Therefore, ds x ds y can be written as dk x dk y / k and theinterface-induced aberration can be written as Ψ ( φ , φ , d ) = d ( n cos φ − n cos φ ) . (14) The spherical-polar form of the complex strength vector afterthe interface is C ( φ , φ , r , θ ) = T ( φ , φ , θ ) · a ( r , θ ) e ik Φ ( r , θ ) / cos φ . (15) By using c ( φ , φ , r , θ ) to brief note T ( φ , φ , θ ) · a ( r , θ ) e ik Φ ( r , θ ) ,Eq. (12) can be finally rewritten as e ( x , y , z , d ) = − ik π k (cid:90)(cid:90) r < R [ c ( φ , φ , θ ) e ik Ψ ( φ , φ , d ) e ik z z / cos φ ] × exp [ − i ( k x x + k y y )] dk x dk y . (16) By using the method developed by Leutenegger et al. [31],we set | c | = r > R , the Debye-Wolf intergral is nowexpressed as the Fourier transform of the field distribution afterthe interface formed by the two material, ultimately resulting in e ( x , y , z , d ) = − ik π k F [ c ( φ , φ , θ ) e ik Ψ ( φ , φ , d ) e ik z z / cos φ ] . (17) As inspired by Leutenegger et al. [31], we used the chirpedZ-transform (CZT) algorithm [35] for the Fourier transforma-tion. This algorithm (i) allows breaking the relationship betweenthe sampling points (M) over the aperture radius and the min-imal sampling points (N) for fast Fourier transform (FFT), (ii) enables an implicit frequency offset, and (iii) internalizes thezero padding. Applying this generalization, one can adapt thesampling step in the focus field independently of the samplingstep in the input field, introduce an additional shift of the regionof interest, and finally improve the computational efficiency.
B. Representation of the complex strength vector
To determine the complex strength vector c , let us assumethat the incident field is linearly polarized. By choosing thecorresponding Cartesian coordinate and letting the y- and z-components of the incident electric field as zero, the incidentelectric strength vector can be written as E ( ) = E . (18) According to Török et al. [34], the transform relation be-tweeen the incident vector and the refracted vector after theinterface can be expressed by a refraction operator T . This oper-ator has a matrix form of T = A ( φ ) a a a a a a − a a a . (19) with a = τ p cos θ cos φ + τ s sin θ , a = ( τ p cos φ − τ s ) cos θ sin θ , a = τ p cos θ sin φ , a = τ s cos θ + τ p cos φ sin θ , a = τ p sin θ sin φ , a = τ p cos φ ,where τ s and τ p are the Fresnel transmission coefficients for s-polarized and p-polarized light, respectively. The function A ( φ ) is an apodization function that depends on the lens. Moreoverwhen the system obeys Abbe’s sine condition, i.e, is aplanatic,then A ( φ ) = f l (cid:112) cos φ (20) where f is the focal length of the lens in vacuum and l is anamplitude factor. It is assumed that the Abbe’s sine condition isverified in this step, as it is usually fullfiled in a corrected micro-scopic lens or in a single spherical lens with the stop located atthe lens.To represent the aberration phase factor and to scale it inwavelength λ , the aberration factor e ik Φ ( r , θ ) is rewritten as e i π W ( r , θ ) . To describe a circular lens-induced wavefront aberra-tion and calculate the deviation of the wavefront from an idealspherical shape, one widely used method relies on the Zernikepolynomials.By using the vendors provided lens data, the aberration func-tion can be calculated in the form of superposition of Zernikepolynomials through any homemade ray tracing program orcommercial software such as Zemax or Code V. In this paper, weused Zemax [36] to calculate the aberration function and usedthe notation convention defined by Noll [37] (so-called Zernikestandard polynomials): W = ∑ i = c i Z i ( ρ , θ ) , (21) Laser (pulse) (a)
ND Sample Focus lens MONA 0.85 Camera (InGaAs)NDMotorizedstage Tube lens (b)
Fig. 2.
Illustration of in-bulk propagation imaging process. (a)Experimental setup, MO: microscope objective, ND: neutraldensity filters. (b) Intensity distribution near the focal regionreconstructed from the stack of recorded images.where c i is the orthonormal Zernike coefficient computed byZemax, Z i is the corresponding Zernike standard polynomialand dimensionless radius ρ = r / r max normalized to the radiusof the entrance pupil.Finally, from Eq. (15), Eq. (18) and Eq. (21), the complexstrength vector after the interface can be written as: c = f l (cid:112) cos φ × a a a a a a − a a a e i π W e i π W
00 0 e i π W E . (22) By substituting Eq. (22) into Eq. (12), a CZT algorithm basedhighly-efficient propagation model is complete, it allows fastin-focus fields calculation by taking into account both the lens-induced and interface-induced aberrations.
3. EXPERIMENTAL ARRANGEMENT
A. Propagation imaging
Propagation imaging methods have been widely used to in-vestigate the linear or nonlinear optical effects in a medium[19, 24, 38, 39]. In this paper, to have a systematic compari-son with the numerical results, we introduce non-destructivemeasurements on the in-bulk volumetric intensity distributionswhich rely on (i) the focusing of the laser beam at the exit surfaceof the sample with desired thickness, and (ii) an inverted mi-croscope working in transmission for imaging the beam profilein the xy plane for various positions of the focusing objectivealong the propagation direction z . The experimental set-up isschematically depicted in Fig. 2.The measurements can be explicitly divided into two steps.The first step consista in focusing the laser beam with identicalcharacteristics (beam size, spectrum, phase distribution, etc.) asin the simulations. Here, we have chosen a femtosecond pulselaser to avoid the parasitic interference caused by residual mul-tiple reflections. In this ultrashort pulse regime, before focusing,the pulse energy is kept low (typically, a few tens of picojoules)to avoid nonlinear propagation effects. The second step consistsin imaging the exit surface of the sample with the inverted micro-scope working in transmission along the z-axis. This microscopeis composed of an infinity-corrected objective lens whose NAis larger than that of the focusing elements, a tube lens and an image sensing array. In this paper, the NA of the imaging ob-jective lens is 0.85 and its focal plane is precisely adjusted atthe exit surface of the sample under white light illuminationwith a translation stage. When there is no sample (focus in air),this plane can be arbitrarily chosen. Thanks to a precision stage(Physik Instrumente, M-126DC1), the laser intensity evolutionalong the z-axis is recorded by alternating 100-nm movementsof the focus lens (corresponding to n × n ), and image acquisitionsby the camera. The stack of images is then post-processed forreconstructing the fluence distribution as follows. Firstly, due tothe jitter of the laser, some recorded images may show a muchhigher maximum intensity compare to both the preceding andthe following one. These rare outlier images are replaced by theaverage of the preceding and the following images. Then, in thecorner of each image, the noise is calculated as the average ofthe pixel amplitude, and subsequently subtracted from all pixels.The image containing the maximum gray level is used for thenormalization of the whole stack. Two images before and afterthis latter image are defined for evaluating and correcting theresidual tilt of the collecting objective lens with respect to theincoming beam. Finally, the beam propagation is reconstructedby displaying a cross-section (along x or y) of each image at thecenter of the beam.To demonstrate the universality of our method, propagationimaging experiments have been carried out with four differentfocusing elements and two laser platforms. The first platformis based on an erbium-doped fiber laser (Raydiance Inc, SmartLight 50) with a wavelength λ of 1555 nm and a pulse durationof 860 fs. The second one is based on a prototype TRUMPFTruMicro 2030 Femto Edition laser with a wavelength λ of 1030nm and a pulse duration of 265 fs. The average power stabilityof both lasers is < 1%. B. Materials
Two singlet lenses (Thorlabs, LA1951-C and C240TME-C) and a50 × microscope objective lens (Mitutoyo, Apo NIR) are testedwith the Raydiance laser platform. The radius of the inputGaussian beam profile is 5.2-mm (at 1/ e ). A 20 × microscopeobjective (Mitutoyo, Apo NIR) is tested with the TruMicro 2030platform, the input beam profiles are shaped through ampli-tude modulation (slit) or phase modulation (phase plate). Torepresent the residual aberrations of the two singlet lenses, thenonzero terms of the Zernike standard coefficients are calculatedaccording to our experimental conditions and listed in the Ap-pendix A (Table A1). The laser beams are focused in air ( n =1)and in crystalline silicon (c-Si, n =3.475 at λ =1555 nm). Detailedinformation on the materials used in the experiment are listedin Table 1.
4. RESULTS AND DISCUSSIONS
In this section we systematically compare our numerical resultswith the experiments. Under the in-air focus condition, theexperimental results benchmarked our numerical model. Thefurther numerical investigations of the in-silicon focus condition,in return, pointed out a limitation of the experimental method.
A. Gaussian beam focused by an aspheric lens in air
To achieve in-bulk processing, tightly focused laser beams areoften required. A cost-efficient solution to get the diffraction-limited high quality tight focusing is to use an aspheric lens. Inthis paper, an aspheric lens (Thorlabs C240TME-C) with NA =
Vendor, Lens f [mm] NA Φ [mm] 1/ e radius [mm] λ [nm] Focused-in medium Presened inThorlabs, C240TME-C 8.0 0.50 8.0 5.2 1555 Air Fig. 3Thorlabs, LA1951-C 25.3 - 25.4 5.2 1555 Air Fig. 4Mitutoyo, 20 × Plan Apo NIR 10.0 0.40 8.0 User-defined 1030 Air Fig. 5Mitutoyo, 50 × Plan Apo NIR HR 4.0 0.65 5.2 5.2 1555 Si Fig. 6
Table 1.
Detailed information of the lenses used in this paper. f is the focal length, NA is the numerical aperture, Φ is the clearaperture diameter, λ is the center wavelength of the incident beam.0.5 is used for the first demonstration. As the 1/ e radius ofthe input beam is 5.2-mm, the lens is overfilled by a factor of1.3. Therefore, to calculate the aberrations induced by this lens,the apodization factor G (refers to the rate of decrease of thebeam amplitude as a function of radial pupil coordinate) is setas √ A ( ρ ) = e − G ρ , where ρ is thenormalized pupil coordinate. Under this approximation, theaberrations induced by the lens are calculated and representedby the standard Zernike coefficients. The nonzero coefficientsup to the 37th term are listed in Appendix A, Table A1. As afirst demonstration, the focused-in medium is air and the corre-sponding refractive index is 1. Based on Eq. (17) and Eq. (22),the normalized intensity distribution near the focus is calculatedand presented in Fig. 3(c-d). By applying the method presentedin section 2(A), the corresponding experimental results are ob-tained and presented in Fig. 3(a-b). The normalized longitudinalintensity distributions are shown in Fig. 3(a) and (c), and theintensity distributions at the focal plane are shown in Fig. 3(b)and (d). The overall distributions are very similar. To quanti-tatively evaluate how good the simulation results fit with theexperimental one, the normalized experimental intensity profilesare compared to the simulated ones in Fig. 3(e-f), together withthe absolute values for the intensity difference. Based on thesequantitative analyses, the root-mean-square deviation (RMSD)of those two normalized intensities are calculated for each axis.As shown in Fig. 3(e), by setting the intensity maximized posi-tion as the origin, the intensity profiles and their differences areplot from -5 µ m to 5 µ m. The experimental beam width (at 1/ e )is measured as 2.64 µ m, and the simulated one is 2.93 µ m. In thistransverse direction, the aforementioned RMSD is calculatedas 0.0305. We applied the same methods to the results alongthe longitudinal direction. The normalized intensity profilesare compared in Fig. 3(f) along the z-axis. The RMSD of thelongitudinal profiles is calculated as 0.0322. To anchor a refer-ence, we also simulated the intensity distribution near the focuswithout considering the influence of the aberrations (not shownhere). In this case, the RMSD of the normalized profiles alongx, z-axis are 0.0382 and 0.0365. In other words, by taking theaberrations induced by the C240TME-C lens into account, theRMSDs have decreased by 20.2% and 11.8% along the x- andz-axis, respectively. B. Gaussian beam focused by a plano-convex lens in air
The plano-convex lens is one of the simplest converging lensesthat has been widely used to focus collimated light. The spheri-cal aberrations can be minimized due to the asymmetric designby placing the curved surface face toward the collimated beam,however, it cannot be completely eliminated. In this demonstra- (e)
50 µm5 µm 5 µm (a) (b)(c) (d)(f) kk exp. exp.sim. sim. kk I [arb. unit]I [arb. unit]
Fig. 3.
Intensity distribution near the focus of an asphericallens in air with a Gaussian beam as the input. (a) Experimen-tally measured intensity along xz plane; (b) experimentallymeasured intensity at the focal plane; (c) simulated intensityalong xz plane; (d) simulated intensity at the focal plane. (e)Comparison of the intensity profiles at the focal plane alongthe x-axis, (f) along the z-axis. tion, we choose a plano-convex lens (Thorlabs LA1951-C) with f = mm and used the two methods mentioned in Section2 and 4 to evaluate the intensity distribution at its focus. Thesame collimated laser beam that was used in the previous sec-tion is focused by this lens in air. To calculate the lens-inducedaberrations, the beam amplitude is normalized to unity at thecenter of the pupil, at other points of the entrance pupil theamplitude is given by A ( ρ ) = e − G ρ , where G = 4 and ρ is thenormalized pupil coordinate. Under this approximation, theaberrations induced by the lens are calculated and representedby the standard Zernike coefficients. In Appendix A, Table A1the nonzero coefficients are listed up to the 37th term. Finally,the normalized intensity distribution near the focus is calculatedbased on Eq. (17) and Eq. (22) and presented in Fig. 4(c-d). Thecorresponding experimental results are presented in Fig. 4(a-b).Fig. 4(a) and (c) show the normalized longitudinal intensitydistributions and Fig. 4(b) and (d) shown the intensity distribu-tion at the focal plane. To quantitatively evaluate the agreement,in Fig. 4 (e-f), we plot the normalized intensity value along the xand y central segment at the focal plane.As shown in Fig. 4 (e), along the x-axis, the experimentalbeam diameter at 1/ e is measured as 10.6 µ m and the simulatedone is 11.4 µ m. The RMSD of the two profiles is 0.051. Similarcomparisons are applied to the profiles along the y-axis. Theexperimental beam width is measured as 11.3 µ m, the simulatedone is 11.4 µ m and the RMSD is 0.046. From the quantitativecomparison one can notice that, in contrast to the simulationresults (d), the experimental ones (b) exhibit an elliptical feature.This deviation is due to the fact that for the simulations theincident beam has been oversimplified as a circular Gaussianbeam. Apart from that, the simulated results fit well with theexperiment.It is also worth noting that, given that the experimental resultsshown in Fig. 4 (a) are composed by 6000 individual frames, onthe left-hand side of this figure some frames are misaligneddue to the environmental disturbance such as air flows andvibrations. Except for that, one can notice that the aberrationfeatures as well as the overall lengths of the focal regime areidentical. C. Exotic beams focused by an objective lens in air
Microscope objective lenses are also widely used for laser beamfocusing and material processing. With a proper design, theaberrations can be well corrected for the design wavelengthrange. In this section, we choose a 20 × Mitutoyo Plan Apoobjective lens (NA = 0.4) as the focus lens. The center wavelengthof the laser source is 1030-nm. Instead of using a standardGaussian beam as the input, we choose amplitude- and phase-shaped beams for the investigation. To precisely simulate thefocusing conditions, we first record the intensity profile of theexotic beams at the entrance pupil of the objectives with an arrayof image sensors, and then calculate the amplitude distributionby taking the square root of the intensity profile. Based onthose measured “user defined” amplitude profiles and theirpolarization states, we obtain the complex fields on the entrancepupil. Finally, based on Eq. (17) and Eq. (22), the fields near thefocus can be calculated. These calculated intensity distributionsare compared to the experimental results.A circular cross-sectional focusing is often required fortransversal writing of single-mode waveguides or microfluidicchannels. Slit beam shaping is a simple technique that providessuch an isotropic resolution in transverse and vertical direction kk I [arb. unit](c) (d) 00.20.40.60.81.0 exp. exp.sim. kk sim. (e) Fig. 4.
Intensity distribution near the focus of a plano-convexlens in air with a Gaussian beam as the input. (a) Experimen-tally measured intensity along xz plane; (b) experimentallymeasured intensity at the focal plane; (c) simulated intensityalong xz plane; (d) simulated intensity at the focal plane. Com-parisons of the intensity profiles at the focal plane (e) along thex-axis; (f) along the y-axis.[40, 41]. In this section, we use the slit beam shaping as a firstexample to prove that our methods are also applicable for theamplitude shaped beam.As shown in Fig. 5 (a), a linearly-polarized Gaussian beamwith a 1/ e radius of 5.2 mm is cut by a 1.38-mm width slit. Theexperimentally measured intensity distributions obtained byfocusing this beam are shown in Fig. 5(b-d), and the simulatedresults are shown in Fig. 5(e-g). At the focal plane, as shownin Fig. 5 (d), the beam width along the x-axis is measured as2.9 µ m and along the y-axis is 10.9 µ m. These beam width valuesat the focal plane are identical to those obtained in the simula-tions in Fig. 5(g). Similarly, an excellent agreement between theexperiments and the calculations is found in the xy and the yz planes, as shown in (b-c) and (e-f), respectively. As expectedwith such slit beam shaping, the light is focused tighter in the xz plane than in the yz plane. This mainly originates from the lossin the effective numerical aperture along the y-axis due to thefact that the input beam does not overfill the entrance pupil asfor the x-axis. All in all, a near-perfect agreement between theexperimental measurements and the calculations demonstratesthat our proposed model is a powerful tool even for investi-gating laser-matter interaction scenarios where sophisticatedanisotropic beams are employed.However, one should be careful when applying our compu-tational method under extremely asymmetric conditions such asthe line-focus microscopy (LTM), as Wolf and Li have indicated (b) (c) (d)(e) (f) (g)(i) (j) (k)(l) (m) (n)(h) 100µm 15 µm15 µm5 mm5 mm 10 µm100µm 10 µm(a) xy exp.xz exp.yz exp.xysim.xz sim.yz sim.xyexp.xz exp.yz exp.xysim.xz sim.yz sim.xy kkkk kkkk Fig. 5.
Intensity distribution of the input beam at the entrancepupil for (a) a Gaussian beam cut by a slit with a width of1.38-mm, (h) a Gaussian beam accumulating a spiral phase(bottom). Intensity distributions near the focus: (b-d) (i-k)experimentally measured along xz , yz and xy plane; (e-g) (l-n)simulated along xz , yz and xy plane.that their Debye integral representation should be consideredonly beyond a critical value of the Fresnel number [42]. TheFresnel number is a dimensionless number, N = a / λ f , whichreflects the relative contribution of focusing against diffractioneffects for a given aperture radius a , focal length f , and wave-length λ . For more rigorous approaches dedicated to this specificproblem, one can refer to the works of De la Cruz et al. [43] andLou et al. [44]. They have illustrated that, for Fresnel numbersclose to unity, the focus shifts backward, thus leading to astig-matic focusing when the circular symmetry of the input light isbroken. In our case where N ≈
41, such a backward shift couldnot be observed experimentally.In the same way as we have studied the influence of am-plitude shaping, the impact of phase shaping has also beeninvestigated. Over the last decade, the helical wavefront is oneof the most extensively studied complex phase shapes of light.This type of light beams have an azimuthal phase dependence ofexp ( il θ ) , where l is the topological charge and θ is the azimuthangle. The optical vortex has many innovative applications inoptical tweezers [45], atom manipulation [46] and material pro-cessing [47]. When focused, this optical vortex forms a ringinstead of a spot in the focal plane. In this section, experimen-tally, we used a spiral phase plate (SPP) to discretely generatean azimuthal phase distribution of exp ( − i θ ) . The topologicalcharge is -1, and the number of discrete steps is 12. After thephase plate, the laser beam is focused by the same 20 × objectivelens as previously. In the simulation, to obtain the correspondingcomplex field as the input, we multiplied the measured ampli-tude distribution with a discrete spiral phase map that exhibitsthe same phase distribution as the one used in the experiment.The measured input intensity distribution and the calculatedspiral phase map are presented in Fig. 5 (h).The simulations in Fig. 5(l-n) exhibit similar optical vortexfeatures as the experimental measurements in Fig. 5(i-k). From the longitudinal intensity distributions (i-j) and (l-m) one cansee that in both cases the light waves along the propagation axiscancel each other out and the ratios of the outer and inner radiiare identical. At the focal plane, as shown in Fig. 5 (k) and (n),this ratio is measured as 4.3 (experimental) and 4.0 (simulated).While, overall, the experimental and theoretical results are ingood agreement, one can still note some deviations in the spotelliptically and the intensity homogeneity, most likely provokedby imperfections of the spiral phase plate as well as residualmisalignments of this plate. D. Gaussian beam focused by an objective lens in silicon
In the previous sections, the laser beams are focused in air. Inthis section, the in-bulk focus condition is chosen. Unlike inthe aforementioned conditions, large experimental to simulationdeviation are observed in this scenario. These deviations can beascribed to the experimental procedure to acquire the intensityprofile.When a planar interface is present, according to the aberra-tion function [Eq. (14)] as well as noted in experimental observa-tions [48], the interface-induced aberration is a function of thefocus depth, the NA and the refractive index of the bulk material.By applying the same methodology as described above, we canachieve the experimental and simulated results of the longitu-dinal intensity distributions near the focus. When a Gaussianbeam ( λ : 1555-nm, 1/ e radius: 5.2-mm) is focused into a 5-mmthick crystalline silicon (c-Si) sample ( n = 3.475 at 1555-nm) byan objective lens with NA of 0.65, the experimental result is pre-sented in Fig. 6(a) and the corresponding simulation is shown inFig. 6(b). Contrary to the previous cases, this comparison showssignificant deviations, especially in terms of the overall length.The main reason for this deviation is that, given the employedpropagation imaging measurements, the experimental resultsdo not correspond to the intensity distribution at 5-mm focaldepth, but it is the stacking of intensity profiles at the exit surfaceof the 5-mm thick sample for a focusing objective lens movedstep-by-step from -1.5 to +0.5 mm (the position 0 correspondingto the exit surface).One could wonder why not directly measuring the inten-sity distribution in the bulk of the material. Unfortunately, itis nearly impossible to measure the real intensity distributionexperimentally. This problem cannot be solved by fixing theposition of the focusing objective lens and moving the recordingmicroscope since, in this case, the recorded images would bestrongly affected by the aberrations provoked by the refractiveindex mismatch at the back surface of the sample. The only wayto avoid the artificial result is to compensate the depth relatedaberration after each movement. Considering 200 movementsduring one measurement, it is unrealistic to compensate it withany simple means such as a correction collar. However, theexperimental result shown in Fig. 6(a) can still be exploited forbenchmarking our model.In order to carry out simulations in a situation that is com-parable to the experiments, we first calculate the correspondingintensity profile at the exit surface of the 5-mm sample for differ-ent focal depths (from 3.5 to 5.5 mm). The focal depth incrementstep size is 5 µ m. After the calculation, we stack the intensityprofiles from the left to the right with the focus depth decrementfrom 5.5 mm to 3.5 mm. With these two steps, the experimentalimage acquisition procedure and the corresponding intensitydistribution are simulated and displayed in Fig. 6(c). While thecomparison between the experiments (a) and the simulation for afixed focus at 5-mm focal depth (b) is mediocre, the experiments
500 µm00.20.40.60.81.0 k
15 µm exp.sim.simulated experiments d=5mmd=3.5-5.5mm (a)(b) kk (c) Fig. 6.
Intensity distributions near the focus, (a) Experimen-tally measured intensity along xz plane; (b) simulated inten-sity along xz plane with focal depth of 5-mm; (c) simulatedintensity at the exit surface of the 5-mm sample with focaldepth varying from 3.5-mm to 5.5-mm.and the simulated experiments (c) are in excellent agreement.One should emphasize that the simulation at 5-mm focaldepth (c) corresponds to the real intensity distribution, whichcould not be measured experimentally in a simple way. If oneaims at measuring the intensity distribution experimentally, ex-tra attention must be paid to the potentially unrealistic characterof the results.To quantitatively illustrate the deviation of the experimentalmeasured results and the actual in-bulk intensity distribution,we compare in Fig. 7 the actual simulated intensity distributionfor a fixed focal depth to the simulated experimental results fordifferent numerical apertures and silicon thicknesses.Given that this deviation is caused by the aberrations inducedby the planar interface of the bulk material, the focusing condi-tions (NA and focal depth) will determine how significant thisdeviation is. In order to illustrate this relationship and high-light the limitations of propagation imaging experiments, theRMSDs under different numerical apertures and sample thick-nesses are calculated in Fig. 7(a). As shown in the 3D bar chart,this deviation becomes negligible for small numerical aperturesand sample thicknesses. In Fig. 7(b-d) comparisons of the cross-section profiles along the z-axis are illustrated for three selectedcases. When the NA is decreased to 0.26 and the focal depthis decreased to 1 mm, the experimentally measured results arealmost identical to the actual in-bulk intensity distribution. TheRMSD for this case is 0.009. This value is even smaller than thedeviation of the experimental measured and simulated inten-sity distributions in air, for which there is no influence of theinterface-induced aberration. Therefore, in this case, it is safeto use the propagation imaging measurements for representingthe actual intensity distribution represent the actual intensitydistribution and further apply this propagation imaging method z [µm]Intensity [a.u.] (a)(b)(c)(d) bcd SimulationSimulated experimentSimulationSimulated experimentSimulationSimulated experimentNA = 0.65d = 5mmNA = 0.40d = 2.5mmNA = 0.26d = 1.0mm
Fig. 7. (a) The RMSDs of the simulated actual intensity distri-bution and the simulated experimental result under differentNAs and focal depths. Comparison of the cross-section inten-sity profiles near the focus along the z-axis with (b) NA = 0.65and focal depth of 5 mm; (c) NA = 0.4 and focal depth of 2.5mm; (d) NA = 0.26 and focal depth of 1 mm.for investigating more complex nonlinear propagation problems[19].
5. CONCLUSION
In this paper, we introduced a benchmarked model and an ex-perimental tool for analyzing the field distribution when a laserbeam is focused in air or in bulk materials. The vectorial anal-ysis model considered the lens-induced as well as the planarinterface-induced aberrations. The in-bulk propagation imagingsetup provides a 300-nm longitudinal resolution and diffractionlimited lateral resolution. Using the tools introduced in thispaper, one can deal with a wide variety of focus conditions inwhich arbitrary input fields, non-aplanatic lenses, and in-bulkfocus might be involved.While the numerical simulation tool is applicable, we havealso pointed out that for high numerical apertures and signifi-cant focal depths, the experimental results acquired by propa-gation imaging methods that are similar to the one describedin this paper might lead to large deviations against the actualintensity distribution in the desired focal depth. One should calculate the RMSD or at least check the RMSD chart beforeapplying these methods to any linear or nonlinear propagationimaging experiments.We plan to continue the development of InFocus [1] in thespirit of open-source and would be pleased to find collabora-tors. We anticipate that our proposed model and correspondingsoftware InFocus can be utilized in countless laser processing ap-plications involving various wavelengths, beam shapes, phasedistributions and focusing optics.
APPENDIX A: ZERNIKE STANDARD COEFFICIENTS OFTHE SINGLE LENS
This appendix lists the nonzero terms of the Zernike standardcoefficients calculated at the exit pupil of the single lens LA1951-C and C240TME-C. For both lenses, the properties of the incidentlaser beams are the same, i.e., a wavelength of 1555 nm and a1/ e radius of 5.2 mm. The refractive index of the medium afterthe lens is 1.LA1951-C C240TME-C Polynomials Z Z √ ( r − ) Z √ ( r − r + ) Z √ ( r − r + r − ) z √ ( r − r + r − r + ) Table A1.
Nonzero terms of the Zernike standard coefficientscalculated at the exit pupil of a plano-convex lens ThorlabsLA1951-C and an aspherical lens C240TME-C.
Acknowledgment.
This research has been supported by the Bun-desministerium für Bildung und Forschung (BMBF) through the NU-CLEUS project, grant no. 03IHS107A, as well as the glass2met project,grant no. 13N15290.
Disclosures.
The authors declare no conflicts of interest.
REFERENCES
1. Q. Li, “InFocus (version 1.0),” https://github.com/QF06/InFocus. Ac-cessed: 2021-01-18.2. Y. Bellouard, A. Said, M. Dugan, and P. Bado, “Fabrication of high-aspect ratio, micro-fluidic channels and tunnels using femtosecondlaser pulses and chemical etching,” Opt. Express , 2120–2129(2004).3. R. Osellame, V. Maselli, R. M. Vazquez, R. Ramponi, and G. Cerullo,“Integration of optical waveguides and microfluidic channels both fabri-cated by femtosecond laser irradiation,” Appl. Phys. Lett. , 231118(2007).4. V. Maselli, J. R. Grenier, S. Ho, and P. R. Herman, “Femtosecond laserwritten optofluidic sensor: Bragg grating waveguide evanescent probingof microfluidic channel,” Opt. Express , 11719–11729 (2009).5. F. He, Y. Cheng, Z. Xu, Y. Liao, J. Xu, H. Sun, C. Wang, Z. Zhou,K. Sugioka, K. Midorikawa et al. , “Direct fabrication of homogeneousmicrofluidic channels embedded in fused silica using a femtosecondlaser,” Opt. Lett. , 282–284 (2010).6. K. M. Davis, K. Miura, N. Sugimoto, and K. Hirao, “Writing waveguidesin glass with a femtosecond laser,” Opt. letters , 1729–1731 (1996). 7. M. Chambonneau, Q. Li, M. Chanal, N. Sanner, and D. Grojo, “Writingwaveguides inside monolithic crystalline silicon with nanosecond laserpulses,” Opt. Lett. , 4875–4878 (2016).8. I. Pavlov, O. Tokel, S. Pavlova, V. Kadan, G. Makey, A. Turnali, Ö. Yavuz,and F. Ilday, “Femtosecond laser written waveguides deep inside sili-con,” Opt. Lett. , 3028–3031 (2017).9. W. Gebremichael, L. Canioni, Y. Petit, and I. Manek-Hönninger,“Double-track waveguides inside calcium fluoride crystals,” Crystals , 109 (2020).10. X. Wang, X. Yu, M. Berg, B. DePaola, H. Shi, P. Chen, L. Xue, X. Chang,and S. Lei, “Nanosecond laser writing of straight and curved waveg-uides in silicon with shaped beams,” J. Laser Appl. , 022002 (2020).11. M. Chambonneau, D. Richter, S. Nolte, and D. Grojo, “Inscribing diffrac-tion gratings in bulk silicon with nanosecond laser pulses,” Opt. Lett. , 6069–6072 (2018).12. J. Zhang, M. Geceviˇcius, M. Beresna, and P. G. Kazansky, “Seeminglyunlimited lifetime data storage in nanostructured glass,” Phys. Rev. Lett. , 033901 (2014).13. A. Crespi, R. Ramponi, R. Osellame, L. Sansoni, I. Bongioanni, F. Scia-rrino, G. Vallone, and P. Mataloni, “Integrated photonic quantum gatesfor polarization qubits,” Nat. Commun. , 1–6 (2011).14. K. Lammers, M. Ehrhardt, T. Malendevych, X. Xu, C. Vetter, A. Al-berucci, A. Szameit, and S. Nolte, “Embedded nanograting-basedwaveplates for polarization control in integrated photonic circuits,” Opt.Mater. Express , 2560–2572 (2019).15. S. Richter, F. Zimmermann, A. Tünnermann, and S. Nolte, “Laser weld-ing of glasses at high repetition rates–fundamentals and prospects,”Opt. & Laser Technol. , 59–66 (2016).16. G. Zhang, J. Bai, W. Zhao, K. Zhou, and G. Cheng, “Interface modifica-tion based ultrashort laser microwelding between sic and fused silica,”Opt. Express , 1702–1709 (2017).17. K. Cvecek, S. Dehmel, I. Miyamoto, and M. Schmidt, “A review on glasswelding by ultra-short laser pulses,” Int. J. Extrem. Manuf. , 042001(2019).18. E. Penilla, L. Devia-Cruz, A. Wieg, P. Martinez-Torres, N. Cuando-Espitia, P. Sellappan, Y. Kodera, G. Aguilar, and J. Garay, “Ultrafastlaser welding of ceramics,” Science , 803–808 (2019).19. M. Chambonneau, Q. Li, V. Y. Fedorov, M. Blothe, K. Schaarschmidt,M. Lorenz, S. Tzortzakis, and S. Nolte, “Taming ultrafast laser filamentsfor optimized semiconductor–metal welding,” Laser & Photonics Rev. p.2000433 (2020).20. R. Meyer, L. Froehly, R. Giust, J. Del Hoyo, L. Furfaro, C. Billet, andF. Courvoisier, “Extremely high-aspect-ratio ultrafast Bessel beam gen-eration and stealth dicing of multi-millimeter thick glass,” Appl. Phys.Lett. (2019).21. A. Couairon and A. Mysyrowicz, “Femtosecond filamentation in trans-parent media,” Phys. Reports , 47–189 (2007).22. L. Bergé, S. Skupin, R. Nuter, J. Kasparian, and J.-P. Wolf, “Ultra-short filaments of light in weakly ionized, optically transparent media,”Reports on Prog. Phys. , 1633 (2007).23. E. G. Gamaly, L. Rapp, V. Roppo, S. Juodkazis, and A. V. Rode,“Generation of high energy density by fs-laser-induced confined mi-croexplosion,” New J. Phys. , 025018 (2013).24. V. Y. Fedorov, M. Chanal, D. Grojo, and S. Tzortzakis, “Accessing ex-treme spatiotemporal localization of high-power laser radiation throughtransformation optics and scalar wave equations,” Phys. Rev. Lett. ,043902 (2016).25. P. K. Sahoo, T. Feng, and J. Qiao, “Dynamic pulse propagation mod-elling for predictive femtosecond-laser-microbonding of transparentmaterials,” Opt. Express , 31103–31118 (2020).26. J. Liu, B. Xu, and T. C. Chong, “Three-dimensional finite-differencetime-domain analysis of optical disk storage system,” Jpn. J. Appl. Phys. , 687 (2000).27. A. Couairon, O. G. Kosareva, N. Panov, D. Shipilo, V. Andreeva,V. Jukna, and F. Nesa, “Propagation equation for tight-focusing bya parabolic mirror,” Opt. Express , 31240–31252 (2015).28. D. E. Shipilo, I. A. Nikolaeva, V. Y. Fedorov, S. Tzortzakis, A. Couairon,N. A. Panov, and O. G. Kosareva, “Tight focusing of electromagnetic fields by large-aperture mirrors,” Phys. Rev. E , 033316 (2019).29. P. Varga and P. Török, “Focusing of electromagnetic waves byparaboloid mirrors. i. theory,” JOSA A , 2081–2089 (2000).30. P. Varga and P. Török, “Focusing of electromagnetic waves byparaboloid mirrors. ii. numerical results,” JOSA A , 2090–2095(2000).31. M. Leutenegger, R. Rao, R. A. Leitgeb, and T. Lasser, “Fast focus fieldcalculations,” Opt. Express , 11277–11291 (2006).32. J. Lin, O. G. Rodríguez-Herrera, F. Kenny, D. Lara, and J. C. Dainty,“Fast vectorial calculation of the volumetric focused field distribution byusing a three-dimensional Fourier transform,” Opt. Express , 1060(2012).33. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems,ii. structure of the image field in an aplanatic system,” Proc. Royal Soc.London. Ser. A. Math. Phys. Sci. , 358–379 (1959).34. P. Török, P. Varga, Z. Laczik, and G. Booker, “Electromagnetic diffrac-tion of light focused through a planar interface between materials ofmismatched refractive indices: an integral representation,” JOSA A ,325–332 (1995).35. J. L. Bakx, “Efficient computation of optical disk readout by use of thechirp z transform,” Appl. Opt. , 207 (1976).38. C. Pasquier, P. Blandin, R. Clady, N. Sanner, M. Sentis, O. Utéza, Y. Li et al. , “Handling beam propagation in air for nearly 10-fs laser damageexperiments,” Opt. Commun. , 230–238 (2015).39. A. Wang, A. Das, D. Grojo et al. , “Ultrafast laser writing deep inside sili-con with thz-repetition-rate trains of pulses,” Research , 8149764(2020).40. Y. Cheng, K. Sugioka, K. Midorikawa, M. Masuda, K. Toyoda,M. Kawachi, and K. Shihoyama, “Control of the cross-sectional shapeof a hollow microchannel embedded in photostructurable glass by useof a femtosecond laser,” Opt. Lett. , 55–57 (2003).41. M. Ams, G. Marshall, D. Spence, and M. Withford, “Slit beam shapingmethod for femtosecond laser direct-write fabrication of symmetricwaveguides in bulk glasses,” Opt. Express , 5676–5681 (2005).42. E. Wolf and Y. Li, “Conditions for the validity of the debye integralrepresentation of focused fields,” Opt. Commun. , 205–210 (1981).43. A. R. De la Cruz, A. Ferrer, J. Del Hoyo, J. Siegel, and J. Solis, “Mod-eling of astigmatic-elliptical beam shaping during fs-laser waveguidewriting including beam truncation and diffraction effects,” Appl. Phys. A , 687–693 (2011).44. K. Lou, S. Granick, and F. Amblard, “How to better focus waves byconsidering symmetry and information loss,” Proc. Natl. Acad. Sci. ,6554–6559 (2018).45. M. Padgett and R. Bowman, “Tweezers with a twist,” Nat. Photonics ,343–348 (2011).46. K. Ladavac and D. G. Grier, “Microoptomechanical pumps assembledand driven by holographic optical vortex arrays,” Opt. Express ,1144–1149 (2004).47. C. Hnatovsky, V. G. Shvedov, W. Krolikowski, and A. V. Rode, “Materialsprocessing with a tightly focused femtosecond laser vortex pulse,” Opt.Lett. , 3417–3419 (2010).48. Q. Li, M. Chambonneau, M. Chanal, and D. Grojo, “Quantitative-phasemicroscopy of nanosecond laser-induced micro-modifications insidesilicon,” Appl. Opt.55