1. Introduction

Optical beams possessing spatial variations in amplitude, phase, and/or polarization have supported an increasing diversity of applications in recent years [1]. As the technology required to produce this “structured light” has improved (e.g., spatial light modulators), researchers have explored optical manipulation of particles [2,3], improved microscopy [4], and developed new communications [5] and cryptography schemes [6] via spatial control of beam properties. Of particular relevance to this work are beams possessing tailored spatial variation in their state of polarization, i.e., “vector beams”. Such beams are said to be “classically entangled” in the sense that their spatial and polarization distributions are non-separable [7]. Various methods for producing such beams have been developed [8], [9], [10] and subsequently applied to microscopy [11], spectroscopy [12], communications [13], and even material processing [14].

In this work we predict a new phenomenon associated with vector beams, namely a curved trajectory that results from spatial variations in the state of linear polarization. We show that such beams will “self-bend” upon exiting an apparatus without requiring any variations in the properties of the medium through which the beam travels.

To arrive at our conclusion, we view beam propagation as a transport problem. Under this model, the complex electric field amplitude is governed by the Transport of Intensity Equation (TIE) [15,16], expressing conservation of optical intensity, while the phase gradient is seen to obey an expression resembling the momentum balance commonly used in mechanics [17]. It is this transport formulation of the eikonal equation that is the focus of this work. To date, optical transport models have been applied only to scalar propagation problems which cannot, by definition, account for spatial variations in polarization. However by extending the transport model to include such variations, we show there appears a new term comprising spatial polarization derivatives that is seen to produce a transverse change (curvature) in the beam path, even in free space. This term clearly suggests a new mechanism for controlled beam bending as we demonstrate both theoretically and experimentally.

We begin with the familiar Helmholtz equation for propagation of monochromatic waves, assume a complex vector solution in the form of a linearly-polarized wave, and allow the state of linear polarization to depend on transverse position in the beam. After invoking the paraxial assumption and making the appropriate definition of terms, we arrive at a system of three partial differential equations that comprise the transport model. The problem is then solved by leveraging the Lagrangian viewpoint, known in optics as the method of characteristics [18,19]. By switching to Lagrangian coordinates, we are able to solve analytically for the polarization-induced, curved path. We consider both the diffraction-free case (bending only) and the case with diffraction included (bending and spreading of the beam).

Lastly, we verify these predictions in experiment. We use a spatial light modulator (SLM) to produce the needed polarization distribution and then record the relative shift in beam position as a function of propagation distance out to 34 meters. Here, relative shift refers to the transverse distance between the two beam profiles - one with uniform polarization distribution (SLM off) and one with quadratic polarization distribution (SLM on). The effect is not large, amounting to a shift of approximately 4mm at 30 meters (similar in magnitude to that observed in Airy beam bending), but it may have serious implications for beams travelling over much larger terrestrial or interstellar distances.

We emphasize that the physics described here differ fundamentally from Airy beam bending, where the intensity distribution shifts during propagation but where the center of mass remains unchanged [20]. In fact, we will show that the Airy beam bending is predicted by an entirely different term in the transport model, specifically the term governing standard diffraction. The Airy bending effect has been leveraged in optical particle manipulation [21] and microscopy [22], for example. Additional work by Zhou et al. altered the spatial distribution of Airy beam phase to produce Airy vortex beams [23] while more recent work by Zhao et al. produced Airy beams that carry spatially distributed polarization states along the path of the main Airy beam lobe [7]. By combining both the parabolic trajectories and the spatial phase and polarization variations, one is afforded additional degrees of freedom in applications such as optical trapping and optical tweezers [1].

In contrast, the parabolic path we predict, and subsequently confirm, is entirely a consequence of the spatial polarization variations and not the intensity profile. In fact, the beam used in this work possesses a standard Gaussian intensity distribution. Moreover, unlike the Airy beams, it is the centroid of our vector beam that accelerates rather than the center of intensity of a particular lobe. This new mechanism for creating and controlling the path of an optical field may offer improvements in the aforementioned application areas and may also open new areas of inquiry.

The analysis and discussion of this new effect draws from a number of diverse fields including standard electromagnetic theory, transport of intensity approaches, Lagrangian dynamics and coordinates, Stokes analysis of polarization transformations, and Pancharatnam-Berry phase. Therefore, we have included a set of Appendices offering brief overviews of various topics of relevance.

2. Derivation of the eikonal transport equations

2.1 Background

We begin the analysis in the general case of an isotropic, weakly inhomogeneous medium and, later in the calculation, show that the effect is present even in free space. For a monochromatic beam propagating in the z-direction, the complex electric field amplitude vector is governed by the Helmholtz equation

Substituting (2) into (1) and equating the imaginary parts of the result yields

The only way for Eq. (3) to be satisfied is if (a) $\partial \rho /\partial z+\nabla _X\cdot (\rho \vec {v})=0$ and (b) $D\gamma /Dz=\partial \gamma /\partial z+\nabla _X\gamma \cdot \vec {v}=0$ where the notation $D(\cdot )/Dz$ denotes the total derivative (see Appendix 6). The expression (a) is the TIE equation alluded to previously (see again, [15], [16]) while (b) states that, within the Lagrangian coordinate system (see section 2.3), the polarization angle distribution across the face of the beam does not change with propagation distance [25]. We will leverage this property later in the derivation.

2.2 Eikonal equation

Two mechanisms will be shown to contribute to $\vec {v}$. The first is diffraction, elaborated on in section (3.2). The second is the main topic of this paper, namely, a second-order spatial distribution in state of polarization across the transverse beam profile which derives naturally from the beam eikonal. Gathering the real portions of (1), defining $\eta =(n^2-1)/2$, and again invoking the paraxial assumption requires that

Taking the transverse gradient of the eikonal, and defining $\vec {\omega }=k_0^{-1}\nabla _X\gamma$, we arrive at the system of equations

Importantly, however, by considering the vector electric field model (2) and allowing a spatially varying polarization gradient we observe an additional “potential” term in the eikonal governing the beam path (6c) namely, $(\vec {\omega }\cdot \nabla _X)\vec {\omega }$. This term can produce a freely accelerating beam independent of the intensity distribution. Specifically, this model predicts that through proper choice of spatial polarization distribution, it is possible to cause the beam to follow a curved trajectory, even in the absence of variations in the properties of the medium (that is, even with $\nabla _X\eta =0$). Furthermore, because the polarization angle distribution is independent of $z$ (Eq. 6b) the effect will persist over the the entire propagation distance and, by Eq. (6a), will preserve the total beam intensity. To conclude this section, we point out that one can also arrive at the expression governing optical path (6c) via the principle of least-action. We will briefly explore this approach in section (4).

2.3 Effect of a polarization gradient

Consider first the case where $\gamma$ varies only in $\hat {y}$ and where there are no index fluctuations or intensity-dependent (diffraction) effects. While the first component of the vector Eq. (6c) states that the beam path will remain unchanged in $\hat {x}$, the second component becomes

