Efficient calculation of highly focused electromagnetic Schell-model beams

Marcos Aviñoá; Rosario Martínez-Herrero; Artur Carnicer

doi:10.1364/OE.432677

1. Introduction

Three-dimensional electromagnetic field distributions generated in the focal region of a high numerical aperture (NA) focused system has been extensively investigated in the last years [1–28]. Non-paraxial fields have demonstrated very useful in many fields, for instance in confocal microscopy, second and third harmonic generation problems, plasmonic effects, optical tweezers, or optical security, among others. See, for instance, [29–37].

The structure of highly focused beams has a strong dependence on the polarization state of the input beam. Such dependence on polarization has been exploited to generate beams with interesting properties [1,7,8,11,14,18–20,22,24,26,27]. An interesting example is light needles [8,26,28], in which the cross-section of the beam remains confined to a region the size of its wavelength for distances of many wavelengths away from the focal plane of the system.

For incident coherent and totally polarized paraxial beams, the calculations of the focal field structure can be carried out using Fourier transforms. Despite Discrete Fourier Transforms (DFT) have a time complexity of $O(N^{2})$, where $N$ is the total number of points, the so-called Fast Fourier Transform Algorithms (FFT) [38,39], display a time complexity of $O(N\log N)$. However, limitations arise when working with partially coherent and partially polarized incident beams: the statistical properties of light must be taken into account and instead of working with field amplitudes cross-correlation matrices should be taken into account [40]. As a matter of fact, this increases the dimensionality of the problem from 2D to 4D. FFTs could still be performed in theory, but in practice memory limitations impose a hard limit on these algorithms. For instance, when trying to calculate the irradiance of a $256\times 256$ pixels image, each component of the electric field requires at least 32 GiB of RAM of single precision floating point numbers. These memory requirements make most calculations unachievable under most common computer set-ups. To circumvent this problem, modal expansions of the field can be used to reduce the dimensionality from 4D to 2D, at the expense of a higher number of 2D operations [40–42]. Recently, this method has been successfully applied to the study of highly focused electromagnetic beams [43–45]. Indeed, this is a powerful method however, it might be limited by the convergence speed. In theory, each field requires an infinite number of modes. If the weights of each of these modes do not fall off sufficiently fast, the modal representation would require a high number of terms and thus, a large amount of memory. Several other methods relying on the decomposition

In this work we propose a fast method of calculating the focal structure of a highly focalized, partially coherent and partially polarized beam without using a modal decomposition. The incident beam is assumed to be an electromagnetic Schell field. This kind of paraxial fields is one of the more remarkable models of partially coherent and partially polarized beams, mainly due to their simple mathematical form and the wide range of applications [46–49]. Taking advantage of their particular properties, it is possible to reduce the 4D problem to a set of six 2D convolutions, independently on the shape of the beam. Interestingly, this approach has been successfully applied in paraxial propagation conditions [50–52].

The paper is organized as follows. In section 2 we discuss the theory behind the tight focusing of electromagnetic beams both in the coherent and partially coherent and partially polarized regimes. Section 3 introduces the family of electromagnetic Schell-model sources, the properties of which are exploited to reduce the dimensionality of the problem. In section 4 we introduce the algorithm for the fast calculation of focalized electromagnetic Schell source beams. To demonstrate our method, in section 5 we analyze the focalization of a partially coherent beam with no closed analytical form. We also examine its computation time and make some extrapolations to other algorithms. Finally, we summarize the results in section 7.

2. Focal structure of light beams

The electric field distribution at any point in the focal region of a high numerical aperture focusing system is given by the well known Richards-Wolf integral [53]

(1)$$\boldsymbol{E}(x, y, z) = A \iint\limits_{-\infty}^{+\infty} P(\theta) D(\theta, \phi) \boldsymbol{E}_0(\theta, \phi) \mathrm{e}^{ik\sin\theta (x\cos\phi+y\sin\phi)} \mathrm{e}^{{-}ikz\cos\theta} \sin\theta d\theta d\phi,$$

where $A$ is a constant related to the focal length and wavelength, $k$ is the wavenumber, $P(\theta )$ denotes the so called the apodization function (obtained from energy conservation and geometric considerations) and $D(\theta , \phi )$ is the region where $\boldsymbol {E}_0(\theta , \phi )$ differs from zero. Angles $\theta$ and $\phi$ are the coordinates at the exit pupil. The geometry of the problem is schematically presented in Fig. 1. The relation of $\boldsymbol {E}_0(\theta , \phi )$ with the radial and azimuthal components of the incident paraxial beam is

