Formation of Inverse Energy Flux in the Case of Diffraction of Linearly Polarized Radiation by Conventional and Generalized Spiral Phase Plates

Abstract

Recently, there has been increased interest in the shaping of light fields with an inverse energy flux to guide optically trapped nano- and microparticles towards a radiation source. To generate inverse energy flux, non-uniformly polarized laser beams, especially higher-order cylindrical vector beams, are widely used. Here, we demonstrate the use of conventional and so-called generalized spiral phase plates for the formation of light fields with an inverse energy flux when they are illuminated with linearly polarized radiation. We present an analytical and numerical study of the longitudinal and transverse components of the Poynting vector. The conditions for maximizing the negative value of the real part of the longitudinal component of the Poynting vector are obtained.

1. Introduction

The ability to shape light fields with the desired complex structure, so-called structured laser beams, plays an important role in the field of laser manipulation for implementing optical tweezers with advanced functionality. Conventional optical tweezers are unique tools that use a strongly focused Gaussian laser beam for the trapping and three-dimensional confinement of nano- and microparticles. Structured optical tweezers, also called holographic optical tweezers (HOT), provide more possibilities for laser manipulation due to the possibility of controlling the local distribution of the optical forces acting on the particles located in such a light field. For example, trapped particles can be guided along with the desired two- or three-dimensional trajectory or rotated with a desired frequency depending on the parameters of the used HOT. Recently, laser beams with an inverse energy flux have attracted researchers’ attention. Such light fields allow the optically trapped particle to be pulled towards the radiation source, i.e., in the direction opposite to field propagation. Researchers have studied the formation of such advanced optical tweezers using light fields with phase or polarization singularities. Currently, the most popular element for generating light fields with phase singularities is probably the conventional spiral phase plate (SPP) [1], an optical element with a complex transmission function exp(imφ), where φ is the polar angle and m is the topological charge of the formed vortex beam. The wide popularity and spread of this element have led to a growth of studies on its various modifications and generalizations to obtain new effective elements for forming light fields with phase singularities with new properties. In particular, if conventional SPPs with an integer value of m form vortex beams with an annular intensity distribution [2,3], then SPPs with a fractional value of m [4] allow more complex intensity distributions [5] and can be used for additional information coding [6]. A power-exponent phase plate (PEPP) [7,8] with the transmission function described as exp[i2πm(φ/2π)^s], where m and s are arbitrary numbers and φ is an azimuthal angle, can be used for shaping a spiral intensity distribution. In 2020, we showed that a generalized spiral phase plate (GSPP) [9,10] with the transmission function exp[ig(φ)], where g(φ) is an arbitrary smooth monotonic function, always forms a spiral intensity distribution with a phase singularity point in the focal plane of a lens. Variations in the phase function g(φ) affect the size and the twisting degree of the generated light spiral, thereby determining the phase gradient along the light curve, which is an important light-field characteristic for laser manipulation [11,12,13]. The transverse energy flow density, which is proportional to the product of the field intensity and phase gradient [14,15], allows guiding the trapped particle in a certain direction. However, the longitudinal distribution of the energy flow density is a key for realizing light fields with an inverse energy flux. The inverse energy flux phenomena [16,17,18,19,20] is related to “tractor beams” [21,22,23] in optical trapping and manipulation and to the detection of “invisibility cloaks” [24]. It should be noted that just one inverse energy flux does not guarantee the existence of the pulling optical force—the optical force can be negative due to the negative Poynting vector component, but its effect is often blocked by the gradient-phase force [25,26,27,28].

Earlier, the possibility of forming an inverse energy flux in the focal region for cylindrical vector beams (CVBs) of various orders was investigated, as was a linearly polarized beam with a cosine amplitude distribution [17,18]. Despite the existence of many methods for forming CVBs, linearly polarized laser sources remain the most readily available, affordable, and widespread. Therefore, in this work, we investigate the possibility of forming inverse energy flux in the focal region for a linearly polarized beam under diffraction by both conventional and generalized SPPs. Due to the generation of non-annular light distributions, GSPPs provide more opportunities to control the distribution of optical forces acting on the trapped particle. In particular cases, we carry out analytical and numerical studies of both the longitudinal and transverse components of the Poynting vector, including both the real and the imaginary part. Recently, the imaginary part of this component has attracted attention from researchers [29,30,31,32], since it was shown to contribute to the transverse force acting on optically trapped Mie particles [29,30].

2. Methods

Consider the focused electromagnetic field in the focal region in the framework of the Debye theory [33,34,35]:

(1)$\begin{array}{c}\left[\begin{array}{c}E(r,\gamma ,z)\\ H(r,\gamma ,z)\end{array}\right]=\\ -\frac{if}{\lambda }{\displaystyle \underset{0}{\overset{{\mathsf{\theta }}_{\max }}{\int }}{\displaystyle \underset{0}{\overset{2\pi }{\int }}T(\mathsf{\theta })F(\mathsf{\theta },\mathsf{\phi })\left[\begin{array}{c}{Q}_E(\mathsf{\theta },\mathsf{\phi })\\ Q_H(\mathsf{\theta },\mathsf{\phi })\end{array}\right]\exp \left\{ik\left[r\sin \mathsf{\theta }\cos (\mathsf{\phi }-\gamma )+z\cos \mathsf{\theta }\right]\right\}}}\sin \mathsf{\theta }d\mathsf{\theta }d\mathsf{\phi },\end{array}$

where (r, γ, z) are the cylindrical coordinates in the focal region, (θ, φ) are the spherical angular coordinates of the focusing system’s output pupil, θ_max is the maximum value of the azimuthal angle related to the system’s numerical aperture, F(θ, φ) is the complex amplitude of the input field, T(θ) is the pupil’s apodization function ($T(\mathsf{\theta })=\sqrt{\cos \mathsf{\theta }}$ for aplanatic systems), k = 2π/λ is the wavenumber, λ is the wavelength, f is the focal length, and (c_x(θ, φ), c_y(θ, φ))^T is the vector of polarization coefficients of the input field,