This is a fundamental result of the paper. It shows that, provided $d\gamma /dy\neq 0$ and $d^2\gamma /dy^2\neq 0$, then $Dv/Dz\neq 0$, that is, the beam centroid "accelerates" in the $\hat {y}$-direction as it propagates in the $z-$direction, even in free space. The non-zero second-order gradient in $\gamma$, imposed by the experimental set-up, guarantees both conditions. Furthermore, although the velocity term $v$ in (4) would vanish identically in homogeneous free space where $\phi (x,y,z)$ must be constant, in the case considered here the beam itself acquires a spatially-dependent geometric phase $\phi (x,y,z)$ (Panchratnam-Berry phase) during the polarization preparation process (see Section 4.3).

To solve the eikonal (5), we change to Lagrangian coordinates, e.g., $x\rightarrow x(z,x_0,y_0)\equiv x_z$, whereby the transverse locations are functions of propagation distance and the starting location (see Appendix 6). Both initial and Lagrangian locations are continuous quantities, defined everywhere in the optical field they are being used to describe. With this transformation, the total derivative becomes an ordinary derivative and the “velocities” become $u_z=dx_z/dz,~v_z=dy_z/dz$ (see, for example, [24]). Equation (6b) then becomes $d\gamma (y_z)/dz=0$ meaning that $\gamma (y_z)=\gamma (y_0)$ is a constant function of the initial (Eulerian) coordinates. Thus, Eq. (7) governing the optical beam path becomes

Finally, note that under the conditions specified here, $|y_0|\leq a/2$ and, hence, (8) in Lagrangian coordinates becomes

The model governing the propagation path of an optical beam with quadratically-distributed linear polarization is thus mathematically equivalent to that of a projectile travelling in the $+\hat {z}$-direction while in the presence of a constant gravitational field acting in $\hat {y}$. In both cases (optics and mechanics) the observed positions are appropriately modeled in a curvilinear, Lagrangian coordinate system defined by the problem physics that is, by the“least action” geodesics [33] (see section 4. for the action associated with the model described herein). Motion in a gravitational field can be described geometrically by such geodesics in space-time. Likewise, in our experiment the SLM array creates an optical field in which the light from each pixel simply propagates along its corresponding geodesic curve (See Fig. 7 in Appendix 6).

When viewed in the rectilinear Cartesian system of the laboratory, however, the beam appears to propagate as if it were acted upon by a constant force. Just as in the case of a particle moving in a gravitational field, the “bending” effect observed in our experiment is thus properly described as a geometric effect, without the need to invoke the concept of forces. In fact, the beam follows Lagrangian geodesics precisely because no forces act on it. Shortest-path lines (geodesics) in the curvilinear coordinate system are not shortest path lines (straight lines) in the Cartesian system.

In a sense, the effect demonstrated here is reminiscent of the famous measurement of 1919 by Dyson, Eddington and Davidson [34] that provided support to the theory of general relativity (GR). In that measurement, light from a distant star followed a curved trajectory in passing near the Sun. When viewed within the context of GR, however, the light simply propagated along a geodesic in the otherwise flat spacetime coordinate system that had become locally curved by the mass of the Sun. In the bending effect described here, each portion of the beam also follows a geodesic path in a curvilinear Lagrangian coordinate system, and where the polarization gradient in the optical field produced by the SLM induces the curvature (in mathematical analogy to the curvature induced by mass in GR).

3. Experiment

3.1 Overview

To experimentally verify Eq. (11) we generated a beam with the desired polarization angle profile by using a reflective liquid crystal on silicon SLM (HoloEye-Pluto Telco) and the optical setup shown in Fig. 1(a). Specifically, a linearly polarized beam (prepared using the first set of waveplates and polarizer) impinges first upon the phase mask located on the right half of the SLM, which applies a spatially-dependent phase retardance to rotate the beam’s polarization pixel-by-pixel. After the first SLM incidence, a uniform polarization rotation is applied using another half-wave plate (HWP) with fast-axis oriented at 247.5 deg with respect to the laboratory x-axis. Then, the beam impinges upon a second phase mask located on the left half of the SLM, applying another spatially-dependent polarization rotation. A map of the pixel locations on the SLM at which the two passes occur is shown in Fig. 1(b). This two-pass SLM scheme, combined with the intervening HWP, allows for the generation of nearly arbitrary spatial polarization distributions (See Section 3.3 and Appendix 9). The required retardances for each phase mask were calculated using the known input and targeted output polarization distribution (10) and the theory described in [35] and [27]. A Stokes polarimeter, consisting of rotating QWP, polarizer, and an infrared (IR) camera (Xenics Bobcat 640E), was used to measure the 2D polarization distribution across the transverse profile of the beam.

 

Fig. 1. (a) Diagram of the optical setup adapted from [35] used to generate the required $\gamma (y_0)$. (b) A map of the pixel locations on the SLM at which the first and second pass occur as viewed looking onto the face of the SLM and into the beam. (c) Targeted (red line) and measured $\gamma (y_0)$ across the center of the beam’s transverse profile (sliced at $x=0~mm$) (black points). The inset shows the full measured 2-D polarization angle distribution of the beam (arrows) overlaid upon its normalized intensity profile. M: mirror, QWP: quarter-wave plate, HWP: half-wave plate, P: linear polarizer, NPBS: non-polarizing beam-splitter, l: lens.

Download Full Size | PDF

Figure 1(c) shows the measured polarization angle distribution across the center of the beam (sliced at $x=0~mm$) after propagating $z=5~cm$ and the targeted polarization angle distribution given by (10) with $a=4.17~mm$. The beam’s normalized Gaussian intensity profile and the 2D polarization angle distribution are shown in the inset. The beam was almost completely linearly polarized with an eccentricity of $>0.98$ across its entire transverse profile.

The transverse displacement of the beam center was then measured as a function of distance using the IR camera mounted on a wheeled cart. Values of the transverse displacement were obtained based on the known physical dimensions of the IR camera array. In order to propagate the beam beyond a few meters without using a multitude of mirrors, the beam was directed out of the optical laboratory and down the inner hallway of the building. To reduce the effects of air turbulence, 10 separate 2D transverse profiles of the beam were captured at each $0.5m$ increment along the propagation path (10 each of the beam with phase masks on and off for a total of 1200+ profiles). Each profile was fit to a Gaussian beam profile and the peak value recorded.

Figure 2(a) shows the beam’s transverse profiles (and fits) at several example propagation distances while Fig. 2(b) shows the mean transverse displacement of the beam’s maximum intensity peak along with the intervals spanning the 10 observations associated with each distance. With the SLM turned on (polarization gradient applied), the measured transverse displacements agree extremely well with the curve predicted using (11) with $y_0=0~mm$ and $a=4.17~mm$. Also shown is the path taken by the beam when the polarization gradient is removed (SLM phase masks off). In this case, the beam propagates without bending as predicted. The slight deviations of the measured transverse displacements from theory at larger propagation distances (>15m) can be attributed to beam pointing errors due to air turbulence, e.g. from the opening and closing of doors down the hallway, as well as slight shifts in the position of the IR camera and cart in-between measurements of the beam with phase masks on and off.

 

Fig. 2. (a) Example transverse profiles of the beam at several propagation distances, $z$, and (b) The mean transverse displacement of the beam’s maximum intensity peak (black data points) and the curve predicted using Eq. (11) with $y_0$=0 and $a=4.17~mm$ (red dashed line) plotted as a function of distance. Note that when the polarization gradient is removed (phase masks “off”) the beam’s path is unaltered.