(2)$$\boldsymbol{E}_0 (\boldsymbol{r}) = (\boldsymbol{E}_i\cdot \boldsymbol{e}_\phi)\boldsymbol{e}_\phi+ (\boldsymbol{E}_i\cdot \boldsymbol{e}_\theta ')\boldsymbol{e}_\theta,$$

where the paraxial beam $\boldsymbol {E}_i(\theta , \phi )$ is

(3)$$\boldsymbol{E}_i(\theta, \phi) = \left( E_{ix}(\theta, \phi), E_{iy}(\theta, \phi), 0 \right),$$

and the unit vectors are given by

(4)$$\Bigg\lbrace \begin{aligned} & \boldsymbol{e}_\theta' = (\cos \phi,\ \sin\phi,\ 0)\\ & \boldsymbol{e}_\phi = (-\sin \phi,\ \cos\phi,\ 0)\\ & \boldsymbol{e}_\theta = (\cos\phi \cos\theta,\ \sin\phi\cos\theta,\ -\sin\theta).\end{aligned}$$

Fig. 1. Geometry of the electromagnetic diffraction problem.

Download Full Size | PDF

Here, we consider an apodization function of the form $P(\theta ) = \sqrt {\cos \theta }$, corresponding to an optical system satisfying Abbe’s sine condition. Furthermore, let us consider the change of variables $\rho = \sin \theta$ which, upon introducing into Eq. (1) we obtain

(5)$$\boldsymbol{E}(x, y, z) = A \iint\limits_{-\infty}^{+\infty} D(\rho, \phi) \frac{\boldsymbol{E}_0(\rho, \phi)}{\left(1-\rho^{2}\right)^{1/4}} \mathrm{e}^{ik\rho(x\cos\phi + y\sin\phi)}\rho d\rho d\phi.$$

If we introduce the $2\times 3$ matrix $\hat {\boldsymbol {N}}(\rho , \phi )$

(6)$$\hat{\boldsymbol{{N}}}(\rho, \phi) = \begin{pmatrix} \sin^{2}\phi + \cos^{2} \phi \sqrt{1-\rho^{2}} & \sin\phi \cos\phi \left(\sqrt{1-\rho^{2}}-1\right) & \rho\cos\phi \\ \sin\phi \cos\phi \left(\sqrt{1-\rho^{2}}-1\right) & \cos^{2}\phi + \sin^{2} \phi \sqrt{1-\rho^{2}} & \rho\sin\phi \end{pmatrix},$$

we may write $\boldsymbol {E}_0(\rho , \theta ) = \boldsymbol {E}_i(\rho , \phi ) \hat {\boldsymbol {{N}}}(\rho , \phi )$. Finally, by letting $\boldsymbol {\rho } = (\rho \cos \phi , \rho \sin \phi )$ and $\boldsymbol {r}=(x, y)$ Eq. (5) can be rewritten in the much compact form

(7)$$\boldsymbol{E}(\boldsymbol{r}, z) = A\iint\limits_{-\infty}^{+\infty} \frac{D(\boldsymbol{\rho})}{\left(1-\rho^{2}\right)^{1/4}} \boldsymbol{E}_i(\boldsymbol{\rho}) \hat{\boldsymbol{N}}(\boldsymbol{\rho})\mathrm{e}^{ik\boldsymbol{\rho}\cdot\boldsymbol{r}} \mathrm{e}^{{-}ikz\sqrt{1-\rho^{2}}} d\boldsymbol{\rho}.$$

2.1 Partially coherent and partially polarized beams

Equation (7) describes a focused coherent, totally polarized electromagnetic beam. However, electromagnetic waves have in general an intrinsic randomness associated with them. The description given by Eq. (7) can be insufficient to capture all properties of the field and correlation functions should be used [40]. Specifically, since we are dealing with a single monochromatic component of the electromagnetic radiation, we introduce the cross spectral density matrix (CSDM)

(8)$$\hat{\boldsymbol{W}}(\boldsymbol{r}_1, \boldsymbol{r_2}, z) = \langle \boldsymbol{E}^{{\dagger}}(\boldsymbol{r}_1, z) \boldsymbol{E}(\boldsymbol{r}_2, z)\rangle.$$

Using the focal field in Eq. (7) and introducing it into Eq. (8) we obtain the following expression

(9)$$ \begin{aligned} \hat{W}\left(\boldsymbol{r}_{1}, \boldsymbol{r}_{2}, z\right)=|A|^{2} & \iiiint_{-\infty}^{+\infty} \frac{D\left(\boldsymbol{\rho}_{1}\right) D\left(\boldsymbol{\rho}_{2}\right)}{\left(1-\rho_{1}^{2}\right)^{1 / 4}\left(1-\rho_{2}^{2}\right)^{1 / 4}} \hat{\boldsymbol{N}}^{\dagger}\left(\boldsymbol{\rho}_{1}\right) \hat{\boldsymbol{W}}_{i}\left(\boldsymbol{\rho}_{1}, \boldsymbol{\rho}_{2}\right) \hat{N}\left(\boldsymbol{\rho}_{2}\right) \\ & \times \mathrm{e}^{-i k\left(\boldsymbol{r}_{1} \cdot \boldsymbol{\rho}_{1}-\boldsymbol{r}_{2} \cdot \boldsymbol{\rho}_{2}\right)} \mathrm{e}^{i k z\left(\sqrt{1-\rho_{1}^{2}}-\sqrt{1-\rho_{2}^{2}}\right)} d \boldsymbol{\rho}_{1} d \boldsymbol{\rho}_{2} \end{aligned} $$

where $\dagger$ is the Hermitian conjugate and $\hat {\boldsymbol {W}}_i(\boldsymbol {\rho }_1, \boldsymbol {\rho }_2)$ is the $2\times 2$ CSDM of the input field. With this, the irradiance of the field can be determined by calculating the trace of the CSDM and evaluating it at $\boldsymbol {r} = \boldsymbol {r}_1 = \boldsymbol {r}_2$,

(10)$$I(\boldsymbol{r}, z) = \operatorname{Tr} \left[ \hat{\boldsymbol{W}}(\boldsymbol{r}, \boldsymbol{r}, z)\right].$$

In general, Eq. (9) involves a 4D Fourier Transform with 4D arrays. Memory limits the maximum number of pixels for sampling the field. As explained above, squared arrays of $256$ pixels and single precision complex numbers, would require at least $32$ GiB of memory for a single field component array.

3. Electromagnetic Schell-model sources

The CSDM for an electromagnetic Schell incident beams is described by

(11)$$\hat{\boldsymbol{W}}_i(\boldsymbol{\rho}_1, \boldsymbol{\rho}_2) = \hat{\boldsymbol{\tau}}^{{\dagger}}(\boldsymbol{\rho}_1) \hat{\boldsymbol{W}}_S (\boldsymbol{\rho}_1-\boldsymbol{\rho}_2)\hat{\boldsymbol{\tau}}(\boldsymbol{\rho}_2),$$

where

(12)$$\hat{\boldsymbol{W}}_S(\boldsymbol{\rho}_1-\boldsymbol{\rho}_2) = \iint\limits_{-\infty}^{+\infty} \hat{\boldsymbol{H}}^{{\dagger}} (\boldsymbol{\sigma})\hat{\boldsymbol{H}}(\boldsymbol{\sigma}) \mathrm{e}^{ik\boldsymbol{\sigma}\cdot (\boldsymbol{\rho_1}-\boldsymbol{\rho_2})} d\boldsymbol{\sigma},$$

and $\hat {\boldsymbol {H}}(\boldsymbol {\sigma })$ is an arbitrary matrix which determines the coherence and polarization properties of the incident field, along with the matrix $\hat {\boldsymbol {\tau }}$ [54,55]. Note that the CSDM is the matrix product of three terms, one depending on the difference of coordinates and two which depend on one set of them. Now, we introduce matrix $\hat {\boldsymbol {N}}_0(\boldsymbol {\rho }, z)$ as

(13)$$\hat{\boldsymbol{N}}_0(\boldsymbol{\rho}, z) = \frac{D(\boldsymbol{\rho})}{\left(1-\rho^{2}\right)^{1/4}} \mathrm{e}^{{-}ikz\sqrt{1-\rho^{2}}} \hat{\boldsymbol{\tau}}(\boldsymbol{\rho})\hat{\boldsymbol{N}}(\boldsymbol{\rho}).$$

Then, the CSDM of the focused field [Eq. (9)] can be rewritten as

(14)$$ \hat{W}\left(\boldsymbol{r}_{1}, \boldsymbol{r}_{2}, z\right)=|A|^{2} \iiiint_{-\infty}^{+\infty} \hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{1}, z\right) \hat{\boldsymbol{W}}_{S}\left(\boldsymbol{\rho}_{1}-\boldsymbol{\rho}_{2}\right) \hat{N}_{0}\left(\boldsymbol{\rho}_{2}, z\right) \mathrm{e}^{-i k\left(r_{1} \cdot \boldsymbol{\rho}_{1}-\boldsymbol{r}_{2} \cdot \boldsymbol{\rho}_{2}\right)} d \boldsymbol{\rho}_{1} d \boldsymbol{\rho}_{2} . $$

