Experimental estimation of the longitudinal component of a highly focused electromagnetic field

The detection of the longitudinal component of a highly focused electromagnetic beam is not a simple task. Although in recent years several methods have been reported in the literature, this measure is still not routinely performed. This paper describes a method that allows us to estimate and visualize the longitudinal component of the field in a relatively simple way. First, we measure the transverse components of the focused field in several planes normal to the optical axis. Then, we determine the complex amplitude of the two transverse field components: the phase is obtained using a phase recovery algorithm, while the phase difference between the two components is determined from the Stokes parameters. Finally, the longitudinal component is estimated using the Gauss’s theorem. Experimental results show an excellent agreement with theoretical predictions.


Scientific Reports
| (2021) 11:17992 | https://doi.org/10.1038/s41598-021-97164-z www.nature.com/scientificreports/ or in information optics 21 . In the present paper we propose a method to estimate the longitudinal components in the focal region of a highly focused beam based on the information contained on complex amplitudes of the corresponding transverse part. This paper is organized as follows: first, we describe the theoretical background for the longitudinal component estimation. Thereafter, in "Methods", we describe the adapted algorithm used to retrieve the phase of an electromagnetic field, and explain the experimental setup. The estimation of the longitudinal component can be found in "Results and discussion" and finally, our conclusions are summarized in the "Concluding remarks" section.

Longitudinal and transverse components of an electromagnetic field
The electromagnetic field in free space must satisfy Maxwell's Equations, specifically, the Gauss' theorem where E(r) is the electric field and r is the position vector. Time dependence is dropped since in this work we only consider quasi-monochromatic waves. For plane waves, the Gauss' theorem represents the transverse condition for the electromagnetic field: the polarization direction of the beam is perpendicular to the direction of propagation. However, non-homogeneous fields can be understood as being composed by a set of plane waves traveling in different directions. Therefore, the direction of propagation is not perfectly defined, and we cannot strictly talk about transverse waves.
Without loss of generality, we consider the propagation of electromagnetic waves with respect to a reference axis, say the z axis. We split the electromagnetic field E(r) in the following way where E ⊥ (r) and E z (r) are the transverse and parallel components to the z axis respectively, and e z is the unit vector in the direction of the z axis. Introducing Eq. (2) into Eq. (1) we obtain the following identity with ∇ ⊥ = e x ∂ ∂x + e y ∂ ∂y and e x , e y , e z is a orthogonal triad of unit vectors. If we consider each component of the electromagnetic field as composed by a superposition of plane waves 57 , where k ⊥ = (k x , k y ) and k z are the transverse and longitudinal wave-vectors, respectively, satisfying k 2 = k 2 ⊥ + k 2 z : In the present paper, we only consider the case k z real due to the experiments we are carrying out.