(2)$\begin{array}{l}{Q}_E(\mathsf{\theta },\mathsf{\phi })=\left[\begin{array}{cc}A(\mathsf{\theta },\mathsf{\phi })& C(\mathsf{\theta },\mathsf{\phi })\\ C(\mathsf{\theta },\mathsf{\phi })& B(\mathsf{\theta },\mathsf{\phi })\\ -D(\mathsf{\theta },\mathsf{\phi })& -E(\mathsf{\theta },\mathsf{\phi })\end{array}\right]\left(\begin{array}{c}{c}_x(\mathsf{\theta },\mathsf{\phi })\\ c_y(\mathsf{\theta },\mathsf{\phi })\end{array}\right),\\ Q_H(\mathsf{\theta },\mathsf{\phi })=\left[\begin{array}{cc}C(\mathsf{\theta },\mathsf{\phi })& -A(\mathsf{\theta },\mathsf{\phi })\\ B(\mathsf{\theta },\mathsf{\phi })& -C(\mathsf{\theta },\mathsf{\phi })\\ -E(\mathsf{\theta },\mathsf{\phi })& D(\mathsf{\theta },\mathsf{\phi })\end{array}\right]\left(\begin{array}{c}{c}_x(\mathsf{\theta },\mathsf{\phi })\\ c_y(\mathsf{\theta },\mathsf{\phi })\end{array}\right),\end{array}$

with

(3)$\begin{array}{l}A(\mathsf{\theta },\mathsf{\phi })=1-{\cos }^2\mathsf{\phi }(1-\cos \mathsf{\theta }),\\ B(\mathsf{\theta },\mathsf{\phi })=1-{\sin }^2\mathsf{\phi }(1-\cos \mathsf{\theta }),\\ C(\mathsf{\theta },\mathsf{\phi })=-\sin \mathsf{\phi }\cos \mathsf{\phi }(1-\cos \mathsf{\theta }),\\ D(\mathsf{\theta },\mathsf{\phi })=\cos \mathsf{\phi }\sin \mathsf{\theta },\\ E(\mathsf{\theta },\mathsf{\phi })=\sin \mathsf{\phi }\sin \mathsf{\theta }.\end{array}$

The complex Poynting vector [29,30,31,32] is proportional to the following expression:

(4)$P\simeq E^{*}\times H.$

The real part of expression (4) is a standard Poynting vector, and the imaginary part is interesting, as it relates to exciting non-standard forces arising at the interaction between electric- and magnetic-induced dipoles [29].

In this work, we focus on the longitudinal component of the Poynting vector, which is proportional to the following expression [36]:

(5)$P_z\simeq E_x^{*}{H}_y-E_y^{*}{H}_x.$

Note that Re[P_z] is the most important for our research, since it allows us to determine the areas of reverse energy flow. In addition, the azimuthal component of the Poynting vector Re[P_γ], which is related to the longitudinal component of the angular momentum, is also interesting [37,38]:

(6)$j_z\simeq \mathrm{Re}{\left[r\times P\right]}_z=\mathrm{Re}\left[x{P}_y-y{P}_x\right]=r\mathrm{Re}\left[P_y\cos \gamma -P_x\sin \gamma \right]=r\mathrm{Re}\left[P_{\gamma }\right].$

Further, we consider the relative magnitudes of negative and positive values, so constants and common factors can be ignored. To obtain explicit analytical expressions, consider an input field with a superimposed narrow annular aperture with a central angle θ₀ and width Δ [39,40]:

(7)$F(\mathsf{\theta },\mathsf{\phi })=F({\mathsf{\theta }}_0,\mathsf{\phi })=S\left({\mathsf{\theta }}_0\right)F_0(\mathsf{\phi })$

where F₀(φ) is an arbitrary function and