By considering the following change of variables

(15)$$\bigg\lbrace \begin{aligned} & \boldsymbol{\rho}_1-\boldsymbol{\rho}_2 = \boldsymbol{\Delta}\\ & \boldsymbol{\rho}_1+\boldsymbol{\rho}_2 = 2\boldsymbol{\rho}_0, \end{aligned}$$

Equation (14) becomes

(16)$$ \begin{aligned} \hat{\boldsymbol{W}}\left(\boldsymbol{r}_{1}, \boldsymbol{r}_{2}, z\right)=|A|^{2} & \iiiint_{-\infty}^{+\infty} \hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2, z\right) \hat{\boldsymbol{W}}_{S}(\boldsymbol{\Delta}) \hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2, z\right) \\ & \times \mathrm{e}^{-\frac{i k \Lambda}{2} \cdot\left(r_{1}+\boldsymbol{r}_{2}\right)} \mathrm{e}^{-i k \boldsymbol{\rho}_{0} \cdot\left(\boldsymbol{r}_{1}-\boldsymbol{r}_{2}\right)} d \boldsymbol{\rho}_{0} d \boldsymbol{\Delta} . \end{aligned} $$

In this form, the CSDM is still determined by a 4D Fourier Transform. The analysis of each term of the matrix can help to bring it into a simpler form.

(17)$$ \begin{aligned} \left[\hat{\boldsymbol{W}}\left(\boldsymbol{r}_{1}, \boldsymbol{r}_{2}, z\right)\right]_{i j}=|A|^{2} \sum_{p=1}^{2} \sum_{q=1}^{2} & \iint_{-\infty}^{+\infty}\left[\hat{\boldsymbol{W}}_{S}(\boldsymbol{\Delta})\right]_{q p} \mathrm{e}^{-\frac{i k \Delta}{2} \cdot\left(\boldsymbol{r}_{1}+\boldsymbol{r}_{2}\right)} \\ & \times\left\{\iint_{-\infty}^{+\infty}\left[\hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2\right)\right]_{i q}\left[\hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2\right)\right]_{p j}\right.\\ &\left.\times \mathrm{e}^{i k \boldsymbol{\rho}_{0} \cdot\left(\boldsymbol{r}_{2}-\boldsymbol{r}_{1}\right)} d \boldsymbol{\rho}_{0}\right\} d \boldsymbol{\Delta} \end{aligned} $$

where $i, j = 1, 2, 3$. Note that each term in the $(p, q)$ sum has the form of the 2D Fourier transform of a 4D convolution, as the two exponential terms in the integrals depend each one on the sum and the difference of focal plane coordinates respectively.

Analyzing each term of the CSDM, and bringing our attention to measurements of the CSDM across points $\boldsymbol {r}_1 = \boldsymbol {r}_2 = \boldsymbol {r}$, Eq. (16) reads

(18)$$ \begin{aligned} [\hat{W}(\boldsymbol{r}, \boldsymbol{r}, z)]_{i j}=|A|^{2} \sum_{p=1}^{2} \sum_{q=1}^{2} & \iint_{-\infty}^{+\infty}\left[\hat{\boldsymbol{W}}_{S}(\boldsymbol{\Delta})\right]_{q p} \mathrm{e}^{-i k \boldsymbol{\alpha} \cdot \boldsymbol{r}} \\ & \times\left\{\iint_{-\infty}^{+\infty}\left[\hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2\right)\right]_{i q}\left[\hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2\right)\right]_{p j} d \boldsymbol{\Delta}\right\} d \boldsymbol{\rho}_{0}, \end{aligned} $$