Methods
Phase recovery algorithm. To recover the longitudinal component of an electromagnetic field by means of Eq. (12), it is necessary to determine the complex amplitudes of the transverse components. However, only the irradiance of the field can be recorded without the use of specific techniques such as holography. However, thanks to the mathematical properties of the electromagnetic fields, a variety of methods to recover the phase by means of irradiance measurements have been developed 58,59 . Among them, iterative algorithms such as the Hybrid Input-Output 50 , might be suitable for the estimation of the longitudinal component.
In order to obtain a fair estimation of the phase of the beam, we record the modulus of the electromagnetic field at P different planes perpendicular to the z−axis z j (j = 1, . . . , P), with i = x, y . Note we should determine the phase associated to each polarization component independently and thus, the procedure should be repeated twice. Complex amplitudes U j = A j e iφ j are obtained by means of the following procedure: 1. Assign an initial phase estimation for the first modulus, U 1 = A 1 e iφ 1 . 2. Propagate complex amplitude U 1 to plane z 2 using a Fresnel propagation. The modulus is discarded, and the phase is assigned to the experimental measure of the irradiance. 3. Repeat the previous step until we reach the final plane, P 4. Back-propagate U P to the first plane. This gives us the next estimation for the initial phase e iφ 1 . The error between the recovered modulus and the experimental one is measured at this step. This process is repeated until the error measure arrives at the prescribed value and/or the measure stagnates.
This kind of iterative algorithms have two main drawbacks: slow rate of convergence and stagnation at local minima of the error function 59,60 . We address each of these problems separately in the following subsections.
Propagation method and local minima evasion. The propagation of the electromagnetic field is performed using the angular spectrum of plane waves and the free space transfer function (see Eqs. (4) and (10)) 55,57 . Note that other propagation methods could be used as long as the size of the window is not modified 61,62 . The relationship between two planes separated a distance z is where Û (k x , k y ; 0) is the spectrum at the first plane.
Recorded moduli might contain a certain amount of noise. This means that the calculated spectra from the experimental recordings have non-zero high frequency components. Moreover, most of the recorded amplitudes have values near to the lower end of the dynamic range of the sensor. These two effects combined can produce the iterative algorithm to prioritize the fitting of noise instead of the actual beam values.
To avoid this effect, we limit the extent of the field spectra. The limiting frequency is estimated by considering the Fourier transforms of the recorded irradiance A 2 j distributions. Because the maximum frequency extent of the irradiance is twice the cutoff frequency of the complex amplitude, in our calculation we set the radius of the support region by considering the frequency for which the Fourier Transform of the irradiance falls to a value close to zero.
Acceleration of the convergence speed. Iterative algorithms, although robust, display slow rates of convergence. To overcome this limitation, we use the ad hoc acceleration procedure developed by Biggs and Andrews 63 . This algorithm is independent of the exact shape of the phase recovery method and can be used with any iterative method.
The acceleration algorithm is implemented as follows. We define the following parameters www.nature.com/scientificreports/ where x k is the estimation for the initial phase at the kth iteration, h k is the difference between the current phase and the previous one, y k is the estimation for the next value of the phase at the kth iteration, and α k is the kth iteration acceleration parameter, explained below. The next value of the phase x k+1 is determined using g k , a parameter defined as the difference between ψ(y k ) and y k ; ψ(·) represents a single iteration of the phase recovering algorithm. Note we work with complex exponentials to ensure a robust convergence: the phase has a branching point at 2π which would create meaningless gradients or oscillations around this point.
The definition of the acceleration parameter is 63 where the summation is taken over all values of the elements of g k and * denotes complex conjugate. We have introduced the symbol R , which discards the possible imaginary part, as we are working with complex valued functions. For consistency and to ensure convergence, the acceleration parameter is forced to reside inside the interval 0 < α k < 1 . At each iteration, we start with the current estimation of the phase φ i . Then, we compute the acceleration parameter and gradient to predict the closest possible phase to the next iteration, y k . To this estimation, we apply our algorithm, ψ(y k ) , to obtain the next point phase estimation, φ k+1 .
Experimental implementation. The experimental system used in this work is divided in two parts. First, a beam generator able to produce highly focused fields with arbitrary irradiance and phase distribution (Fig. 1, left blue box). Second, a beam analyzer used to retrieve the transverse Stokes images of the produced beam ( Fig. 1, right blue box). A collimated beam is obtained after lens L 1 placed at a distance equal to focal length from the fiber end of a pig-tailed laser (Thorlabs LP520-SF15@ 520 nm ). The beam is modulated with a translucent twisted-nematic liquid crystal display (Holoeye HEO 0017), with a pixel pitch of 32µm . Linear polarizers LP 1 and LP 2 and quarter wave plate QWP 1 are set to achieve a phase-mostly modulation response 64 able to produce computer generated holograms according to the Arrizon's Double-Pixel Hologram approach 65 . In order to produce radially polarized beams, we set a vortex retarder VR after polarizer LP 2 . This element might be replaced with a QWP or simply removed to generate circularly or linearly polarized beams. Lens L 2 and L 3 form a telecentric (4f) system. A spatial filter is placed at the back focal plane of L 2 . This element is required to remove high orders of the diffracted beam produced by the encoded hologram. An extended description on the use of the spatial filter in combination with Arrizon's holograms can be found in 12,28 . The use and generation of holograms with this method in the context of high NA optical systems, as well as the intrinsic limitations due to the sampling of the modulation points, can be found in 34 . Then, microscope objective MO 1 (Nikon Plan Fluorite N40X-PF with NA= 0.75 ) is used to focus the tailored beam and produce a highly focused beam.
The generated beam is imaged on a CCD camera (Stingray with a 14 bit depth and a pixel pitch of 3.75 µm ) by means of microscope objective MO 2 (Nikon with NA = 0.8) mounted on a movable stage, which is driven by a motorized device (Newport LTA-HL) with a uni-directional repeatability of ±100 nm . Note that MO 2 should www.nature.com/scientificreports/ have a larger NA than the MO 1 in order to be able to collect all the beam. In this way, a set of observation planes separated by 2 µm are recorded. The magnification (M) and resolution of the beam analyzer part (Fig. 1, right blue box) is measured by imaging a 1951 USAF resolution test placed in front of MO 2 , resulting in M=50x and a spatial sampling of 75 nm. As described in "Longitudinal and transverse components of an electromagnetic field" and "Methods" sections, the longitudinal component of the electric field can be inferred from the two complex amplitudes of the transverse electric field, E ⊥ = (E x , E y ) [see Eq. (12)]. This means that the phase recovery algorithm should be independently used for the set of recorded irradiances |E x | 2 and |E y | 2 and thus, phases φ x and φ y can be obtained. Note that analyzers are placed in a zone where the beam can be considered under paraxial conditions. Therefore, LP 3 acts as a projector 66 and can be used to select the proper direction of polarization of E ⊥ ; the use of quarter wave plate QWP 2 is explained in the next paragraph.
Note that there is an arbitrary constant phase factor for each component that has no effect on the propagation of each amplitude alone, i.e. the relative phase delay between the two transverse components is arbitrary. However, since this relative phase is relevant for the longitudinal component estimation, we experimentally retrieve it by recording the corresponding Stokes images. A set of six polarimetric images is taken for each observation plane using quarter wave plate QWP 2 and linear polarizer LP 3 in front of the camera plane. Denoting by β, θ the phase delay introduced by the quarter wave plate and the rotation of the polarizer axis respectively, the Stokes measurements I α,β are: I 0,0 ; I 0,90 ; I 0,45 ; I 0,135 ; I π/2,45 ; I π/2,135 . Finally, the relative phase delay δ is evaluated by whereas the amplitudes |E x | and |E y | are