Download Full Size | PDF

3.2 Inclusion of diffraction

The influence of diffraction on the beam profile is clearly observed in the profiles of Fig. 2(a). Including this effect in the analytical model requires solving

Figure 3(a) shows the normalized intensity profiles of the experimentally generated beam plotted as a function of propagation distance as well as the beam profile obtained via Eqs. (13). The theory thus correctly captures both the beam bending as well as beam spreading due to diffraction. Diffraction is therefore seen to be largely decoupled from the polarization-induced curvature in this case, and it’s effects can be accurately modeled using standard Gaussian beam propagation models.

 

Fig. 3. (a) Normalized y-axis intensity slices taken from single images of the beam vs. propagation distance. Also plotted for reference are the paths predicted using Eq. (11) for $y_0=0,\pm 750\mu m,\pm 1500\mu m$ and $a=4.17~mm$ (white dashed lines). (b) The corresponding normalized beam profile predicted by theory (Eq. (13a)) using the measured initial beam waist $\sigma _0=0.0013~m$.

Download Full Size | PDF

3.3 Description of the polarization transformations

In this section we describe the 3-step polarization transformation executed by the experimental apparatus. The analysis is performed in Stokes space and illustrates how the required spatial distribution in linear polarization state is realized experimentally (for additional details, see Appendix 9). An arbitrary state of polarization corresponding to a point on the unit Poincare sphere is represented by unit-length vector $\hat {V}$. Each polarization transformation is assumed to be lossless and is represented by an orthogonal matrix $R_{\hat {\beta }}(\theta )$ producing rotation through angle $\theta$ on the Poincare sphere about an axis defined by unit vector $\hat {\beta }= [\beta _1 \,\,\,\,\, \beta _2 \,\,\,\,\, \beta _3]^T$. Throughout, we use the standard Euler-Rodrigues form of rotation matrix (see Appendix 8. ).

In the apparatus shown in Fig. 1(a), the combination QWP-HWP-P in between the laser and the SLM prepares the entire beam in initial state $\hat {V}_0=[0\,\,\,\,\,1\,\,\,\,\,0]^T$. The fast axis of the SLM is oriented along the laboratory horizontal, $\hat {\beta }_f(SLM)=[1\,\,\,\,0\,\,\,\,0]^T$, and the fast axis of the HWP is oriented at 247.5 deg with respect to horizontal so $\hat {\beta }_f(HWP)=(1/\surd \, 2)[-1\,\,\,\,1\,\,\,\,0]^T$. Hence the slow axes are given by $\hat {\beta }_s(SLM)=[-1\,\,\,\,0\,\,\,\,0]^T$ and $\hat {\beta }_s(HWP)=(1/\surd \, 2)[1\,\,\,\,-1\,\,\,\,0]^T$. Let $\alpha (x,y)$ denote the rotation due to the first pass on the SLM and let $\delta (x,y)$ denote the rotation due to the second pass on the SLM, both on a pixel-by-pixel $(x,y)$ basis. The state of polarization at the output of the apparatus is denoted $\hat {V}_{3}(x,y)$ where

Due to physical limitations on the operation of the SLM, arbitrary values of $\alpha$ and $\delta$ cannot be obtained. Instead, the range of $\gamma$ values must be divided into three regions with physically-realizable values of $\gamma$, $\alpha$ and $\delta$ in each region as shown in Fig. 4 and described in detail in Appendix 9.

 

Fig. 4. (a) Required values of rotation by the first pass on the SLM ($\alpha (deg)$) and the second pass on the SLM ($\delta$) needed to yield a given value of $\gamma (deg)$. Analytic expressions for $\alpha$ and $\delta$ as a function of $\gamma$ are given in Appendix 9. (b) Plots of $\alpha$, $\delta$ and $\gamma$ as a function of $y_0(mm)$ where $\gamma$ and $y_0(mm)$ are related by Eqn. (10).

Download Full Size | PDF

An example of a polarization transformation represented on the Poincare sphere is shown in Fig. 5 for $\gamma =60\,deg$. To obtain the values of $\alpha$ and $\delta$ required to produce a given value of $\gamma$ we used the prescriptions provided in [27] and [35].

 

Fig. 5. Example of the 3-step polarization process. The red curve is the transformation by the first pass on the SLM, the blue curve corresponds to the HWP, and the green curve to the second pass on the SLM. The purple line indicates the slow axis of the HWP. The initial state is ${0}=[0\,\,1\,\,0]^T$ and the final state is $[\cos 2\gamma \,\,\,\, \sin 2\gamma \,\,\,\, 0]^T$ where, in this example, $\gamma =60^0$.

Download Full Size | PDF

4. Discussion

4.1 General remarks

This investigation has yielded a number of unusual and unexpected results for the propagation of optical vector beams. In this section we present a preliminary discussion of some of these results.

Curvature of the beam trajectory is clearly a coherent effect if for no other reason than defining a polarization distribution across the beam’s wavefront implicitly requires spatial coherence across the beam. Absent such coherence, it is not clear what the distribution function $\gamma (x,y)$ would mean. This also suggests that at large distances where the beam may experience a loss of transverse spatial coherence and/or polarization the effect would disappear.

Recall also that in forming the eikonal (5) we also derived the equality

Multiplying (6a) by the scalar factor $\epsilon _0/2$ ($\epsilon _0$ is the vacuum permittivity) we obtain a continuity equation in the transverse plane for the average electromagnetic field energy density $\mathcal {E}=\epsilon _0|E|^2/2$ of a monochromatic beam,

The interpretation of Eq. (6b) requires some care. This equation does not say that $\partial \gamma /\partial z=0$ but rather that the total derivative $D\gamma /Dz=0$. Consider the light exiting the SLM from the pixel located at $(x_0,y_0)$ having linear state of polarization $\gamma (x_0,y_0)$ as defined in (10). As the beam propagates, the small parcel of light from $(x_0,y_0)$ travels along a curve parameterized by $z$ and defined by Lagrangian coordinates $(x[x_0,x_0,z],y[x_0,y_0,z])$ where $y[x_0,y_0,z]$ is the $y_L$ of (11). Equation (6b) says that the $\gamma$ of this piece of the beam does not change as it follows its trajectory. Strictly speaking this means that $\gamma$ is fixed as measured by a polarimeter having its input always normal to the beam direction. As a practical matter, the curvature of the beam in this case is so small that errors due to keeping the polarimeter aligned instead along the $z-$axis are small.

From (15) we see that the strength of the quadratic bending is controlled by the quantities $k_0$ and $a$ and hence the beam curvature can be increased by decreasing the optical wavelength and by decreasing the size of the exit aperture $a$. In the latter case, decreasing $a$ clearly increases diffractive broadening where the diffraction angle increases as $1/a$ but that deleterious effect is more than compensated by the $1/a^3$ increase in curvature. Although it is possible to fabricate sub-wavelength apertures that would reduce $a$ by orders of magnitude, a similarly-sized SLM would be required to impose the necessary polarization gradient. Nonetheless, even reducing the aperture to a feasible 2mm would result in >4cm of curvature at the 34m propagation distance we used in this work. It may also be possible to impose the polarization gradient on an Airy beam intensity profile in an effort to mitigate diffraction while maintaining the polarization-induced bending. We leave this question for future consideration.