which is just the 2D Fourier Transform of the function $\left [\hat {\boldsymbol {W}}(\boldsymbol {\Delta })\right ]_{qp}$ multiplied by a 2D convolution. Equation (18) provides the relationship between the properties of coherence-polarization of the incident beam and the elements of the polarization matrix of the focused field. This equation appears to be specially suitable for the numerical evaluation of the focused field since it only involves 2D arrays. This means that the calculation of the elements of the CSDM can be performed in a reasonable short amount of time.

4. Algorithm implementation

Now, Eq. (18) is written in an explicit form, suitable for calculations. The CSDM in the focal region consists of nine terms, corresponding to the correlations of the $x, y$ and $z$ components of the electric field, which in turn are the sum of the four elements of $\hat {\boldsymbol {W}}_S(\boldsymbol {\Delta })$. Note that $i\leq j$, as the off-diagonal terms are complex conjugate of each other, $\hat {\boldsymbol {W}}_{ji} = \hat {\boldsymbol {W}}_{ij}^{\star }$, because of the Hermiticity of the CSDM. This reduces the number of terms from nine to six. Unfortunately, although the terms in $\hat {\boldsymbol {W}}_S(\boldsymbol {\Delta })_{p\neq q}$ are complex conjugate of one another, their calculation cannot be simplified in general as matrices $\hat {\boldsymbol {N}}_0$ are not square matrices. Expanding Eq. (18), we obtain:

(19)$$ \begin{aligned} [\hat{\boldsymbol{W}}(\boldsymbol{r}, \boldsymbol{r}, z)]_{i \leq j}=|A|^{2} \iiiint \int_{-\infty}^{+\infty}\{&\left[W_{S}(\boldsymbol{\Delta})\right]_{11}\left[\hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2\right)\right]_{i 1}\left[\hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2\right)\right]_{1 j}+\\ &\left[W_{S}(\boldsymbol{\Delta})\right]_{22}\left[\hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2\right)\right]_{i 2}\left[\hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2\right)\right]_{2 j}+\\ &\left[W_{S}(\boldsymbol{\Delta})\right]_{12}\left[\hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2\right)\right]_{i 1}\left[\hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2\right)\right]_{2 j}+\\ &\left.\left[W_{S}(\boldsymbol{\Delta})\right]_{21}\left[\hat{\boldsymbol{N}}_{0}^{\dagger}\left(\boldsymbol{\rho}_{0}+\boldsymbol{\Delta} / 2\right)\right]_{i 2}\left[\hat{\boldsymbol{N}}_{0}\left(\boldsymbol{\rho}_{0}-\boldsymbol{\Delta} / 2\right)\right]_{1 j}\right\} \\ & \times \mathrm{e}^{-i k \boldsymbol{\Delta} \cdot r} d \boldsymbol{\rho}_{0} d \boldsymbol{\Delta} . \end{aligned} $$

Equation (19) can be numerically evaluated by the use of FFTs and the convolution theorem. As the final Fourier Transform factors outside all convolution sums, we just need to perform it once. An outline of the program structure is presented in Algorithm (1).

The algorithm consists of the computation of six independent terms, composed by the Fourier Transform of four different convolutions. This means that $6\times (4\times 2+1) = 54$ Fourier Transformations are required to completely determine the CSDM. As this is the most complex calculation we are performing, we neglect the time complexity of the rest of the operations. With this, we state that the time complexity of our algorithm is

(20)$$T(N) = \mathrm{O}(54 N\log N),$$

where $N = n_y\times n_x$ is the number of pixels of the amplitude of the electromagnetic field and the so called "big O" $O(\cdot )$ gives the asymptotic time dependence of the algorithm on the number of operations. As the number of operations does not depend on a modal decomposition of the field, we call our method the Mode-less convolution algorithm (Moleca). The code can be found in [56].

5. Numerical example

To assess our algorithm, we analyze a beam with no known closed analytical solution for Eq. (16). The characterization of the beam depends on the shape of the terms of Eq. (11). Let us consider an incident beam with $\hat {\boldsymbol {W}}_S(\boldsymbol {\Delta })$

(21)$$\hat{\boldsymbol{W}}_S(\boldsymbol{\Delta}) = \hat{\boldsymbol{P}} \mathrm{L}_n\left(\frac{\Delta^{2}}{2\mu^{2}}\right) \mathrm{e}^{-\frac{\Delta^{2}}{2\mu^{2}}} ,$$

where $\hat {\boldsymbol {P}}$ is a constant matrix given by

(22)$$\hat{\boldsymbol{P}} = \begin{pmatrix} 1 & 1\\ 1 & 1 \end{pmatrix},$$

$\mathrm {L}_n$ is the Laguerre polynomial of order $n$ and $\mu$ is the so-called correlation length, which provides the distance where the correlations in the field vanishes. Matrix $\hat {\boldsymbol {\tau }}(\boldsymbol {\rho })$ takes the form

(23)$$\hat{\boldsymbol{\tau}} = \begin{pmatrix} \cos\phi & 0\\ 0 & \sin\phi \end{pmatrix}\frac{\rho}{w_0}\mathrm{e}^{-\frac{\rho^{2}}{w_0}}.$$

Here $w_0 = R_{\mathrm {EP}}$ is the width of the beam and $R_{\mathrm {EP}}$ is the radius of the entrance pupil of the microscope objective. The choice made of the above matrices ensures that the CSDM of the incident beam is well-defined [54,55]. Specifically, with the chosen parameters we are dealing with a radially polarized, Laguerre-Gaussian Schell type incident beam. The irradiance distribution and degree of coherence [57,58] at the entrance pupil are represented in Fig. 2.

Fig. 2. (a) Irradiance distribution and (b) degrees of coherence of the beam at the entrance pupil of the microscope objective.

Download Full Size | PDF