Results and discussion
Two different beams have been generated to estimate their corresponding longitudinal component. The first one is a radially polarized beam with a vortex phase : where ρ = (x, y, 0) , φ is the azimuth coordinate, and f = 5 mm is the focal length of MO1. NA e is the effective pupil size of the beam which is determined according to the present size of the beam. The second one is a linearly polarized (1,1)-Hermite-Gauss beam: The exact size of the beam at the EP of objective lens MO 1 is difficult to assess. Note that the aperture of each optical element limits the size of the beam as it propagates and the shape of the beam is modified by means of the computer generated hologram displayed on the liquid crystal display. Moreover, the EP of the objective is not physically accessible nor measurable. As the spectra of each component is limited by the physical size of the EP of the microscope, all generated beams are band limited. Then, as the Fourier Transform of the intensity can be demonstrated to contain up to twice the frequency content of the amplitude, we use its size to determine the effective pupil size that each beam uses. This calculation gives us NA e = 0.406 and NA e = 0.379 for the radially polarized vortex and the Hermite-Gauss beams respectively. We have NA= 0.75 for MO1, for reference, meaning that the beams do not completely fill the EP of the objective. The recorded transverse irradiance distributions of these beams after propagation through the optical system, using the electromagnetic propagation theory of Richards and Wolf 57 , is presented in Fig. 2. Figure 3 shows the Stokes parameters for two experimentally measured planes (first and second rows) for the radially polarized vortex beam case. The distance between these two planes is 2 µm . Smaller distances proved to converge to a spherical wave independently on the shape of the beam. Note that the position z 0 is close to the focal plane, but cannot be easily set due to the difficulty to determine the focal plane with enough precision. We will numerically estimate its position as the plane where the energy of the beam is tightly concentrated. We observe that the polarization state of the beam is a complex combination of radial and circular polarizations distributed along the width of the beam. Moreover, we note that the exact polarization state of the beam changes as it propagates, specially near the center of the beam. This can be attributed to the spiral phase of the beam, which curls and uncurls while changing the phase difference between components. Figure 4 (first and second rows) shows beam amplitudes √ I 0,0 and √ I 0,90 at planes z 0 and z 0 + 2 µm . Using the phase retrieval algorithm ("Phase recovery algorithm" section), we obtained the corresponding phases φ x and φ y . However, due to recovering each phase separately, the origin of phases for both φ x and φ y might be different. We have determined this constant random phase difference δ 0 at the maximum of irradiance of the beam max I 0,0 + I 0,90 to be the one given by the Stokes parameters (Fig. 3), δ 0 = arctan S 3 /S 2 . With this information, the beam is propagated up to he focal plane. It can be observed that the focal plane is also the position where the phase has its simplest form, as a series of concentric ramification points. We must also remark that the phase difference in the focal plane is φ y − φ x = 0.062 rad, in contrast with its theoretical value, φ y − φ y = π/2 rad.  To compensate for this experimantally observed discrepancy, we include the phase difference between theory and experiment into the simulations of the beam in Eq. (19).
As the Hermite-Gaussian beam is linearly polarized, calculation of the Stokes parameters is unnecessary, and the amplitude is readily obtained from the irradiance. Figure 5 shows the x and y amplitudes of the experimentally measured beam and the synthetically refocused at the focal plane. The x component, although present, is very weak in comparison with the y component and does not affect the shape of the total irradiance. Its presence can be attributed to a failure of the analyzer to completely absorb the strong y component.