Equation (6c) for the total derivative of the generalized velocity $\vec v$, reproduced here for convenience,

One can also place these results in the context of energy minimization. Following the works of [38], [39], and [40], if the optical intensity obeys the constraint (6a), an appropriate action is given by

4.2 Beam momentum

We can continue the interpretation of the continuity equations in terms of the transport of various quantities corresponding to fluid-mechanical quantities. For example, if we define an optical linear momentum density by

Next we note that the vector beam created by the apparatus in this experiment is linearly-polarized and hence does not possess spin angular momentum. The phase distribution contains no singularities and, hence, the beam does not possess orbital angular momentum of the type typically attributed to vector OAM beams. Given the quadratic trajectory of the beam in laboratory coordinates it is tempting to assign a mechanical angular momentum to the beam. Using the linear momentum density (25) we could, if desired, define an optical angular momentum density by the cross product

4.3 Spatially-dependent Pancharatnam-Berry phase

The eikonal (5) linking optical path and polarization distribution was derived from the inhomogeneous wave equation along with appropriate assumptions surrounding propagation. However, it is not immediately obvious how the polarization preparation process outlined above can provide the initial non-uniform phase front and thus produce a non-zero “velocity” (4) in free-space. From the transport analysis above we see that the beam propagation can be determined without explicit knowledge the spatial dependence of the overall phase $\phi (x,y,z)$ but it is necessary that $\phi (x,y,z)$ have non-zero spatial derivatives. In this section we provide a physical justification for the existence of a spatially-dependent overall phase, even in free-space, due entirely to the spatial dependence of the pixel-by-pixel polarization transformation process. Details of the calculation are found in Appendix 9.

It is known that any polarization transformation is accompanied by a change in overall phase $\phi$ [28,29], a quantity referred to as the Pancharatnam-Berry (PB) phase or, more generally, as the geometric phase. This phase $\phi$ is proportional to the oriented geometric area $\vec {\Omega }$ enclosed on the surface of the Poincare sphere by the piecewise-continuous "circuit" defining the polarization transformation (Fig. 6). The circuit in this case comprises an ordered set of vertices and oriented curves connecting the vertices corresponding to each individual polarization transformation in the sequence of transformations. Since in general the final SOP will not be identical to the initial SOP the circuit must be closed mathematically (not experimentally) by performing a parallel transport from the final back to the initial state. This final transformation is defined by the shortest geodesic curve linking the final state back to the initial state which, in our experiment is along an equatorial arc. Hence the PB phase acquired in a transformation from an initial to a final polarization state of depends on the transformation path linking the two states. Since the transformation paths in the apparatus are spatially dependent, occurring on a pixel-by-pixel basis, the PB phase $\phi (x,y,z)$ of light exiting the apparatus must also be spatially dependent.

 

Fig. 6. (a) Surface area $\Omega$ on the Poincare sphere enclosed by a piece-wise smooth curve comprising an ordered set of vertices $\hat {V}_{n}$ and connecting curves (polarization transformations).

Download Full Size | PDF

One analytical approach to calculating the PB phase is to employ the local Gauss-Bonnet theorem, and this approach is described in detail in Appendix 9.

5. Summary

In this paper we have predicted an effect whereby a vector beam can self-bend in a manner that depends on the spatial distribution of the linear polarization angle. Using Lagrangian coordinates in conjunction with a “transport” formulation of optical propagation dynamics, the degree of bending was predicted analytically as a function of propagation distance. We then demonstrated the accuracy of these predictions in experiment. While the magnitude of the effect is small (scales as $k_0^{-2}$ as does Airy beam bending) it is certainly measurable, particularly over long distances, and provides a re-configurable approach to beam steering. The model also clearly suggests that other state-of-polarization distributions will yield still other, quantitatively different beam trajectories. More generally, a new mechanism for curving an optical beam in a locally flat spacetime may have implications in other areas of physics. For example, the effect may play a role at cosmological scales where recent evidence [41] indicates that millimeter-wave radiation from black holes exhibits both partial polarization and polarization gradients.

6. Appendix: Review of Lagrangian coordinates

In this appendix we review very briefly the distinction between Eulerian and Lagrangian approaches to transport problems. As with any physics problem a first step is to decide on a means for specifying positions. The almost-universal approach takes the Eulerian viewpoint in which quantities of interest are written as functions of a fixed set of coordinates, e.g., $(x,y,z)$, defined over the entire problem domain. For example, the electric field in (2) is represented in the Eulerian approach where the amplitude, phase, and polarization angle are all written as functions of Eulerian position $(x,y,z)$.