The parameters used for the calculation are: NA = 0.95, $N = 512\times 512$ sampling points inside a padded array of $N_p = 1536\times 1536$ and $z=0$ (focal plane); three different coherence lengths have considered: $\mu \rightarrow \infty$, $\mu = w_0/4$ and $\mu =w_0/50$. To assess the polarization at the focal plane of the objective, we consider the 3D Degree of Polarization (DoP) [59]:

(24)$$\mathrm{DoP} = \left[ \frac{3}{2}\left(\frac{\operatorname{Tr}(\hat{\boldsymbol{W}}^{2})}{(\operatorname{Tr} \hat{\boldsymbol{W}})^{2}}-\frac{1}{3}\right)\right]^{1/2}.$$

The results are presented in Fig. 3. Computer and library specifications are summarized in Table (1). The first row of Fig. 3 corresponds to a coherent beam, with $\mu \rightarrow \infty$. We observe that the focalized beam produces a doughnut-like distribution transverse irradiance, while a circularly symmetric spot for the $z$ irradiance. Under these conditions, the latter contains more energy than the former and the total irradiance is mostly driven by the shape of the $z$ component. The DoP is constant and equal to one throughout the region, as it should be for a coherent and totally polarized beam.

Fig. 3. Focalization of a partially coherent radially polarized beam in the EP of a microscope. From top to bottom, the rows correspond to $\mu \rightarrow \infty$, $\mu = w_0/4$ and $\mu = w_0/50$. From left to right, the columns are the intensities in the transverse and longitudinal directions and the DoP through the region.

Download Full Size | PDF

Table 1. Relevant technical specifications of the computer used to perform the calculations.

View Table

The second row contains the results for $\mu = w_0/4$. It can be seen that the decrease in the values of the cross-correlations increases the spread of the irradiance, reducing the amount of light capable to be focalized at a single point. The transverse component remains similar, but broadened. The longitudinal component has completely changed its shape: the central maximum has become a minimum while the outer ring has completely disappeared. Interestingly, with a smaller coherence length the DoP in the region with the higher energy density is reduced, achieving a value $\mathrm {DoP} = 0.33$ for its lowest point. Outside this region, the values increase to $\mathrm {DoP} = 0.78$, but the irradiance is practically null. As the limit $\mathrm {DoP} = 0.5$ corresponds to a 2D unpolarized electromagnetic field [59], the field becomes completely unpolarized in the 2D sense and highly unpolarized in the 3D sense.

Finally, the last row corresponds to $\mu = w_0/50$. The beam notably fails to focalize in any meaningful way, while the irradiances for the transverse and longitudinal components become almost indistinguishable, the latter still containing slightly more energy. The DoP becomes almost homogeneous throughout the region, with a minimum value of $\mathrm {DoP} = 0.43$ and a maximum value of $\mathrm {DoP} = 0.47$. Remarkably, at any point the DoP is higher than the minimum value of the $\mu = w_0/4$ case: the field is less polarized in the latter case than in the former.

To provide more insight on the beam we are analyzing, we computed the irradiance at different planes by modifying the propagation distance $z$ with respect to the focal plane. We observe interesting behaviors of the field considered in Eqs. (22)–(23). Figure 4(a) shows the transverse irradiance for the three coherence lengths considered at the center of the beam, $\boldsymbol {r}=0$. The coherent beam, $\mu \rightarrow \infty$ remains null throughout the region considered. The mostly incoherent beam, $\mu = w_0/50$ remains with a mostly constant irradiance through the considered region. On the other hand, the partially coherent beam, $\mu = w_0/4$, has a local minimum in the focal plane and two maxima as the beam propagates outside this plane.

Fig. 4. Analysis of the irradiance at $\boldsymbol {r}=0$: (a) Transverse irradiance, (b) axial irradiance (c) total irradiance and (d) DoP of the beam at $\boldsymbol {r}=0$. Each of the irradiances is normalized to their respective maximum total irradiance.

Download Full Size | PDF

Figure 4(b) shows the irradiance corresponding to the longitudinal component of the electric field. For the coherent case it contains a global maximum in the focal plane, and rapidly decreases as it propagates. For the quasi incoherent beam, the irradiance shows a similar slow variation as in the transverse irradiance but with double the energy. The partially coherent beam contains a similar amount of energy in both the longitudinal and the transverse component, but with a shallower minimum in the focal plane.

The total irradiance is plotted in Fig. 4(c). We observe that the coherent beam sharply focalizes in a region of the order of the wavelength; the mostly incoherent beam fails to focalize because its intensity remaining constant, and the partially coherent display a non-zero minimum in the focal plane.

The DoP has radically different shapes for the different coherence lengths for the input beam, as shown in Fig. 4(d). As the shape of the polarization matrix at the center of the beam is diagonal, the DoP depends only on the relationship between the transverse and longitudinal irradiances. For the coherent beam, it remains constant and equal to one throughout the region, as it has only one nonzero value in its polarization matrix. For the case of the mostly incoherent beam, it also remains mostly constant but with a lower value, around $\mathrm {DoP} = 0.48$. Interestingly, for the partially coherent case, the $\mathrm {DoP}$ is higher than the mostly incoherent at the focal plane, but quickly drops its value below it.

6. Comparative of computation time