Longitudinal component estimation.
Provided that the complex amplitude of the transverse vector E ⊥ (k ⊥ ; z = 0) is known, Eq. (12) can be used to estimate the longitudinal component E z (k ⊥ ; z = 0) . Figure 6 summarizes the longitudinal component outcomes for the two beams considered. We observe a remarkable agreement between theoretical and experimentally estimated distributions: both show a similar size and shape.
As previously discussed, our beam does not strictly fulfill Eq. (19) due to an unforeseen phase difference introduced by some optical element in the experimental setup. This results in Fig. 6a, b, where two protruding lobes      Figure 6c shows the intensity across the diagonal depicted in red, demonstrating a good agreement between experiment and theory. The experimental distribution is somewhat wider, has stronger secondary maxima and is slightly asymmetric with respect to the diagonal. These effects can be attributed to the imperfections in the hologram at the EP of the optical system, mainly due to the double pixel encoding technique used. Since the Hermite-Gauss beam is linearly polarized [ Fig. 6d-g], there is no need to carry out the Stokes analysis as in the previous case. For this reason, it is noticeable an excellent agreement between theoretical and experimental results. Their size, principal maxima and secondary maxima coincide and are almost equal. Some dissimilarities can be observed in the irradiance in the xy plane, Fig. 6f, as the experimental beam is not perfectly symmetrical with respect to the y = 0 (red) line. The profiles shown in Fig. 6g are also very similar as well. As in the case of the vortex beam, discrepancies can be explained mainly due to the double pixel encoding used in the generation of the hologram. We do not have access to a continuous set of modulation points, but a limited number of them, as we discussed in 34,67 . Moreover, they are not uniformly distributed in the unit circle, which might cause a non-uniform error distribution across the modulation points 12,34 .

Concluding remarks
In this work we have described a method to visualize the longitudinal component of a highly focused beam. The method is based on the characterization of the complex amplitude of the two transverse components of the focused beam in various planes: whereas the irradiance is recorded by means of a camera, the phase is estimated using phase recovery techniques. We realized it is also required to consider the relative phase between the two transverse components. This extra measurement can be carried out by determining the Stokes parameters of the beam. Finally, this information allows us to estimate the longitudinal component with the help of the Gauss theorem. The results obtained in the two cases analyzed show an excellent agreement between theory and experiments.