In contrast, the Lagrangian approach treats positions as functions that satisfy the physics of the problem. In classical geometric optics, for example, the method of characteristics defines the path taken by optical rays as the coordinate system for displaying optical intensity (see e.g. [19]). In general, to find this path one minimizes a cost function known as the “action”. For our problem the action is minimized by the coordinates that solve (13) [24]. Given starting point, $(x_0,y_0,z_0)$, the path is denoted $(x(s,x_0,y_0), y(s,x_0,y_0), z(s,x_0,y_0)$ where $s$ is a monotonically-increasing, independent variable that parameterizes the curve. For example $s$ could be time or it could be arc length. A location in Lagrangian coordinates is now a specified function of both the parameter $s$ and the starting location $(x(s,x_0,y_0), y(s,x_0,y_0), z(s,x_0,y_0))$. In standard introductions to this topic quantities of interest are typically parameterized by time, $t$. In our case, however, by assuming the beam was monochromatic, we eliminated time from the problem immediately in going from the full wave equation to the Helmholtz Eq. (1). The beam propagation in this problem is assumed to be paraxial along the $z-$axis, so we can parameterize the ray trajectory of a small volume of the light beam instead by $z$ as $(x(z,x_0,y_0), y(z,x_0,y_0), z)$.

The distinction between the two approaches can be illustrated by the example shown in Fig. 7 which depicts a trajectory originating at $z=0$ and propagating in $z$ along the curved trajectory defined by (2). Points in the Eulerian system (right) are labelled by the fixed coordinates $(y,z)_E$ which contain no information about the propagating light. By contrast, the Lagrangian system coordinates $(y[y_0,z],z)_L$ are defined by the initial conditions $y_0$ and the physics of propagation (obtained by solving 13). It is for this reason that the Lagrangian approach is sometimes said to correspond to to a viewpoint that “rides along” with the fluid flow or, in this case, along with the rays of light. We note that in moving from the Eulerian to Lagrangian description, the total derivatives of Eqs. (6) become ordinary derivatives thereby greatly facilitating an analytical solution.

 

Fig. 7. Illustrating the distinction between Lagrangian (left) and Eulerian (right) viewpoints. In the former case the coordinates are themselves solutions to the governing equations and hence encode the problem physics. By contrast in the Eulerian viewpoint the coordinate system is independent of the physics. For problems involving transport of a quantity (e.g., optical intensity) the Lagrangian point of view can greatly facilitate analysis for reasons discussed in the text.

Download Full Size | PDF

Consider the beam phase $\phi (y,z)$ in Eulerian coordinates, or $\phi (y[y_0,z],z)$ in Lagrangian coordinates and suppose we wish to calculate the change in phase for a small change in $z$. In the Eulerian viewpoint this is simply $\partial \phi (y,z)/\partial z$. However, in the Lagrangian view, since $y$ is itself a function of $z$ we must apply the chain rule, and we obtain instead the quantity

(6.1)$$\dfrac{\partial\phi}{\partial z} + \dfrac{\partial\phi}{\partial y}\dfrac{\partial y}{\partial z} \equiv \dfrac{D\phi}{D z}.$$

That is, since in the Lagrangian view the points are constrained to move along the actual trajectory, the change in $\phi$ due to a small change in $z$ now has two contributions: one from the intrinsic variation in $\phi$ with $z$ and one from the change in the local shape of the trajectory as a function of $z$. Derivatives of this type arise frequently in transport problems, hence the Lagrangian view is particularly applicable for such problems. The derivative in (6.1), usually denoted with capital "D“, is referred to as the total derivative but is also referred to in the literature variously as material derivative, substantial derivative, hydrodynamic derivative, and derivative-following-the-motion [42].

7. Appendix: Diffraction solution

In this Appendix we provide the solution for the beam propagation in the presence of diffraction. In this case, the full intensity profile requires us to solve Eqs. (13). Letting $\epsilon =k_0\sigma ^2$ be twice the Raleigh length, (13b) becomes

(7.1)$$\frac{d^2x_z}{dz^2}=\frac{x_z\epsilon^2}{(z^2+\epsilon^2)^2}$$

which, assuming collimated light ($dx_0/dz=0$) and initial position $x_0$, yields the solution

(7.2)$$\begin{aligned} x_z&=\frac{x_0\sqrt{z^2+\epsilon^2}}{\epsilon}\\ &=x_0\left(\frac{\sigma_z}{\sigma}\right). \end{aligned}$$

In the $\hat {y}$ direction, Eq. (13c) can be similarly solved by first noting that that the terms involving the polarization angle gradient are constant w.r.t. $z$. Thus, defining

(7.3)$$c={-}\frac{1}{k_0^2}\frac{d\gamma(y_0)}{dy_0}\frac{d^2\gamma(y_0)}{dy_0^2}$$

we can write (13c) as

(7.4)$$\frac{d^2y_z}{dz^2}=c+\frac{y_z\epsilon^2}{(z^2+\epsilon^2)^2}.$$

which has the exact solution

(7.5)$$\begin{aligned} y_z&=\frac{c z^2}{2} +\frac{y_0\sqrt{z^2+\epsilon^2}}{\epsilon}+\frac{c\epsilon^2}{2}+\frac{c\epsilon\sqrt{z^2+\epsilon^2}}{2}\left\{(2G-1)+\phantom{\frac{\sqrt{z^2} }{\sqrt{\epsilon^2}}}\right.\\ & \left. \tan^{{-}1}\left(\frac{z}{\epsilon}\right)\log\left[\frac{iz+\epsilon+i\sqrt{z^2+\epsilon^2}}{z-i\epsilon+\sqrt{z^2+\epsilon^2}}\right] -iLi_2\left(\frac{z-i\epsilon}{\sqrt{z^2+\epsilon^2}}\right)+iLi_2\left(\frac{-z+i\epsilon}{\sqrt{z^2+\epsilon^2}}\right)\right\} \end{aligned}$$

where $G$ is Catalan’s constant (equal to $\approx 0.916$) and $Li_n(\cdot )$ is the polylogarithm function of order $n$. From the exact solution, it can be observed that the first two terms constitute the polarization-induced bending and free-space diffraction respectively. The remaining terms comprise a coupling between these two effects, but have little influence on the beam trajectories. For small $z$ the term in brackets sums to nearly zero and for large distances, grows linearly at a rate $c\epsilon z$ (the term in brackets is bounded on the interval $[-1,2G-1]$). These terms are therefore dominated for both short and long distances by the quadratic growth in beam displacement and linear growth in beam width proportional to $y_0 z/\epsilon$. For $z=0$, of course, one recovers the initial conditions $y_0$.

Using our polarization profile (10), the constant becomes

(7.6)$$c=\frac{\pi^2 (a-y_0)}{a^4k_0^2}$$

and we have that the transverse displacement at a distance $z$ becomes (in meters)

(7.7)$$\begin{aligned} y_z&=\frac{z^2\pi^2}{2k_0^2a^4}\left(a-y_0\right)+\frac{y_0\sqrt{z^2+k_0^2\sigma_0^4}}{k_0\sigma_0^2}\\ &=\mu_z\left(1-\frac{y_0}{a}\right)+y_0\left(\frac{\sigma_z}{\sigma_0}\right) \end{aligned}$$

where we have introduced the function

(7.8)$$\mu_z=\frac{z^2\pi^2}{2k_0^2a^3}.$$

The expressions (7.2) and (7.7) along with the initial intensity profile are then sufficient to form (13a). Given the initial Gaussian intensity profile, we form the inverse of the determinant of the Jacobian

(7.9)$$\det|J(x_z,y_z)|^{{-}1}=\frac{a\sigma_0^2}{\sigma_z(a\sigma_z-\mu_z\sigma_0)}.$$

Multiplying by the initial intensity profile, and writing the initial coordinates $x_0,~y_0$ in terms of $x_z,~y_z$ by inverting (7.2) and (7.7) gives finally

(7.10)$$\rho(x_z,y_z)=\frac{a\sigma_0^2}{\sigma_z\left(a\sigma_z-\mu_z\sigma_0\right)}e^{-\frac{x_z^2}{\sigma_z^2}+\frac{a^2(y_z-\mu_z)^2}{(a\sigma_z-\mu_z\sigma_0)^2}}.$$

8. Appendix: Euler-Rodrigues rotation matrix

In this Appendix we provide, for convenience, the general Euler-Rodrigues form of rotation matrices [43]. In real 3-space, a positive rotation (in a right-hand sense) through angle $\theta$ about an axis defined by the vector $\hat {\beta }=[\beta _{1}\,\,\beta _2\,\,\beta _3]^T$ is denoted $R_{\hat {\beta }}(\theta )$ and is given by

(8.1)$$R_{\hat{\beta}}(\theta)=(\cos\theta)I_3+(\sin\theta)\beta^{{\times}}+(1-\cos\theta)\left(\hat{\beta}\hat{\beta}^T\right)$$

where $I_3$ is the $3\times 3$ identity matrix,

(8.2a)$$\beta^{{\times}}= \begin{bmatrix} 0 & -\beta_3 & \beta_2 \\ \beta_3 & 0 & -\beta_1 \\ -\beta_2 & \beta_1 & 0 \end{bmatrix},$$

(8.2b)$$\hat{\beta}\hat{\beta}^T= \begin{bmatrix} \beta_1^2 & \beta_1\beta_2 & \beta_1\beta_3 \\ \beta_2\beta_1 & \beta_2^2 & \beta_2\beta_3 \\ \beta_3\beta_1 & \beta_3\beta_2 & \beta_3^2 \end{bmatrix}.$$

Combining all the terms,

(8.3)$$R_{\hat{\beta}}(\theta)= \begin{bmatrix} \cos\theta+\beta_1^2\left(1-\cos\theta\right) & \beta_1\beta_2\left(1-\cos\theta\right) -\beta_3\sin\theta & \beta_1\beta_3\left(1-\cos\theta\right)+\beta_2\sin\theta \\ \beta_1\beta_2\left(1-\cos\theta\right) +\beta_3\sin\theta & \cos\theta+\beta_2^2\left(1-\cos\theta\right) & \beta_2\beta_3\left(1-\cos\theta\right)-\beta_1\sin\theta\\ \beta_1\beta_3\left(1-\cos\theta\right)-\beta_2\sin\theta & \beta_3\beta_2\left(1-\cos\theta\right)+\beta_1\sin\theta & \cos\theta+\beta_3^2\left(1-\cos\theta\right) \end{bmatrix}.$$

9. Appendix: Details of the polarization calculations including geometric phase

In this Appendix we provide details of the polarization calculations including the Pancharatnam-Berry phase (geometric phase) calculation for each of the three regions of $\gamma$ values.

9.1. Polarization Transformations

In this section we provide the mathematical description of the 3-step polarization transformation performed by the apparatus. The analyses are performed in Stokes space and we employ the following conventions and notations.

• (1) We adopt the $+i\omega t$ convention [44] where the beam propagates in the positive $z-$direction with spatio-temporal phase factor $\exp (i(\omega t-kz))$.
• (2) A polarization state is represented by a Stokes column 3-vector with basis states defined by the following particular indexing of the Pauli spin matrices

  (9.1)$$\begin{matrix} \sigma_1= \begin{bmatrix} 1 & 0\\ 0 & -1 \end{bmatrix} & \sigma_2= \begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix} & \sigma_3= \begin{bmatrix} 0 & -i \\ i & 0 \end{bmatrix}. \end{matrix}$$