(8)$S\left({\mathsf{\theta }}_0\right)=\{\begin{array}{c}1,\hspace{0.17em}\hspace{0.17em}\left|\mathsf{\theta }-{\mathsf{\theta }}_0\right|\le \Delta /2,\\ 0,\hspace{0.17em}\hspace{0.17em}else.\end{array}$

In this case, the field components in Equation (4) are calculated in the focal plane (z = 0) by the following formulas:

(9)$\begin{array}{l}\left(\begin{array}{c}{E}_x(r,\gamma ,z)\\ E_y(r,\gamma ,z)\\ E_z(r,\gamma ,z)\end{array}\right)\approx t\left({\mathsf{\theta }}_0,z\right){\displaystyle \underset{0}{\overset{2\pi }{\int }}F({\mathsf{\theta }}_0,\mathsf{\phi })\left[\begin{array}{cc}A({\mathsf{\theta }}_0,\mathsf{\phi })& C({\mathsf{\theta }}_0,\mathsf{\phi })\\ C({\mathsf{\theta }}_0,\mathsf{\phi })& B({\mathsf{\theta }}_0,\mathsf{\phi })\\ -D({\mathsf{\theta }}_0,\mathsf{\phi })& -E({\mathsf{\theta }}_0,\mathsf{\phi })\end{array}\right]\left(\begin{array}{c}{c}_x({\mathsf{\theta }}_0,\mathsf{\phi })\\ c_y({\mathsf{\theta }}_0,\mathsf{\phi })\end{array}\right)}G\left(r,(\mathsf{\phi }-\gamma )\right)d\mathsf{\phi },\\ \left(\begin{array}{c}{H}_x(r,\gamma ,z)\\ H_y(r,\gamma ,z)\\ H_z(r,\gamma ,z)\end{array}\right)\approx t\left({\mathsf{\theta }}_0,z\right){\displaystyle \underset{0}{\overset{2\pi }{\int }}F({\mathsf{\theta }}_0,\mathsf{\phi })\left[\begin{array}{cc}C({\mathsf{\theta }}_0,\mathsf{\phi })& -A({\mathsf{\theta }}_0,\mathsf{\phi })\\ B({\mathsf{\theta }}_0,\mathsf{\phi })& -C({\mathsf{\theta }}_0,\mathsf{\phi })\\ -E({\mathsf{\theta }}_0,\mathsf{\phi })& D({\mathsf{\theta }}_0,\mathsf{\phi })\end{array}\right]\left(\begin{array}{c}{c}_x({\mathsf{\theta }}_0,\mathsf{\phi })\\ c_y({\mathsf{\theta }}_0,\mathsf{\phi })\end{array}\right)}G\left(r,(\mathsf{\phi }-\gamma )\right)d\mathsf{\phi },\end{array}$

where $t({\mathsf{\theta }}_0,z)=-i\exp (ikz\cos {\mathsf{\theta }}_0)\sqrt{\cos {\mathsf{\theta }}_0}\sin {\mathsf{\theta }}_0\cdot \Delta \cdot f/\lambda$, $G(r,\mathsf{\phi })=\exp (ikr\sin {\mathsf{\theta }}_0\cos \mathsf{\phi })$.

3. Results

As an input field F(θ₀, φ), consider a plane x-linearly polarized laser beam with the polarization vector (c_x(θ₀, φ), c_y(θ₀, φ))^T = (1, 0)^T, passed through a GSPP [9,10]:

(10)$F_0(\mathsf{\phi })=\exp \left[ig(\mathsf{\phi })\right],$

where g(φ) is an arbitrary smooth function.

In this case, Equation (9) is simplified and takes the following form:

(11)$\begin{array}{l}\left(\begin{array}{c}{E}_x(r,\gamma ,z)\\ E_y(r,\gamma ,z)\\ E_z(r,\gamma ,z)\end{array}\right)\approx t\left({\mathsf{\theta }}_0,z\right){\displaystyle \underset{0}{\overset{2\pi }{\int }}\exp \left[ig(\mathsf{\phi })\right]\left(\begin{array}{c}A({\mathsf{\theta }}_0,\mathsf{\phi })\\ C({\mathsf{\theta }}_0,\mathsf{\phi })\\ -D({\mathsf{\theta }}_0,\mathsf{\phi })\end{array}\right)}G\left(r,(\mathsf{\phi }-\gamma )\right)d\mathsf{\phi },\\ \left(\begin{array}{c}{H}_x(r,\gamma ,z)\\ H_y(r,\gamma ,z)\\ H_z(r,\gamma ,z)\end{array}\right)\approx t\left({\mathsf{\theta }}_0,z\right){\displaystyle \underset{0}{\overset{2\pi }{\int }}\exp \left[ig(\mathsf{\phi })\right]\left(\begin{array}{c}C({\mathsf{\theta }}_0,\mathsf{\phi })\\ B({\mathsf{\theta }}_0,\mathsf{\phi })\\ -E({\mathsf{\theta }}_0,\mathsf{\phi })\end{array}\right)}G\left(r,(\mathsf{\phi }-\gamma )\right)d\mathsf{\phi }.\end{array}$

Equation (11) shows that the components of the electromagnetic field E_y and H_x are equal; therefore, the real part of the longitudinal component of the Poynting vector defined by Equation (5) is proportional to the following value:

(12)$\mathrm{Re}\left[P_z\right]\simeq \mathrm{Re}\left[E_x^{*}{H}_y\right]-{\left|E_y\right|}^2.$

For an x-linearly polarized field, |E_y|² is, as a rule, small. Therefore, a significant inverse flow can be associated with the first term in Equation (12). Obviously, this term essentially depends on the form of the function defined by Equation (10). Next, we analyze the properties of this function to form a reverse energy flux on the optical axis.

3.1. Theoretical Analysis Based on the Stationary Phase Method

Since the wavenumber k for the optical wavelength range is large, we can apply the stationary phase method [41,42], the error of which is defined as O(k^−1/2). Using this method, Equation (9) can be calculated approximately based on the following relation:

(13)${\displaystyle \underset{a}{\overset{b}{\int }}\exp \left[\pm ip(\mathsf{\phi })\right]q(\mathsf{\phi })d\mathsf{\phi }}\approx \exp \left[\pm ip({\mathsf{\phi }}_0)\right]q({\mathsf{\phi }}_0)\sqrt{\frac{\pm i2\pi }{p^{″}({\mathsf{\phi }}_0)}},$

where φ₀ is a stationary point.

If there are several stationary points, then there will be a sum of terms defined by Equation (13). In our case, the problem is facilitated by the fact that, in all three integrals defined by Equation (11), function p(φ) is the same and is equal to the following:

(14)$p(\mathsf{\phi })=g(\mathsf{\phi })+kr\sin {\mathsf{\theta }}_0\cos (\mathsf{\phi }-\gamma ).$

Using Equation (13) with Equation (11) gives the following equations for the case of one stationary point:

(15)$\begin{array}{l}{E}_x^{*}\approx t^{*}\left({\mathsf{\theta }}_0,z\right)\exp \left[-ip({\mathsf{\phi }}_0)\right]A({\mathsf{\theta }}_0,{\mathsf{\phi }}_0)\sqrt{\frac{-i2\pi }{p^{″}({\mathsf{\phi }}_0)}},\\ E_y\approx t\left({\mathsf{\theta }}_0,z\right)\exp \left[+ip({\mathsf{\phi }}_0)\right]C({\mathsf{\theta }}_0,{\mathsf{\phi }}_0)\sqrt{\frac{+i2\pi }{p^{″}({\mathsf{\phi }}_0)}},\\ H_y\approx t\left({\mathsf{\theta }}_0,z\right)\exp \left[+ip({\mathsf{\phi }}_0)\right]B({\mathsf{\theta }}_0,{\mathsf{\phi }}_0)\sqrt{\frac{+i2\pi }{p^{″}({\mathsf{\phi }}_0)}}.\end{array}$

Then, from Equation (12), the following equation for Re[P_z] is obtained:

(16)$\mathrm{Re}\left[P_z\right]\sim {\left|t\left({\mathsf{\theta }}_0,z\right)\right|}^2\frac{2\pi }{\left|p^{″}({\mathsf{\phi }}_0)\right|}\left(A({\mathsf{\theta }}_0,{\mathsf{\phi }}_0)B({\mathsf{\theta }}_0,{\mathsf{\phi }}_0)-C^2({\mathsf{\theta }}_0,{\mathsf{\phi }}_0)\right)={\left|t\left({\mathsf{\theta }}_0\right)\right|}^2\frac{2\pi }{\left|p^{″}({\mathsf{\phi }}_0)\right|}\cos {\mathsf{\theta }}_0.$

This number is positive, meaning that, in the presence of a single stationary point (taking into account the approximation of the method), the presence of regions with a negative energy flux is unlikely.

If there are two stationary points, then the field components are expressed by equations similar to Equation (15), but there will be two terms corresponding to the first and second points. Omitting the details of the calculations, we here present the final expressions:

(17)$\begin{array}{l}{E}_x^{*}{H}_y=\frac{2\pi }{\left|p_1^{″}\right|}{A}_1{B}_1+\frac{2\pi }{\left|p_2^{″}\right|}{A}_2{B}_2+\frac{2\pi }{\sqrt{\left|p_1^{″}\right|\cdot \left|p_2^{″}\right|}}\left(A_1{B}_2+A_2{B}_1\right)M_1+\frac{2\pi i}{\sqrt{\left|p_1^{″}\right|\cdot \left|p_2^{″}\right|}}\left(-A_1{B}_2+A_2{B}_1\right)M_2,\\ {\left|E_y\right|}^2=\frac{2\pi }{\left|p_1^{″}\right|}{C}_1^2+\frac{2\pi }{\left|p_2^{″}\right|}{C}_2^2+\frac{2\pi }{\sqrt{\left|p_1^{″}\right|\cdot \left|p_2^{″}\right|}}2C_1{C}_2{M}_1,\end{array}$

where

(18)$\begin{array}{l}{M}_1=\{\begin{array}{l}+\cos (p_1-p_2),p_1^{″},p_2^{″}havethesamesign,\\ -\sin (p_1-p_2),p_1^{″}>0,p_2^{″}<0,\\ +\sin (p_1-p_2),p_1^{″}<0,p_2^{″}>0,\end{array}\\ M_2=\{\begin{array}{l}+\sin (p_1-p_2),p_1^{″},p_2^{″}havethesamesign,\\ +\cos (p_1-p_2),p_1^{″}>0,p_2^{″}<0,\\ -\cos (p_1-p_2),p_1^{″}<0,p_2^{″}>0.\end{array}\end{array}$

Then, after taking the real part (dropping the common factor |t(θ₀)²|), we get the following:

(19)$\begin{array}{l}\mathrm{Re}\left[P_z\right]\sim \mathrm{Re}\left[E_x^{*}{H}_y\right]-{\left|E_y\right|}^2=\\ \frac{2\pi }{\left|p_1^{″}\right|}\cdot \cos {\mathsf{\theta }}_0+\frac{2\pi }{\left|p_2^{″}\right|}\cdot \cos {\mathsf{\theta }}_0+\frac{2\pi }{\sqrt{\left|p_1^{″}\right|\cdot \left|p_2^{″}\right|}}\cdot \left[2\cos {\mathsf{\theta }}_0+{(1-\cos {\mathsf{\theta }}_0)}^2{\sin }^2({\mathsf{\phi }}_1-{\mathsf{\phi }}_2)\right]\cdot M_1\end{array}$

Thus, the longitudinal component of the Poynting vector depends on the ratio of the signs of the second derivatives at stationary points. It is easy to see that, if there are exactly two stationary points and the function p(φ) is smooth, then the signs will necessarily be different, since the two adjacent extrema must be of different types.

From Equation (14), we can obtain the following expression for finding stationary points:

(20)$p^{\prime }(\mathsf{\phi })=g^{\prime }(\mathsf{\phi })-kr\sin {\mathsf{\theta }}_0\sin (\mathsf{\phi }-\gamma )=0.$

When are there exactly two stationary points? Taking into account the fact that sin(x) = 0 twice on the period, we see that two roots are provided at condition krsinθ₀ > max|g′(φ)| for any function g(φ). Obviously, it is necessary to consider the region outside the optical axis (r > r0). The case for the optical axis is considered separately in Section 3.2. Note also that the stationary phase method will give a significant error if the stationary points are closely located. Therefore, in order to ensure a sufficient distance between stationary points, we impose a stronger condition krsinθ₀ > 2max|g′(φ)|.

For θ₀→π/2, the first term in Equation (20) is much lower than the second one, and we get two roots spaced apart by an angle close to 180 degrees. It is impossible to ignore the first term, since, in this case, the two roots are exactly π radians apart, which leads to the absence of negative regions (P_z ≥ 0) according to Equation (19).

3.2. Determination of Conditions for Inverse Energy Flow on the Optical Axis

If the goal is to obtain inverse energy flow on the optical axis, then the conditions for this are different. In this case, Equation (20) takes the simple form g′(φ) = 0. Thus, the g(φ) function is decisive. Figure 1 shows an example of a graph of a function with two extrema (two stationary points).

For definiteness, assume that $p_1^{″}<0$ (the first stationary point is the maximum) and $p_2^{″}>0$ (the second is the minimum). Moreover, to simplify calculations, assume that $\left|p_1^{″}\right|=\left|p_2^{″}\right|=p^{″}$. Then, from Equation (19), we obtain the following:

(21)$\mathrm{Re}\left[P_z\right]\simeq \frac{2\pi }{p^{″}}\cdot \left[2\cos {\mathsf{\theta }}_0+\left[2\cos {\mathsf{\theta }}_0+{(1-\cos {\mathsf{\theta }}_0)}^2{\sin }^2({\mathsf{\phi }}_1-{\mathsf{\phi }}_2)\right]\cdot \sin (p_1-p_2)\right].$

Equation (21) shows that the inequality sin(p₁ − p₂) < 0 must be satisfied to obtain negative values. If so, then sin²(φ₁ − φ₂) should be increasing. For fixed $p^{″}$, the minimum value of P_z is found if we have the following:

(22)$\left[\sin (p_1-p_2)=-1\right]\&\left[\sin ({\mathsf{\phi }}_1-{\mathsf{\phi }}_2)=\pm 1\right].$

In this case, it follows from Equation (21) that we get the following:

(23)$P_{z,\min }=-\frac{2\pi }{p^{″}}\cdot {(1-\cos {\mathsf{\theta }}_0)}^2.$

This value is negative for any θ₀, even in the paraxial case.

Thus, Equation (22) is the main result of this section and defines the condition for the inverse energy flow on the optical axis for a purely phase input field defined by Equation (10).

The conditions defined by Equation (22) mean that the difference between the values of function g(φ) at the maximum and minimum must equal 3π/2 + 2πN, and the distance between the positions of the extrema must equal π/2 or 3π/2. Considering that we need to make $p^{″}$ smaller, it is better to choose a difference of 3π/2 (n = 0) and a distance of 3π/2.

The above reasoning shows a certain inconsistency between the desire to enhance the inverse energy flow by decreasing the value of p′′ (see Equation (23)) inaccuracy increasing of the stationary phase method in this case in this case. Thus, we have no guarantee that the phase satisfying condition defined by Equation (22) will be optimal. Therefore, the final conclusion will be made on the basis of numerical modeling by using exact Equation (9).

3.3. Determination of Conditions for Off-Axis Inverse Energy Flow

The general reasoning given in Section 3.1, as well as previous studies [17,18], shows that off-axis inverse energy flow can be greater than on the axis. This can be theoretically justified by using Equation (20) for the full phase p(φ):

(24)$\begin{array}{l}p(\mathsf{\phi })=g(\mathsf{\phi })+kr\sin {\mathsf{\theta }}_0\cos (\mathsf{\phi }-\gamma )=g(\mathsf{\phi })+\mathsf{\mu }\cos (\mathsf{\phi }-\gamma ),\\ p^{\prime }(\mathsf{\phi })=g^{\prime }(\mathsf{\phi })-\mathsf{\mu }\sin (\mathsf{\phi }-\gamma ),\\ p^{″}(\mathsf{\phi })=g^{″}(\mathsf{\phi })-\mathsf{\mu }\cos (\mathsf{\phi }-\gamma ).\end{array}$

As follows from Equation (24), for μ ≠ 0 (i.e., outside the optical axis), at least at one stationary point, $p^{″}$, becomes less than for μ = 0. In general, it is impossible to obtain specific results, so we consider some examples below.

3.4. Analysis of Special Cases and Simulation Results

This section discusses various specific examples and special cases of GSPP. In special cases, it is much easier to carry out theoretical analysis, so in this section, we consider not only the longitudinal but also the transverse components of the Poynting vector, as well as account for not only the real but also the imaginary part of the Poynting vector, which contributes to the transverse force acting on optically trapped Mie particles [29,30]. In the simulation, the following calculation parameters were used: the annular aperture had a central radius θ₀ = 80° and width Δ = 20°.

3.4.1. Example 1: Classic SPP

First, consider the case of a standard SPP with a phase:

(25)$g(\mathsf{\phi })=m\mathsf{\phi },$

where m is an arbitrary real number.

For integer m, the calculations in Equations (11) and (12) are performed analytically and give the following result:

(26)$\begin{array}{l}{P}_z=\mathrm{Re}\left[P_z\right]+i\mathrm{Im}\left[P_z\right]=\\ =\left[{\pi }^2{(1+\cos {\mathsf{\theta }}_0)}^2{J}_m^2(\rho )-({\pi }^2/2){(1-\cos {\mathsf{\theta }}_0)}^2(J_{m+2}^2(\rho )+J_{m-2}^2(\rho ))\right]-\\ -i\left[{\pi }^2{\sin }^2{\mathsf{\theta }}_0\sin 2\gamma J_m(\rho )\left(J_{m+2}(\rho )-J_{m-2}(\rho )\right)\right],\end{array}$

(27)$\begin{array}{l}{P}_r=i\mathrm{Im}\left[P_r\right]=i{\pi }^2\sin {\mathsf{\theta }}_0\cos 2\gamma [(1+\cos {\mathsf{\theta }}_0)J_m(\rho )\left(J_{m-1}(\rho )-J_{m+1}(\rho )\right)+\\ +0.5(1-\cos {\mathsf{\theta }}_0)\left(J_{m-1}(\rho )J_{m+2}(\rho )-J_{m+1}(\rho )J_{m-2}(\rho )\right)],\end{array}$

(28)$\begin{array}{l}{P}_{\gamma }=\mathrm{Re}\left[P_{\gamma }\right]+i\mathrm{Im}\left[P_{\gamma }\right]={\pi }^2\sin {\mathsf{\theta }}_0\{(1+\cos {\mathsf{\theta }}_0)J_m(\rho )\left[J_{m-1}(\rho )+J_{m+1}(\rho )\right]+\\ +(1-\cos {\mathsf{\theta }}_0)\left[J_{m-2}(\rho )J_{m-1}(\rho )+J_{m+1}(\rho )J_{m+2}(\rho )\right]\}+\\ +i{\pi }^2\sin {\mathsf{\theta }}_0\sin 2\gamma \{(1-\cos {\mathsf{\theta }}_0)\left[J_{m-1}(\rho )J_{m+2}(\rho )-J_{m+1}(\rho )J_{m-2}(\rho )\right]-\\ -(1+\cos {\mathsf{\theta }}_0)J_m(\rho )\left[J_{m-1}(\rho )-J_{m+1}(\rho )\right]\},\end{array}$

where $\rho$ = krsinθ₀.

As seen from Equation (26), the real part of the longitudinal component of the Poynting vector Re[P_z] has radial symmetry and does not change when m changes sign. When the sign of m changes, the imaginary part Im[P_z] changes its sign to the opposite and, for any m, has rotational symmetry of the second order due to the presence of an angular dependence sin(2γ). A negative value for Re[P_z] at the origin occurs only at m = ± 2. For m = 0 and m = ± 1, the results have a simpler form. In particular, for m = 0 $\mathrm{Re}\left[P_z\right]={\pi }^2{(1+\cos {\mathsf{\theta }}_0)}^2{J}_0^2(\rho )-{\pi }^2{(1-\cos {\mathsf{\theta }}_0)}^2{J}_2^2(\rho )$, Im[P_z] = 0, and for m = ±1 $\mathrm{Re}\left[P_z\right]={\pi }^2{(1+\cos {\mathsf{\theta }}_0)}^2{J}_1^2(\rho )-({\pi }^2/2){(1-\cos {\mathsf{\theta }}_0)}^2\left[J_1^2(\rho )+J_3^2(\rho )\right]$, $\mathrm{Im}\left[P_z\right]=\mp {\pi }^2{\sin }^2{\mathsf{\theta }}_0\sin 2\gamma \cdot J_1(\rho )\left[J_1(\rho )+J_3(\rho )\right]$.

The radial component, as follows from Equation (27), is purely imaginary. Equation (28) shows that the real part of the azimuthal component Re[P_γ] is radially symmetric and equal to zero at m = 0.

It is more difficult to obtain analytical expressions for fractional values of m (as a rule, the representation through the sum for integer m is used), and they will be less descriptive. Therefore, we consider only the simulation results for these fractional values. Figure 2 shows the simulation results for standard integer and fractional SPPs with a transmission function of exp(imφ).

The distributions of the intensity of the electric |E|² and magnetic |H|² components of the field in the focal plane (z = 0), the distributions of individual components of the Poynting vector, the real and imaginary parts of the longitudinal component, and the real part of the angular component are shown. The longitudinal component is the main focus of this study, so we consider it comprehensively. The angular component of the Poynting vector is closely related to the angular momentum of the field, which is often formed by using SPP.

The simulation results for integer m fully correspond to analytical Equations (26) and (28). In particular, the real part of P_z has radial symmetry, and the imaginary part P_z has a structure corresponding to sin(2γ). Fractional orders of m lead to distortion of the mentioned symmetries. However, for Im[P_z], there is complete symmetry of negative and positive values for any m, which leads to a zero mean value of the imaginary part. Nevertheless, locally, it can affect optically trapped Mie particles [29,30] and for m = ± 1, the local maxima have the greatest value.

The maximum relative negative value for Re[P_z] is observed when m = 2 and it is located at the origin. Note that, for integer m, the distribution of the real part P_z is structurally similar to that for radial polarization of the mth order [17]. In both cases, the generated distributions are radially symmetric and there is a negative value on the axis only for m = 2. However, there are some differences: for radial polarization, the distribution changes slightly when the sign of m changes; when m = 1, there are no negative values. For the imaginary part, the difference is fundamental: for radial polarization, Im[P_z] = 0, and for linear polarization for classic SPPs, Im[P_z] ≠ 0 for any m ≠ 0 (see Figure 2). Obviously, for negative values of m, the results are similar; however, Im[P_z] and Re[P_γ] change signs (compare the lines in Figure 2 corresponding to m = 2 and m = −2).

In Reference [18], superposition 0.5[exp(imφ) + exp(–imφ)] = cos(mφ) was investigated. The case m = ± 2 was also noted as special in terms of maximum negative values Re[P_z] on the optical axis. However, References [17,18] also pointed out the meaning of using an illuminating beam cos(mφ) at |m| > 2, since this allows us to increase the total area and energy of negative regions. The use of fractional values of m allows for continuous changes in the properties of these areas, including changing their shape. Below, Example 2 shows that, for sin(mφ), there is a similar result up to rotation.

3.4.2. Example 2: Amplitude Spiral Plate with a Phase Shift

Consider a spiral plate defined as follows:

(29)$F_0(\mathsf{\phi })=\sin \left(m\mathsf{\phi }+{\mathsf{\phi }}_0\right).$

Add to it a phase corresponding to the compensation of the kernel phase defined by Equation (1):

(30)$g(\mathsf{\phi })=-\alpha \cos (\mathsf{\phi }),$

where α is an arbitrary real number.

Then, p(φ) = –αcos(φ) + μcos(φ − γ) equals zero (for γ = 0) on a ring with a radius of r₀ = α/(ksinθ₀). It is similar to phase compensation in the focal plane of the lens. In this case, (γ = 0, α = μ) and $F_0(\mathsf{\phi })$ = 1, and all quantities included in Equation (12) are calculated exactly:

(31)$\begin{array}{l}{E}_x^{*}\sim {\displaystyle \underset{0}{\overset{2\pi }{\int }}A({\mathsf{\theta }}_0,\mathsf{\phi })}d\mathsf{\phi }=\pi (1+\cos {\mathsf{\theta }}_0),\\ E_y^{*}\sim {\displaystyle \underset{0}{\overset{2\pi }{\int }}C({\mathsf{\theta }}_0,\mathsf{\phi })}d\mathsf{\phi }=0,\\ H_y\sim {\displaystyle \underset{0}{\overset{2\pi }{\int }}B({\mathsf{\theta }}_0,\mathsf{\phi })}d\mathsf{\phi }=\pi (1+\cos {\mathsf{\theta }}_0),\\ P_z\simeq E_x^{*}{H}_y-{\left|E_y\right|}^2={\pi }^2(1+\cos {\mathsf{\theta }}_0{)}^2>0.\end{array}$

Thus, in the considered region, the energy flux is positive.

However, the situation can be changed by using a non-unit amplitude, as was used in Reference [18]. Consider an amplitude that co-phases to C(θ₀,φ), that is, equals sin(φ)cos(φ) = 0.5 sin(2φ). This is the light field defined from Equation (29), when m = 2 and φ₀ = 0. In this case we have the following:

(32)$\begin{array}{l}{E}_x^{*}\sim {\displaystyle \underset{0}{\overset{2\pi }{\int }}\sin \mathsf{\phi }\cos \mathsf{\phi }A({\mathsf{\theta }}_0,\mathsf{\phi })}d\mathsf{\phi }=0,\\ E_y^{*}\sim {\displaystyle \underset{0}{\overset{2\pi }{\int }}\sin \mathsf{\phi }\cos \mathsf{\phi }C({\mathsf{\theta }}_0,\mathsf{\phi })}d\mathsf{\phi }=-\frac{\pi }{4}(1-\cos {\mathsf{\theta }}_0),\\ H_y\sim {\displaystyle \underset{0}{\overset{2\pi }{\int }}\sin \mathsf{\phi }\cos \mathsf{\phi }B({\mathsf{\theta }}_0,\mathsf{\phi })}d\mathsf{\phi }=0,\\ P_z\simeq E_x^{*}{H}_y-{\left|E_y\right|}^2=-\frac{{\pi }^2}{16}(1-\cos {\mathsf{\theta }}_0{)}^2<0.\end{array}$

Thus, in this case, the energy flow is inverse. Consider the amplitude defined by Equation (29) and the phase defined by Equation (30) together:

(33)$F_0(\mathsf{\phi })=\sin \left(2\mathsf{\phi }\right)\exp \left[-i\alpha \cos (\mathsf{\phi })\right].$

For a field defined by Equation (33), using the tabular formulas Section 2.5.27 in Reference [43]) gives the following in the explicit form:

(34)$\begin{array}{l}{E}_x^{*}=0,\\ H_y=0,\\ E_y^{*}=-2\pi (1-\cos {\mathsf{\theta }}_0)\left[J_0(\alpha -\mathsf{\mu })-J_4(\alpha -\mathsf{\mu })\right].\end{array}$

Equation (34) shows that, for γ = 0, the energy flux is negative and equal:

(35)$P_z\simeq -a{\pi }^2{(1-\cos {\mathsf{\theta }}_0)}^2{\left[J_0(\alpha -kr\sin {\mathsf{\theta }}_0)-J_4(\alpha -kr\sin {\mathsf{\theta }}_0)\right]}^2.$

Note that the imaginary part Im[P_z] in this case is equal to zero. Substituting α = 0 (there is no phase) gives an expression that coincides with the solution obtained in Reference [18] for the amplitude cos(2φ). Then an inverse flux was observed on rays corresponding to γ = ±π/4. At α ≠ 0, the phase g(φ) = –αcos(φ) plays a role similar to the linear phase with a spatial carrier (analogue of a prism): it shifts the center of the distribution from the radius r = 0 to the radius r₀ = α/(ksinθ₀). Figure 3 shows the corresponding simulation results. An analytical calculation for integer m shows that P_z is real, but P_γ is purely imaginary.

The results for integer m are similar to those obtained in Reference [18] for the amplitude cos(mφ) up to rotation. However, the introduction of an additional phase was not considered in Reference [18]. Such a phase does not change the main characteristics of the Poynting vector, but allows the center of the distribution to be shifted from the optical axis in accordance with Equation (35). This is important because the maximum negative value has so far been observed only at the origin of coordinates [17,18,19].

It can be seen that Im[P_z] and Re[P_γ] are nonzero only for fractional values of m in the amplitude; however, the average value also remains zero. Note also that, in examples with a half-integer frequency, the distribution Re[P_z] has two axes of symmetry and two positive maxima. Perhaps this is because half-integer frequencies retain some properties of integers, in particular, equality sin(m2π) = 0. Therefore, some features may appear, as was observed in parabolic beams [44].

3.4.3. Example 3: Power-Exponent Phase Plate

A PEPP [40,41] has the following phase profile:

(36)$g(\mathsf{\phi })=m{\mathsf{\phi }}^s,$

where m is an arbitrary real number and s is a positive real number.

The light field formed by the plate with phase defined by Equation (36) can be represented by an infinite sum of fields formed by classic integer-order SPPs. Therefore, analytical formulas do not give a visual representation of the properties of the generated field, but References [9,10] showed that the line of maximum intensity has the shape of a spiral. In the paraxial mode, Re[P_z] is proportional to the field intensity.

Figure 4 shows the modeling results for PEPP defined by Equation (36).

Our analysis of the results shown in Figure 4 reveals that, as in Example 1 (classic SPP), the lines of the level of the real parts of P_z and P_γ are approximately parallel to the mentioned spiral: for the imaginary part of P_z, this is expressed less clearly, since sector dependence is imposed, although not as explicitly as in Example 1. Moreover, Im[P_z] no longer has symmetry of negative and positive values; however, the zero mean value of the imaginary part may remain.

3.4.4. Example 4: Light Fields Resulting from the Stationary Phase Method

In Section 3.1, we considered the problem of the formation of an inverse flow on the optical axis within the framework of the stationary phase method. The conditions for maximizing the negative value of Re[P_z] on the optical axis, imposed on the phase function, were formulated at the end of Section 3.2. These conditions, for example, are satisfied by a function of the following form (see Figure 1):

(37)$g_1(\mathsf{\phi })=\frac{3\pi }{4}\sin \left[\frac{2}{3}\left(\mathsf{\phi }+\frac{\pi }{2}\right)\right].$

Figure 5 shows the simulation results for the input field, which is a narrow ring with a phase function g₁(φ) from Equation (37).

The two functions defined by the following polynomials have similar properties:

(38)$g_2(\mathsf{\phi })=\frac{2}{27}{(\mathsf{\phi }-\pi )}^3-\frac{{\pi }^2}{8}\left(\mathsf{\phi }-\pi \right).$

(39)$g_3(\mathsf{\phi })=\frac{8}{9{\pi }^2}{\left(\mathsf{\phi }-\pi \right)}^3-\frac{3}{2}\left(\mathsf{\phi }-\pi \right).$

All three functions have a maximum value at a point φ₁ = π/4 and a minimum value at a point φ₂ = 7π/4. The first and second functions have the same values of the second derivatives at these points, but the difference, g_max − g_min = π³/8, is less than the value of 3π/2 mentioned in Section 3.2, which is in the first function. The third function has a difference of 3π/2, but the second derivatives at the extrema are slightly larger in absolute value than those of the first and second functions. Figure 6 shows graphs of all three functions. The figured shows they are very close and, therefore, will give similar results.

As follows from the calculation results shown in Figure 5, predictions based on the stationary phase method did not give the expected effect: negative values of the energy flux are small and, moreover, are mainly observed off-axis. Obviously, this is due to the approximation of the stationary phase method and, possibly, not the full legitimacy of its use.

Significant changes can be expected if the amplitude is made to equal zero outside the vicinity of the extremum points, but this entails a loss of the incident beam energy. Therefore, it is more efficient to make the second derivative equal zero (although the stationary phase method is no longer directly applicable here), i.e., to use the binary version of the considered phase functions, as shown in the following example.

3.4.5. Example 5: Binary Phase

Reference [18] showed that the distribution of the Poynting vector changes little if the incident beam with amplitude cos(mφ) is replaced by a beam with unit amplitude and binary phase.

Therefore, consider a phase function of the following form:

(40)$g_b(\mathsf{\phi })=\pi \mathrm{sgn}\left[\sin \left(m\mathsf{\phi }+{\mathsf{\phi }}_0\right)\right]/2+h.$

where $\mathrm{sgn}(\cdot )$ is the sign function and h is a real constant. In contrast to Reference [18], m can be a non-integer. Obviously, parameters φ₀ and h do not qualitatively affect the distribution pattern: φ₀ determines its rotation, and h should not change it at all, but h ≠ 0 may be more convenient when carrying out an experiment. The simulation results for binary phase functions defined by Equation (40) for various values of m are shown in Figure 7.

Comparing Figure 3 and Figure 7 reveals that the distributions in the focal plane for the binary phases strongly correlate with distributions for the amplitude in the form of a sine, as it is actually supposed to be. There are two main differences: the spots have slightly changed their shape and the energy efficiency is increased.

4. Discussion

The obtained analytical and numerical results allow us to draw some conclusions about the preferred variants of the input field for amplifying the reverse energy flow in different areas of the focal plane.

Obviously, to ensure the reverse energy flow, Equation (5) must take negative values in a certain area. However, if the specific form of the fields participating in Equation (5) is unknown, it is impossible to make any assumptions. Note that Equation (5) includes the field distributions in the focal region; therefore, it is difficult to guess what the initial field (in the input plane) should be to achieve this goal. To do this, it would be necessary to solve the inverse problem for vector fields taking into account Maxwell’s equations, which is difficult in general form. Therefore, we have limited the consideration to a certain class of input fields. In particular, the field was assumed to be linearly polarized. This made it possible to reduce the general Equation (5) to a particular Equation (12), which is more illustrative for the search for possible solutions. Moreover, to perform analytical calculations, we considered the input field in the form defined by Equation (7), imposing a narrow annular aperture, and thus obtained a dependence on only one variable.

The imposition of such restrictions allowed us to apply the stationary phase method to calculate the integrals and obtain condition defined by Equation (22) for inverse energy flow on the optical axis for a purely phase input field (the amplitude is constant). However, due to the approximation of the stationary phase method, which is especially manifested on the optical axis, the obtained analytical functions for the input field did not give a large value of the inverse energy flow. However, the analysis based on the stationary phase method suggested that the inverse energy flow out of the optical axis may be greater than at the axis.

Due to the limitations of the stationary phase method, we proceeded to the analysis of exact expressions for several types of the input field. Among the considered fields, two types are special: (1) the classical spiral phase plate of the second order g(φ) = 2φ, which provides the maximum inverse energy flow on the optical axis; and (2) a generalized spiral phase plate g(φ) = −αcos(φ + φ₁) with the amplitude sin[2(φ + φ₂)], which forms regions of the inverse energy flow with a larger size and allows this region to be displaced.

5. Conclusions

Structured laser beams with the predetermined amplitude, phase, and polarization distributions play an important role in the field of laser manipulation. The possibility to form light fields with the desired distribution of the Poynting vector representing the directional energy flux of an electromagnetic field allows one to affect optically trapped particles in a controlled manner. In this paper, we demonstrated the possibility to form light fields with an inverse energy flux and control the distributions of the real and imaginary parts of the Poynting vector, using conventional and GSPPs illuminated with linearly polarized laser radiation. Our theoretical analysis based on the stationary phase method provides the conditions for generating inverse energy flux on the optical axis in the case of conventional and generalized SPPs. We also obtained the conditions for maximizing the negative value of the real part of the longitudinal component of the Poynting vector.

We showed that, for a classic SPP exp(imφ), the maximum relative negative value for the real part of the longitudinal component of the Poynting vector is observed when m = 2 and it is located at the origin. In addition, for integer m, the real part of the longitudinal component has radial symmetry, and the imaginary part of the longitudinal component has a structure corresponding to sin(2φ). The use of a special phase allows us to shift the center of the Poynting vector distribution from the optical axis, so the maximum negative value of the inverse energy flux will be out of the origin of coordinates. The use of PEPPs allows us to break the symmetry of negative and positive values of the imaginary part of the longitudinal component of the Poynting vector; however, the zero mean value of the imaginary part may remain. Currently, laser sources with a linearly polarized output radiation are the most readily available, which makes the results obtained in this work attractive for use in the field of laser manipulation.

Author Contributions

Conceptualization, S.N.K.; methodology, S.N.K. and A.V.U.; software, S.N.K.; validation, S.N.K. and A.P.P.; formal analysis, S.N.K. and A.V.U.; investigation, S.N.K. and A.P.P.; writing—original draft preparation, S.N.K. and A.P.P.; writing—review and editing, S.N.K. and A.P.P.; supervision, S.N.K.; funding acquisition, A.P.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Russian Foundation for Basic Research (grant number 20-37-70025) and Russian Federation Ministry of Science and Higher Education (007-GZ/Ch3363/26).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare that there are no conflict of interest related to this article.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Function g1(φ): original function (blue curve) and function reduced to positive values (red curve).

Figure 2. Simulation results for classic SPPs exp(imφ) of various orders m (red color—for x-components, green color—for y-components, and blue color—for z-components). The distributions in the focal plane with a size of 4λ × 4λ are shown.

Figure 3. Simulation results for amplitude spiral plates (red color—for x-components, green color—for y-components, and blue color—for z-components). The distributions in the focal plane with a size of 4λ × 4λ are shown.

Figure 4. Simulation results for PEPPs exp(imϕs) (red color—for x-components, green color—for y-components, and blue color—for z-components). The distributions in the focal plane with a size of 6λ × 6λ are shown.

Figure 5. Simulation results for the input field with phase g1(φ) defined from Equation (37): amplitude (a) and phase (b) distribution of the input field; |E|2 (c) and |H|2 (d) distributions in the focal plane, as well as the distributions of Re[Pz] (e), Im[Pz] (f), and Re[Pγ] (g).

Figure 6. Phase functions: g1(φ) defined by Equation (37) (blue line), g2(φ) defined by Equation (38) (green line), and g3(φ) defined by Equation (39) (red line).

Figure 7. Simulation results for binary phases (red color—for x-components, green color—for y-components, and blue color—for z-components). The distributions in the focal plane with a size of 4λ × 4λ are shown.