Finally, we assess the speed of our algorithm by comparing it with the naive 4D convolution and the decomposition suggested in [45]. Results are present in Fig. (5). We used the following lateral number of points: $n = 64,\ 128,\ 256,\ 1024,$ and $2048$ and, accordingly, the total number of sampling points is $N = n^{2}$. Since the estimated time complexity our proposed Moleca method [Eq. (19) follows Eq. (20), we expect it to have an approximate dependence

(25)$$T_{\mathrm{convolution}}(N) = A{N\log N} = A{n^{2}\log n^{2}} \approx A n^{k}$$

for moderately low values of $n$. $A$ and $k$ can be easily determined by fitting the yellow points in Fig. (5) (blue line): $A = 4.7\cdot 10^{-6}\ \mathrm {s}, \ k = 2.07$ Note the exponent in the lateral number of pixels is slightly higher than 2 because of the $\log n^{2}$ term. Due to the $\mathrm {O}(n^{4})$ dependence of the naive algorithm, results cannot feasibly be reproduced. However, since we know that for our algorithm the number of Fourier Transforms performed is $54 = 6\times (4\times 2+1)$, we can use the constant $A$ as a basis to obtain the raw performance of our computer $C_0$. We define it by dividing the constant $A$ previously obtained by the number of Fourier Transforms performed, resulting in $C_0 = A/54 \approx 9.4\cdot 10^{-8}$. Using this number, the rough estimate for an $n = 2048$ image would approximately take

(26)$$T_{\mathrm{naive}}(n) \approx C_0 n^{4} \approx 17 \ \mathrm{days}$$

which is an exceedingly long computation time. At any point, our method results at least two orders of magnitude faster than the naive method.

Finally, using a similar approach, we have estimated the performance of an alternative method based in a modal decomposition [45]. This method is based on a 2D mode decomposition of the CSDM at the EP of the optical system. Consequently, the time complexity should be of the order of the FFT algorithm, multiplied by the number of modes required in the decomposition. This gives us a time complexity of the form

(27)$$T_{\mathrm{alt}}(n, m) = \mathrm{O}(m n^{2}\log n^{2}) \approx A m n^{k},$$

which is similar to the time complexity of our algorithm, but with a variable number of FFTs. As a compromise between speed and accuracy, in the article they proposed $m=2000$ modes in the decomposition for $N = 512\times 512$. As this depends on the number of sampling points, we have scaled the number $m$ proportionally with the number of sampling points. With this, we obtain the red points in Fig. (5). We see that the modes-based algorithm is faster than Moleca for a small number of points, but its slope increases as the number of sampling points is increased. In fact, for a worse case completely incoherent field, its time complexity could arrive at $T_{\mathrm {alt}}(n, n^{2}) = \mathrm {O}(n^{4}\log n^{2})$ which is worse than the naive algorithm. However, if the number of modes remains lower than $54$, this algorithm outperforms Moleca.

Fig. 5. Calculation times for three different algorithms for the determination of the partially coherent focal structure of an electromagnetic beam.

Download Full Size | PDF

7. Concluding remarks

We proposed a computational procedure for calculating the CSDM $W(\boldsymbol {r}, \boldsymbol {r}, z)$ in the focal region of a high NA optical system, assuming electromagnetic Schell input beams. The method makes the most of the mathematical properties of these kinds of beams so that the 4D integrals involved can be rewritten as 2D Fourier transforms of 2D convolutions. As a result, we obtain the Moleca algorithm. Contrary to the methods in the literature, our approach does not rely on a modal decomposition of the input beam, requiring a fixed and low number of operations to compute the CSDM of highly focused beams. The time complexity is demonstrated to be of the form $O(54N\log N)$, where $54 = 6\times (4\times 2+1)$ is the number of Fourier Transforms required in the algorithm and $N = n_x\times n_y$ is the number of pixels in which we sample the beam. Comparison with known algorithms show the clear advantage of Moleca with respect to the naive, direct 4D convolution method, and a relative advantage with respect to modal decomposition approaches as long as the number of required modes is less than 54.

In summary, our method is able to calculate the polarization matrix of highly focused electromagnetic Schell-source beams. Despite our approach is faster, the modal decomposition method can calculate both the polarization and the coherence properties of the propagated beam even for non-Schell sources.

Funding

Agencia Estatal de Investigación (PID2019-104268GB-C21, PID2019-104268GB-C22).

Disclosures

The authors declare no conflicts of interest.

Data availability

The data generated in Fig. 2 and 3 as well as the code used to generated is included in the example in [60]. The data points in Fig. 5 and 4 are not publicly available.

References

1. R. Dorn, S. Quabis, and G. Leuchs, “Sharper focus for a radially polarized light beam,” Phys. Rev. Lett. 91(23), 233901 (2003). [CrossRef]

2. K. S. Youngworth and T. G. Brown, “Focusing of high numerical aperture cylindrical-vector beams,” Opt. Express 7(2), 77–87 (2000). [CrossRef]

3. M. Leutenegger, R. Rao, R. A. Leitgeb, and T. Lasser, “Fast focus field calculations,” Opt. Express 14(23), 11277–11291 (2006). [CrossRef]

4. Y. Kozawa and S. Sato, “Sharper focal spot formed by higher-order radially polarized laser beams,” J. Opt. Soc. Am. A 24(6), 1793–1798 (2007). [CrossRef]

5. C. Maurer, A. Jesacher, S. Fürhapter, S. Bernet, and M. Ritsch-Marte, “Tailoring of arbitrary optical vector beams,” New J. Phys. 9(3), 78 (2007). [CrossRef]

6. Y. Zhao, J. S. Edgar, G. D. M. Jeffries, D. McGloin, and D. T. Chiu, “Spin-to-orbital angular momentum conversion in a strongly focused optical beam,” Phys. Rev. Lett. 99(7), 073901 (2007). [CrossRef]

7. K. Lindfors, A. Priimagi, T. Setälä, A. Shevchenko, A. Friberg, and M. Kaivola, “Local polarization of tightly focused unpolarized light,” Nat. Photonics 1(4), 228–231 (2007). [CrossRef]

8. H. Wang, L. Shi, C. Lukyanchuk, and C. Chong, “Creation of a needle of longitudinally polarized light in vacuum using binary optics,” Nat. Photonics 2(8), 501–505 (2008). [CrossRef]

9. G. M. Lerman and U. Levy, “Effect of radial polarization and apodization on spot size under tight focusing conditions,” Opt. Express 16(7), 4567–4581 (2008). [CrossRef]

10. M. R. Foreman and P. Török, “Focusing of spatially inhomogeneous partially coherent, partially polarized electromagnetic fields,” J. Opt. Soc. Am. A 26(11), 2470–2479 (2009). [CrossRef]

11. X.-L. Wang, Y. Li, J. Chen, C.-S. Guo, J. Ding, and H.-T. Wang, “A new type of vector fields with hybrid states of polarization,” Opt. Express 18(10), 10786–10795 (2010). [CrossRef]

12. R. Martínez-Herrero and P. M. Mejías, “Angular momentum decomposition of nonparaxial light beams,” Opt. Express 18(8), 7965–7971 (2010). [CrossRef]

13. F. Kenny, D. Lara, O. G. Rodríguez-Herrera, and C. Dainty, “Complete polarization and phase control for focus-shaping in high-na microscopy,” Opt. Express 20(13), 14015–14029 (2012). [CrossRef]

14. J. A. Fan, K. Bao, J. B. Lassiter, J. Bao, N. J. Halas, P. Nordlander, and F. Capasso, “Near-normal incidence dark-field microscopy: Applications to nanoplasmonic spectroscopy,” Nano Lett. 12(6), 2817–2821 (2012). [CrossRef]

15. Y. Dong, F. Wang, C. Zhao, and Y. Cai, “Effect of spatial coherence on propagation, tight focusing, and radiation forces of an azimuthally polarized beam,” Phys. Rev. A 86(1), 013840 (2012). [CrossRef]

16. J. Lin, O. G. Rodríguez-Herrera, F. Kenny, D. Lara, and J. C. Dainty, “Fast vectorial calculation of the volumetric focused field distribution by using a three-dimensional fourier transform,” Opt. Express 20(2), 1060–1069 (2012). [CrossRef]

17. R. Martínez-Herrero, I. Juvells, and A. Carnicer, “On the physical realizability of highly focused electromagnetic field distributions,” Opt. Lett. 38(12), 2065–2067 (2013). [CrossRef]

18. D. Maluenda, R. Martínez-Herrero, I. Juvells, and A. Carnicer, “Synthesis of highly focused fields with circular polarization at any transverse plane,” Opt. Express 22(6), 6859–6867 (2014). [CrossRef]

19. R. Martínez-Herrero, I. Juvells, and A. Carnicer, “Design of highly focused fields that remain unpolarized on axis,” Opt. Lett. 39(20), 6025–6028 (2014). [CrossRef]

20. R. Martínez-Herrero, D. Maluenda, I. Juvells, and A. Carnicer, “Experimental implementation of tightly focused beams with unpolarized transversal component at any plane,” Opt. Express 22(26), 32419–32428 (2014). [CrossRef]

21. T. Bauer, S. Orlov, U. Peschel, P. Banzer, and G. Leuchs, “Nanointerferometric amplitude and phase reconstruction of tightly focused vector beams,” Nat. Photonics 8(1), 23–27 (2014). [CrossRef]

22. T. Bauer, P. Banzer, E. Karimi, S. Orlov, A. Rubano, L. Marrucci, E. Santamato, R. Royd, and G. Leuchs, “Observation of optical polarization m obius strips,” Science 347(6225), 964–966 (2015). [CrossRef]

23. S. Roy, K. Ushakova, Q. van den Berg, S. F. Pereira, and H. P. Urbach, “Radially polarized light for detection and nanolocalization of dielectric particles on a planar substrate,” Phys. Rev. Lett. 114(10), 103903 (2015). [CrossRef]

24. T. Bauer, M. Neugebauer, G. Leuchs, and P. Banzer, “Optical polarization möbius strips and points of purely transverse spin density,” Phys. Rev. Lett. 117(1), 013601 (2016). [CrossRef]

25. R. Martínez-Herrero, D. Maluenda, I. Juvells, and A. Carnicer, “Polarisers in the focal domain: Theoretical model and experimental validation,” Sci. Rep. 7(1), 42122 (2017). [CrossRef]

26. R. Martínez-Herrero, D. Maluenda, I. Juvells, and A. Carnicer, “Synthesis of light needles with tunable length and nearly constant irradiance,” Sci. Rep. 8(1), 2657 (2018). [CrossRef]

27. R. Martínez-Herrero, D. Maluenda, I. Juvells, and A. Carnicer, “Effect of linear polarizers on the behavior of partially coherent and partially polarized highly focused fields,” Opt. Lett. 43(14), 3445–3448 (2018). [CrossRef]

28. H. Hu, Q. Gan, and Q. Zhan, “Generation of a nondiffracting superchiral optical needle for circular dichroism imaging of sparse subdiffraction objects,” Phys. Rev. Lett. 122(22), 223901 (2019). [CrossRef]

29. G. Brakenhoff, P. Blom, and P. Barends, “Confocal scanning light microscopy with high aperture immersion lenses,” J. Microsc. 117(2), 219–232 (1979). [CrossRef]

30. D. Biss and T. Brown, “Polarization-vortex-driven second-harmonic generation,” Opt. Lett. 28(11), 923–925 (2003). [CrossRef]

31. D. Oron, E. Tal, and Y. Silberberg, “Depth-resolved multiphoton polarization microscopy by third-harmonic generation,” Opt. Lett. 28(23), 2315–2317 (2003). [CrossRef]

32. C. J. Sheppard and A. Choudhury, “Annular pupils, radial polarization, and superresolution,” Appl. Opt. 43(22), 4322–4327 (2004). [CrossRef]

33. K. C. Neuman and S. M. Block, “Optical trapping,” Rev. Sci. Instrum. 75(9), 2787–2809 (2004). [CrossRef]

34. K. Serrels, E. Ramsay, R. Warburton, and D. Reid, “Nanoscale optical microscopy in the vectorial focusing regime,” Nat. Photonics 2(5), 311–314 (2008). [CrossRef]

35. Y. Gorodetski, A. Niv, V. Kleiner, and E. Hasman, “Observation of the spin-based plasmonic effect in nanoscale structures,” Phys. Rev. Lett. 101(4), 043903 (2008). [CrossRef]

36. L. Vuong, A. Adam, J. Brok, P. Planken, and H. Urbach, “Electromagnetic spin-orbit interactions via scattering of subwavelength apertures,” Phys. Rev. Lett. 104(8), 083903 (2010). [CrossRef]

37. A. Carnicer, I. Juvells, B. Javidi, and R. Martínez-Herrero, “Optical encryption in the longitudinal domain of focused fields,” Opt. Express 24(7), 6793–6801 (2016). [CrossRef]

38. P. Swarztrauber, Vectorizing the FFTs, in Parallel Computations, (Academic, 1998), pp. 51–83.

39. M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93(2), 216–231 (2005). [CrossRef]

40. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge University, 1995), 1st ed.