• (3) Stokes vectors $\hat {V}$ are represented as points on the unit Poincare sphere and the Cartesian coordinates of a SOP can be written in terms of two angles, say $\lambda$ and $\varphi$, by

  (9.2)$$\begin{matrix} \hat{V}= \begin{bmatrix} s_1\\s_2\\s_3 \end{bmatrix}= \begin{bmatrix} \cos\varphi\cos\lambda\\ \cos\varphi\sin\lambda\\ \sin\varphi \end{bmatrix} & \text{where} & \begin{cases} 0\leq\lambda<2\pi \\ -\pi/2\leq\varphi\leq\pi/2 \end{cases} \end{matrix}.$$
• (4) Each polarization transformation is lossless and is described as a rigid rotation through some angle $\theta$ on the Poincare sphere about an axis defined by vector $\hat {\beta }= [\beta _1 \,\,\,\,\, \beta _2 \,\,\,\,\, \beta _3]^T$.
• (5) Rotations are represented by 3x3 matrices expressed in standard Euler-Rodrigues form [43] A positive, right-hand-sense rotation through angle $\theta$ about an axis $\hat {\beta }$ is denoted $R_{\hat {\beta }}(\theta )$. In the following, the directed axis of rotation $\hat {\beta }$ will always be chosen such that the rotation $\theta$ is right-hand-rule positive.
• (6) In the $+i\omega t$ convention, the operation of a waveplate with slow axis $\hat {\beta }_s$ is a positive, right-hand sense rotation about $\hat {\beta }_s$. However, many manufacturers specify instead the fast axis of a waveplate, $\hat {\beta }_f=-\hat {\beta }_s$. To maintain the convention specified above, if $\hat {\beta }_f$ designates the fast axis of a waveplate producing retardance $\theta$, then the rotation corresponds to $R_{\hat {\beta }_{s}}(\theta )$.

As mentioned in the text, due to physical limitations on the operation of the SLM the range of $\gamma$ values must be divided into three regions with physically-realizable values of $\gamma$, $\alpha$ and $\delta$ in each region as follows, and as shown earlier in Fig. 4.

Region I

(9.3a)$$0\leq\gamma\leq \pi/2$$

(9.3b)$$\alpha =\pi-2\gamma$$

(9.3c)$$\delta=\pi/2$$

Region II

(9.4a)$$\pi/2 <\gamma \leq \pi$$

(9.4b)$$\alpha =2\gamma-\pi$$

(9.4c)$$\delta=3\pi/2$$

Region III

(9.5a)$$\pi <\gamma < 3\pi/2$$

(9.5b)$$\alpha=2\pi+\pi-2\gamma$$

(9.5c)$$\delta=\pi/2$$

9.2. Pancharatnam-Berry (PB) Phase Calculation

Next we discuss an analytical approach to calculating the geometric phase using the local Gauss-Bonnet theorem (GBT) which has two forms depending on the sense in which the circuit is traversed as viewed looking down onto the surface of the Poincare sphere,

i) Counterclockwise circuit

(9.6)$$\Omega=2\pi-\sum_{m=1}^M\int_{\hat{C}_m}\kappa_{g,m}(s) ds-\sum_{n=1}^N\theta_n$$

ii) Clockwise circuit

(9.7)$$\Omega={-}2\pi+\sum_{m=1}^M\int_{\hat{C}_m}\kappa_{g,m}(s) ds-\sum_{l=1}^N\theta_l$$

where

• $\hat {C}_m(s)$ denotes the curve from vertex $\hat {V}_{m-1}$ to vertex $\hat {V}_m$ parameterized by variable $s$, $s_{min}\leq s\leq s_{max}$, where $s_{min}$ and $s_{max}$ depends on the particular curve and the particular $\gamma$ value,
• $\kappa _{g,m}(s)$ is the oriented geodesic curvature of curve $\hat {C}_m(s)$ at $s$, given by $\kappa _{g,m}(s)=\det [\hat {C}_m(s)\,\,\hat {C}'_m(s)$ $\hat {C}''_m(s)]$ and where $\hat {C}'_m(s)$ denotes the derivative of $\hat {C}_m(s)$ with respect to parameter [45–47],
• $\theta _{n}$ is the oriented exterior angle formed at vertex $\hat {V}_n$ between the vectors tangent to $\hat {C}_n$ and $\hat {C}_{n+1}$ and similarly for $\theta _{l}$, as shown in the inset in Fig. 6.

We adopt the following convention:

a) If the curve is traversed counterclockwise, then the area enclosed in the region to the left of the curve is assigned a positive orientation, outward from the sphere center.

b) If the curve is traversed clockwise (as in Fig. 8), then the enclosed area is assigned a negative orientation, inward toward the sphere center.

 

Fig. 8. (a) Surface area $\Omega$ enclosed by a piece-wise smooth curve comprising an ordered set of vertices and connecting curves. The labelling scheme for the PB phase calculation indicating vertices $\hat {V}_n$ and curves $\hat {C}_m$, all on the unit Poincare sphere. Inset shows the definition of exterior angle, using $\theta _3$ as an example, as the angle between tangent vectors at the vertex.

Download Full Size | PDF

Note that if the same circuit is traversed in opposite directions then the sign of $\Omega$ changes. This is true for two reasons. Firstly, upon a change of direction the geodesic curvature $\kappa _{g,m}(s)$ remains unchanged since we always assume non-negative rotations by choosing the appropriate rotation axis. Secondly, the sign of any exterior angle $\theta _n$ changes upon change of direction. Then it is easy to see that, for the same circuit traversed in opposite directions, the areas $\Omega$ in (9.6) and (9.7) have opposite signs. Details of the area calculation, which must be performed piece-wise in each Region I, II and III, are presented in Appendix (9.) but the result is remarkably simple,