41. J. Tervo, T. Setälä, and A. T. Friberg, “Theory of partially coherent electromagnetic fields in the space–frequency domain,” J. Opt. Soc. Am. A 21(11), 2205–2215 (2004). [CrossRef]

42. R. Martínez-Herrero, P. M. Mejías, and F. Gori, “Genuine cross-spectral densities and pseudo-modal expansions,” Opt. Lett. 34(9), 1399–1401 (2009). [CrossRef]

43. X. Xiao and D. Voelz, “Wave optics simulation approach for partial spatially coherent beams,” Opt. Express 14(16), 6986–6992 (2006). [CrossRef]

44. B. J. Davis, “Simulation of vector fields with arbitrary second-order correlations,” Opt. Express 15(6), 2837–2846 (2007). [CrossRef]

45. R. Tong, Z. Dong, Y. Chen, F. Wang, Y. Cai, and T. Setälä, “Fast calculation of tightly focused random electromagnetic beams: controlling the focal field by spatial coherence,” Opt. Express 28(7), 9713–9727 (2020). [CrossRef]

46. F. Gori, “Matrix treatment for partially polarized, partially coherent beams,” Opt. Lett. 23(4), 241–243 (1998). [CrossRef]

47. E. Wolf, “Unified theory of coherence and polarization of random electromagnetic beams,” Phys. Lett. A 312(5-6), 263–267 (2003). [CrossRef]

48. O. Korotkova, X. Chen, and T. Setälä, “Electromagnetic schell-model beams with arbitrary complex correlation states,” Opt. Lett. 44(20), 4945–4948 (2019). [CrossRef]

49. M. W. Hyde, “Generating electromagnetic schell-model sources using complex screens with spatially varying auto- and cross-correlation functions,” Results Phys. 15, 102663 (2019). [CrossRef]

50. Y. Chen, F. Wang, C. Zhao, and Y. Cai, “Experimental demonstration of a laguerre-gaussian correlated schell-model vortex beam,” Opt. Express 22(5), 5826–5838 (2014). [CrossRef]

51. Y. Chen, J. Gu, F. Wang, and Y. Cai, “Self-splitting properties of a hermite-gaussian correlated schell-model beam,” Phys. Rev. A 91(1), 013823 (2015). [CrossRef]

52. F. Wang and O. Korotkova, “Convolution approach for beam propagation in random media,” Opt. Lett. 41(7), 1546–1549 (2016). [CrossRef]

53. B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. structure of the image field in an aplanatic system,” Proc. R. Soc. Lond. A 253(1274), 358–379 (1959). [CrossRef]

54. R. Martínez-Herrero and P. M. Mejías, “Elementary-field expansions of genuine cross-spectral density matrices,” Opt. Lett. 34(15), 2303–2305 (2009). [CrossRef]

55. F. Gori, V. Ramírez-Sánchez, M. Santarsiero, and T. Shirai, “On genuine cross-spectral density matrices,” J. Opt. A: Pure Appl. Opt. 11(8), 085706 (2009). [CrossRef]

56. M. Avi noá, “Moleca (mode-less convolution algorithm),” https://github.com/kramos966/moleca. Accessed: 2021-07-20.

57. T. Setälä, J. Tervo, and A. T. Friberg, “Complete electromagnetic coherence in the space–frequency domain,” Opt. Lett. 29(4), 328–330 (2004). [CrossRef]

58. R. Martínez-Herrero and P. M. Mejías, “Relation between degrees of coherence for electromagnetic fields,” Opt. Lett. 32(11), 1504–1506 (2007). [CrossRef]

59. T. Setälä, A. Shevchenko, M. Kaivola, and A. T. Friberg, “Degree of polarization for optical near fields,” Phys. Rev. E 66(1), 016615 (2002). [CrossRef]

60. K. ramos, “Moleca (Mode-less convolution algorithm) is a numerical method to compute the complex spectral density matrix (CSDM) of a well collimated beam through an optical system with a high Numerical Aperture,” GitHub (2021) [accessed 26 May 2021], https://github.com/kramos966/moleca.

Efficient calculation of highly focused electromagnetic Schell-model beams

Abstract

1. Introduction

2. Focal structure of light beams

2.1 Partially coherent and partially polarized beams

3. Electromagnetic Schell-model sources

4. Algorithm implementation

5. Numerical example

6. Comparative of computation time

7. Concluding remarks

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (5)

Tables (1)

Equations (27)

Optics Express

Marcos Aviñoá	https://orcid.org/0000-0002-6458-3095
Rosario Martínez-Herrero	https://orcid.org/0000-0001-7810-297X
Artur Carnicer	https://orcid.org/0000-0002-4936-5778