(9.8)$$\Omega =A\cos{2\gamma}+B$$

where the constants $A$ and $B$ in each Region are given by

(9.9)$$\begin{matrix} \underline{\text{Region}} & \underline{A} & \underline{B}\\ \text{I} & \pi\left(\dfrac{1}{\sqrt2}-\dfrac{1}{2}\right) & -\dfrac{\pi}{2}\\ \text{II} & \pi\left(\dfrac{1}{\sqrt2}+\dfrac{1}{2}\right) & \dfrac{\pi}{2} \\ \text{III} & \pi\left(\dfrac{1}{\sqrt2}-\dfrac{1}{2}\right) & -\dfrac{\pi}{2}+2\pi. \end{matrix}$$

9.3. Contribution of the geodesic curvatures

The curve $\hat {C}_m(s)$ on the unit Poincare sphere corresponding to the particular transformation from vertex $\hat {V}_{m-1}$ to vertex $\hat {V}_m$ is obtained by operating on vertex $\hat {V}_{m-1}$ with the appropriate $3x3$ rotation matrix from (17), where the rotation matrix is now written in parameterized form. Once the curve $\hat {C}_m(s)$ is defined, the geodesic curvature $\kappa _{g,m}$, is given by (see,for example, [45–47])

(9.10)$$\kappa_{g,m}=\hat{C}_m\cdot\hat{C'}_m\times\hat{C}^{\prime\prime}_{m} =\det\Big[\hat{C}_m\,\,\,\hat{C}'_m\,\,\,\hat{C}^{\prime\prime}_{m}\Big]$$

where

(9.11)$$\begin{matrix} \hat{C}'_m(s)=\dfrac{d{\hat{C}}_m/ds} {|d\hat{C}_m/ds|} &, & \text{and} & \hat{C}^{\prime\prime}_{m}(s)=\dfrac{d\hat{C}'_m/ds}{|d\hat{C}'_m/ds|} \end{matrix}$$

are, respectively, the unit tangent vector and unit derivative of the tangent vector to the curve as a function of $s$. In the determinant in (9.10), the three vectors can, of course, be either all column vectors or all row vectors and the ordering of vectors in (9.10) is consistent with the convention stated above that CCW circuits are assigned positive area and CW circuits are assigned negative area.

We note that, for any geodesic curve, $\hat {C}(s)$ and $\hat {C}''(s)$ are everywhere either parallel or anti-parallel and, hence, $\det [\hat {C}''_m\,\,\hat {C}_m\,\,\hat {C}'_m]=0$. Thus, for the transformation shown in Fig. 8(a), $\kappa _g=0$ for curves $\hat {C}_1$ (along a line of longitude) and $\hat {C}_4$ (along the equator) and these two curves make no geodesic contribution to the area $\Omega$.

Calculation of the geodesic contribution is straightforward but somewhat tedious. It is made more tedious by the fact that, in order to conform to our chosen conventions, each of the three regions defined in (9.2) must be divided into two subregions and treated separately. Hence we will here analyze just one of the six subregions in detail, the subregion of Region I where $45\leq \gamma (deg)\leq 90$ which will be denoted as Region Ib. Extension of the calculations to the remaining subregions is straightforward. Once the oriented surface area $\Omega$ is obtained for a particular transformation, the PB phase is given by $\phi _{PB}=-\Omega /2$. The same procedure is used to find $\phi _{PB}$ for all values of $\gamma$ resulting in the plot shown in Fig. 9(a).

 

Fig. 9. (a) Pancharatnam-Berry phase $\phi _{PB}$ as a function of $\gamma$; (b) $\phi _{PB}$ (blue) plotted versus $y$ with $y$ along the vertical axis. Also plotted is the (normalized) beam power profile (red).

Download Full Size | PDF

Finally, the PB phase $\phi _{PB}(x,y)$ or, more generally, the optical geometric phase, is related to $\Omega (x,y)$ by [28,29]

(9.12)$$\phi_{PB}(x,y)={-}\Omega(x,y)/2.$$

The plots in Fig. 9 shows the PB phase $\phi _{PB}(x,y,0)$ as a function of $\gamma$ and as a function of $y$ for the spatial distribution $\gamma (x,y)$ specified in (10). This non-zero spatial variation in $\phi _{PB}(x,y,0)$ thus contributes a non-zero velocity vector in (4) which now becomes

(9.13)$$\vec{v}_{PB}=k_o^{{-}1}\nabla_X\phi_{PB}\equiv \dfrac{1}{k_o}\left\{\begin{array}{c} 0\\ \partial \phi_{PB}/\partial y \end{array}\right\}.$$

Funding

Office of Naval Research (56-17-01-017).

Acknowledgments

Funding for this work was provided by the Office of Naval Research, project number 56-17-01-017

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. E. Otte and C. Denz, “Optical trapping gets structure: Structured light for advanced optical manipulation,” Appl. Phys. Rev. 7(4), 041308 (2020). [CrossRef]  

2. S. C. Chapin, V. Germain, and E. R. Dufresne, “Automated trapping, assembly, and sorting with holographic optical tweezers,” Opt. Express 14(26), 13095–13100 (2006). [CrossRef]  

3. M. Padgett and R. Bowman, “Tweezers with a twist,” Nat. Photonics 5(6), 343–348 (2011). [CrossRef]  

4. P. A. Santi, “Light sheet fluorescence microscopy: A review,” J. Histochem. Cytochem. 59(2), 129–138 (2011). [CrossRef]  

5. J. Gibson, J. Courtial, M. Padgett, M. Vasnetsov, V. Pas’ko, S. Barnett, and S. Franke-Arnold, “Free-space information transfer using light beams carrying orbital angular momentum,” Opt. Express 12(22), 5448–5456 (2004). [CrossRef]  

6. E. Otte, I. Nape, C. Rosales-Guzman, C. Denz, A. Forbes, and B. Ndagano, “High-dimensional cryptography with spatial modes of light: Tutorial,” J. Opt. Soc. Am. B 37(11), A309–A323 (2020). [CrossRef]  

7. B. Zhao, V. Rodríguez-Fajardo, X.-B. Hu, R. I. Hernandez-Aranda, B. Perez-Garcia, and C. Rosales-Guzmán, “Parabolic-accelerating vector waves,” Nanophotonics 11(4), 681–688 (2022). [CrossRef]  

8. B. Ndagano, H. Sroor, M. McLaren, C. R.-Guzman, and A. Forbes, “Beam quality measure for vector beams,” Opt. Lett. 41(15), 3407–3410 (2016). [CrossRef]  

9. Q. Guo, C. Schlickriede, D. Waang, H. Liu, Y. Xiang, T. Zentgraf, and S. Zhang, “Manipulation of vector beam polarization with geometric metasurfaces,” Opt. Express 25(13), 14300–14307 (2017). [CrossRef]  

10. S. Fu, S. Zhang, and C. Gao, “Bessel beams with spatial oscillating polarization,” Sci. Rep. 6(1), 30765 (2016). [CrossRef]  

11. G. Bautista, J.-P. Kakko, V. Dhaka, X. Zang, L. Karvonen, H. Jiang, E. Kauppinen, H. Lipsanen, and M. Kauranen, “Nonlinear microscopy using cylindrical vector beams: Applications to three-dimensional imaging of nanostructures,” Opt. Express 25(11), 12463–12468 (2017). [CrossRef]  

12. L. Kopf, J. R. D. Ruano, M. Hiekkamaki, T. Stolt, M. J. Huttunen, F. Bouchard, and R. Fickler, “Spectral vector beams for high-speed spectroscopic measurements,” Optica 8(6), 930–935 (2021). [CrossRef]  

13. W. Qiao, T. Lei, Z. Wu, S. Gao, Z. Li, and X. Yuan, “Approach to multiplexing fiber communication with cylindrical vector beams,” Opt. Lett. 42(13), 2579–2582 (2017). [CrossRef]  

14. S. Orlov, V. Vosylius, P. Gotovski, A. Grabusovas, J. Baltrukonis, and T. Gertus, “Vector beams with parabolic and elliptic cross-sections for laser material processing applications,” J. Laser Micro/Nanoeng. 13(3), 280–286 (2018). [CrossRef]  

15. M. R. Teague, “Deterministic phase retrieval: A Greens function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]  

16. D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80(12), 2586–2589 (1998). [CrossRef]  

17. J. D. Achenbach, Wave propagation in elastic solids, Applied Mathematics and Mechanics Series (North-Holland Publishing Company, Amsterdam, 1973).

18. S. Jin, P. Markowich, and C. Sparber, “Mathematical and computational methods for semiclassical schrodinger equations,” Acta Numer. 20, 121–209 (2011). [CrossRef]  

19. H. Liu, S. Osher, and R. Tsai, “Multi-valued solution and level set methods in computational high frequency wave propagation,” Communications in computational physics 1(5), 765–804 (2006).

20. G. A. Siviloglou, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of accelerating airy beams,” Phys. Rev. Lett. 99(21), 213901 (2007). [CrossRef]  

21. J. Baumgartl, M. Mazilu, and K. Dholakia, “Optically mediated particle clearing using airy wavepackets,” Nat. Photonics 2(11), 675–678 (2008). [CrossRef]  

22. S. Jia, J. Vaughan, and X. Zhuang, “Isotropic three-dimensional super-resolution imaging with a self-bending point spread function,” Nat. Photonics 8(4), 302–306 (2014). [CrossRef]  

23. J. Zhou, Y. Liu, Y. Ke, H. Luo, and S. Wen, “Generation of Airy vortex and Airy vector beams based on the modulation of dynamic and geometric phases,” Opt. Lett. 40(13), 3193–3196 (2015). [CrossRef]  

24. J. M. Nichols, T. H. Emerson, L. Cattell, S. Park, A. Kanaev, F. Bucholtz, A. Watnik, T. Doster, and G. K. Rohde, “A transport-based model for turbulence-corrupted imagery,” Appl. Opt. 57(16), 4524–4536 (2018). [CrossRef]  

25. O. Korotkova and E. Wolf, “Changes in the state of polarization of a random electromagnetic beam on propagation,” Opt. Commun. 246(1-3), 35–43 (2005). [CrossRef]  

26. J. M. Nichols, T. H. Emerson, and G. K. Rohde, “A transport model for broadening of a linearly polarized, coherent beam due to inhomogeneities in a turbulent atmosphere,” J. Mod. Opt. 66(8), 835–849 (2019). [CrossRef]  

27. A. Sit, L. Giner, E. Karimi, and J. Lundeen, “General lossless spatial polarization transformations,” J. Opt. 19(9), 094003 (2017). [CrossRef]  

28. S. Pancharatnam, “Generalized theory of interference, and its applications,” Proc. Indian Acad. Sci. 44(5), 247–262 (1956). [CrossRef]  

29. M. V. Berry, “The adiabatic phase and Pancharatnam’s phase for polarized light,” J. Mod. Opt. 34(11), 1401–1407 (1987). [CrossRef]  

30. C. Nore, M. E. Brachet, and S. Fauve, “Numerical study of hydrodynamics using the nonlinear schrodinger equation,” Phys. D 65(1-2), 154–162 (1993). [CrossRef]  

31. S. A. Hojman, F. A. Asenjo, H. M. Moya-Cessa, and F. Sot-Eguibar, “Bohm potential is real and its effects are measurable,” Optik 232, 166341 (2021). [CrossRef]  

32. T. Latychevskaia, D. Schachtler, and H.-W. Fink, “Creating airy beams employing a transmissive spatial light modulator,” Appl. Opt. 55(22), 6095–6101 (2016). [CrossRef]  

33. C. Villani, Optimal Transport: Old and New (Springer-Verlag, Berlin, 2008).

34. C. D. F. W. Dyson and A. S. Eddington, “A determination of the deflection of light by the sun’s gravitational field, from observations made at the total eclipse of may 29, 1919,” Phil. Trans. R. Soc. Lond. A 220(571-581), 291–333 (1920). [CrossRef]  

35. M. Runyon, C. Nacke, A. Sit, M. Granados-Baez, L. Giner, and J. Lundeen, “Implementation of nearly arbitrary spatially varying polarization transformations: an in-principle lossless approach using spatial light modulators,” Appl. Opt. 57(20), 5769–5778 (2018). [CrossRef]  

36. Z. I. Feizulin and Y. A. Kravtsov, “Broadening of a laser beam in a turbulent medium,” Radiophys. Quantum Electron. 10(1), 33–35 (1967). [CrossRef]  

37. P. F. Goldsmith, Gaussian Beam Propagation (IEEE, 1998), pp. 9–38.

38. J. Rubinstein and G. Wolansky, “A variational principle in optics,” J. Opt. Soc. Am. A 21(11), 2164–2172 (2004). [CrossRef]  

39. Y. Brenier, “The least action principle and the related concept of generalized flows for incompressible perfect fluids,” J. Amer. Math. Soc. 2(2), 225–255 (1989). [CrossRef]  

40. J.-D. Benamou and Y. Brenier, “A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem,” Numer. Math. 84(3), 375–393 (2000). [CrossRef]  

41. K. Akiyama, J. C. Algaba, A. Alberdi, et al., “First M87 event horizon telescope results. VIII. magnetic field structure near the event horizon,” Astrophys. J. 910(1), 1–43 (2021). [CrossRef]  

42. R. B. Bird, W. E. Stewart, and E. N. Lightfoot, Transport Phenomena (J. Wiley, 2002).

43. D. Koks, Explorations in Mathematical Physics (Springer, 2006).

44. N. Frigo, F. Bucholtz, G. Cranch, and J. Singley, “On the choice of conventions in polarization evolution calculations and representations,” J. Lightwave Technol. 40(1), 179–190 (2022). [CrossRef]  

45. T. Sochi, Introduction to Differential Geometry of Space Curves and Surfaces (LuLu.com, 2017).

46. M. Spivak, A Comprehensive Introduction to Differential Geometry (Publish Or Perish, Inc., 1999).

47. T. Needham, Visual Differential Geometry and Forms (Princeton University, 2021).