Skew line ray model of nonparaxial Gaussian beam

Many ray-optics models have been proposed to describe the propagation of paraxial Gaussian beam. However, those paraxial ray-optics models are inapplicable to the beams that violate the paraxial condition. In this paper, we present a skew line ray (SLR) based model to represent the propagation properties of nonparaxial Gaussian beam under the oblate spheroidal coordinates. The free-space evolution of complex wavefront of the light beam including amplitude and phase is derived via this model. Our analysis demonstrates that the SLR model is available for both nonparaxial and paraxial conditions, and can be used to precisely describe the propagation of complex wavefront. Furthermore, this model changes the transverse density of rays while propagating. The behavior influences the transverse intensity distribution and makes the optical rays become concentrated towards the center. We believe that this ray-optics model can be further developed to describe other kind of structured beams such as Laguerre-Gauss and Bessel-Gauss beams. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (080.7343) Wave dressing of rays; (350.5500) Propagation. References and links 1. L. Novotny, E. J. Sanchez, and X. S. Xie, “Near-field optical imaging using metal tips illuminated by higherorder Hermite—Gaussian beams,” Ultramicroscopy 71(1-4), 21–29 (1998). 2. A. Tycho, T. M. Jørgensen, H. T. Yura, and P. E. Andersen, “Derivation of a Monte Carlo method for modeling heterodyne detection in optical coherence tomography systems,” Appl. Opt. 41(31), 6676–6691 (2002). 3. J. C. Ricklin and F. M. Davidson, “Atmospheric turbulence effects on a partially coherent Gaussian beam: implications for free-space laser communication,” J. Opt. Soc. Am. A 19(9), 1794–1802 (2002). 4. M. Toyoshima, T. Jono, K. Nakagawa, and A. Yamamoto, “Optimum divergence angle of a Gaussian beam wave in the presence of random jitter in free-space laser communication systems,” J. Opt. Soc. Am. A 19(3), 567–571 (2002). 5. A. Ashkin, “Forces of a Single-Beam Gradient Laser Trap on a Dielectric Sphere in the Ray Optics Regime,” Biophys. J. 61(2), 569–582 (1992). 6. D. J. Stevenson, F. Gunn-Moore, and K. Dholakia, “Light forces the pace: optical manipulation for biophotonics,” J. Biomed. Opt. 15(4), 041503 (2010). 7. M. Lax, W. H. Louisell, and W. B. McKnight, “From Maxwell to paraxial wave optics,” Phys. Rev. A 11(4), 1365–1370 (1975). 8. T. Takenaka, M. Yokota, and O. Fukumitsu, “Propagation of light beams beyond the paraxial approximation,” J. Opt. Soc. Am. A 2(6), 826–829 (1985). 9. H. Laabs, “Propagation of Hermite-Gaussian-beams beyond the paraxial approximation,” Opt. Commun. 147(13), 1–4 (1998). 10. G. P. Agrawal and D. N. Pattanayak, “Gaussian beam propagation beyond the paraxial approximation,” J. Opt. Soc. Am. 69(4), 575–578 (1979). 11. C. G. Chen, P. T. Konkola, J. Ferrera, R. K. Heilmann, and M. L. Schattenburg, “Analyses of vector Gaussian beam propagation and the validity of paraxial and spherical approximations,” J. Opt. Soc. Am. A 19(2), 404–412 (2002). 12. K. Duan, B. Wang, and B. Lü, “Propagation of Hermite-Gaussian and Laguerre-Gaussian beams beyond the paraxial approximation,” J. Opt. Soc. Am. A 22(9), 1976–1980 (2005). 13. B. Z. Wang, Z. G. Zhao, B. D. Lü, and K. L. Duan, “Vectorial Hermite–Laguerre–Gaussian beams beyond the paraxial approximation,” Chin. Phys. 16(1), 143–147 (2007). 14. J. Broky, G. A. Siviloglou, A. Dogariu, and D. N. Christodoulides, “Self-healing properties of optical Airy beams,” Opt. Express 16(17), 12880–12891 (2008). 15. Y. Kaganovsky and E. Heyman, “Wave analysis of Airy beams,” Opt. Express 18(8), 8440–8452 (2010). Vol. 26, No. 3 | 5 Feb 2018 | OPTICS EXPRESS 3381 #313772 https://doi.org/10.1364/OE.26.003381 Journal © 2018 Received 17 Nov 2017; revised 22 Jan 2018; accepted 24 Jan 2018; published 31 Jan 2018 16. M. A. Alonso and M. A. Bandres, “Generation of nonparaxial accelerating fields through mirrors. I: two dimensions,” Opt. Express 22(6), 7124–7132 (2014). 17. M. A. Alonso and M. A. Bandres, “Generation of nonparaxial accelerating fields through mirrors. II: Three dimensions,” Opt. Express 22(12), 14738–14749 (2014). 18. J. Arnaud, “Representation of Gaussian beams by complex rays,” Appl. Opt. 24(4), 538–543 (1985). 19. Y. A. Kravtsov, G. W. Frobes, and A. A. Asatryan, “Theory and applications of complex rays,” Prog. Opt. 39, 1–67 (1999). 20. B. T. Landesman, “Geometrical representation of the fundamental mode of a Gaussian beam in oblate spheroidal coordinates,” J. Opt. Soc. Am. A 6(1), 5–17 (1989). 21. B. Shi, Z. Meng, L. Wang, and T. Liu, “Monte Carlo modeling of human tooth optical coherence tomography imaging,” J. Opt. 15(7), 075304 (2013). 22. M. V. Berry and K. T. McDonald, “Exact and geometrical optics energy trajectories in twisted beams,” J. Opt. A, Pure Appl. Opt. 10(3), 035005 (2008). 23. M. A. Alonso and M. R. Dennis, “Ray-optical Poincaré sphere for structured Gaussian beams,” Optica 4(4), 476–486 (2017). 24. B. T. Landesman and H. H. Barrett, “Gaussian amplitude functions that are exact solutions to the scalar Helmholtz equation,” J. Opt. Soc. Am. A 5(10), 1610–1619 (1988). 25. Z. Ulanowski and I. K. Ludlow, “Scalar field of nonparaxial Gaussian beams,” Opt. Lett. 25(24), 1792–1794 (2000). 26. G. Rodríguez-Morales and S. Chávez-Cerda, “Exact nonparaxial beams of the scalar Helmholtz equation,” Opt. Lett. 29(5), 430–432 (2004). 27. D. A. McNamara, C. W. I. Pistorius, and J. A. G. Malherbe, Introduction to the Uniform Geometrical Theory of Diffraction (ARTECH HOUSE, 1990), Chap. 2. 28. L. A. Vainshtein, “Open resonators with spherical mirrors,” Sov. Phys. JETP 18, 471–479 (1964). 29. M. Dietrich, Light Transmission Optics (Van Nostrand Reinhold, 1972), Chap. 6.


Introduction
Gaussian beams are widely used in optical imaging [1,2], optical communications [3,4] and optical manipulations [5,6].The physical expression of Gaussian beams is of great significance for theoretically investigating the beams' propagation or inter-actions between beam and other objects.Since nonparaxial form of Gaussian beams provides a universal way to describe beams' propagation, it has been studied for years.Lots of approaches have been proposed to describe the propagations of nonparaxial Gaussian beams, such as the perturbation method [7], power-series expansion [8,9], angular spectrum representation [10,11] and integral solution of the wave equation [12,13].All of them involve complicated mathematical algorithms based on wave optics.Notably, geometrical optics based description of beam propagation provides a convenient and intuitive explanation for many aspects of light beams, such as "self-healing" [14,15], "self-accelerating" [16,17] as well as the Gouy and Berry phases (also, geometrical phases).Thus, the geometrical optical field is of continuous interest, adding to the momentum of the rapidly growing field of structured light beams.Meanwhile, a ray-optics model of nonparaxial Gaussian beam offers a rather convenient way to reveal the beam's propagation behaviors.
Ray-optics models of paraxial Gaussian beam have been studied for years, including the complex ray representation of Gaussian beam [18,19], the skew line ray (SLR) model of Gaussian beam [20,21], the Poynting vector model of Laguerre-Gauss beam and Bessel-Gauss beam [22], and the Poincaré sphere ray-optics model of structured Gaussian beams [23].Notably, the SLR model regards the skew lines (straight generatrixes) of one-sheet hyperboloid as the trajectories of light rays.However,those ray-optics models are based on the paraxial beam and fail when paraxial condition is broken.Paraxial beam is the solution of scalar Helmholtz equation on paraxial approximation.In other words, nonparaxial beam is the exact solutions of scalar Helmholtz equation.Note that the exact solutions of scalar Helmholtz equation have been given in oblate spheroidal coordinates [24][25][26].On the paraxial approximation, the zero order of the solutions tends to be a fundamental Gaussian beam.Undoubtedly, this beam can be regarded as a reference of nonparaxial Gaussian beam.Considering that one-sheet hyperboloid is a component of the oblate spheroidal coordinates.In the spirit of SLR model, here we present a SLR model of nonparaxial Gaussian beam.This model is based on skew lines of one-sheet hyperboloid in the oblate spheroidal coordinates.
Since the skew lines are regarded as trajectories of light rays, the geometrical optical field of nonparaxial Gaussian beam, including amplitude and phase, is traveling along these rays.Our mathematical derivations demonstrate that the expression of geometrical optical field, derived by SLR model, is consistent with the result obtained by wave optics for both nonparaxial and paraxial conditions.The free-space evolutions of amplitude and phase of beam are given in terms of rays.Further study shows that the transverse distributions of ray density change with beam propagation in the nonparaxial case, leading to the corresponding changes of intensity distribution in the transverse plane.By contrast, on the paraxial condition, the change of the ray density in the transverse plane is so small that it can be ignored.Therefore, our SLR model enables precise description of the propagation behaviors of nonparaxial light fields without explicit diffraction calculations, and provides new insight into the physics of structured light beams.
The structure of this paper is arranged as follows.In the next section, we give a brief introduction of the nonparaxial Gaussian beam under oblate spheroidal coordinates.In Section 3, the geometrical optical fields of the nonparaxial Gaussian beam are deduced via SLR model.In Section 4, the evolution of SLR model on the paraxial approximation is demonstrated.An explanation of the broken self-similar property on nonparaxial regime is given in terms of rays in Section 5.The relationship between ellipsoid wavefront and rays is revealed in Section 6.The last section includes our conclusions.

Gaussian Beam in oblate spheroidal coordinates
The oblate spheroidal coordinates are orthogonally built on confocal oblate spheroids and one-sheet hyperboloids in three-dimensions.The z-axis is defined as the rotation axis, and the distance between foci is set to be 2D.Therefore, the scalar Helmholtz can be separated as the product of the spheroidal radial functions and spheroidal angular functions [24,26].The set of coordinates ( , , )  ξ η θ describes a unique point in Cartesian coordinates ( , , ) x y z with.
( )( ) ( )( ) .Note that all points with a nonzero constant η form a one-sheet hyperboloid, and all points with nonzero constant ξ form an oblate spheroid.
A family of exact solutions of the scalar Helmholtz equation in the oblate spheroidal coordinates have been derived by Landesman and Barrett [24].The zero-order solution is given by is the wave number, λ is the wavelength in vacuum.On paraxial approximation, Eq. (2) corresponds to a fundamental Gaussian beam.The relationship between parameter D, beam waist radius 0 ω and half divergence angle ε is determined by [24,26] ( ) As mentioned in [24], this solution is a nonparaxial Gaussian beam.The wave front is an oblate ellipsoid (ξ = nonzero constant).When 0 ξ = , the beam has a plane wavefront at 0 z = ., ,0

The skew lines of a hyperboloid
According to Eq. ( 1), one-sheet hyperboloid with 1. 1 As shown in Fig. 1, considering that an arbitrary point ( ) is located in waist plane, the parameter Q η of the hyperboloid, containing point Q, is derived by Eq. ( 4) as 1 .
While there are two skew lines cross the same point Q, we define the QS 1 (orange line) as the negative skew line Q l − of point Q, and QS 2 (blue line) as the positive skew line Q l + .These two skew lines are given as : , , , τ is a parameter.We assume that the beam propagates in a linear, homogeneous and isotropic medium, so that rays are straight lines.The skew lines can be regarded as the trajectories of rays of Gaussian beam.Both the amplitudes and phases of optical fields travel along those rays.Since two skew-line rays are associated with point ( , , ) S x y z , these two rays launch from different points in the waist plane.The geometrical optical field at point S is the superposition of the fields carried by these rays.As shown in Fig. 2, assuming that rays are given by [27] ( ) ( )

Amplitude propagation along the SLR
ϕ ) are the amplitude and phase at initial point Q 1 (Q 2 ) of ray, respectively.According to the symmetry of fundamental Gaussian beam, , and Phase l ϕ accumulates with the ray's traveling length l.Variable J denotes the determinant of Jacobian matrix and is determined by 27].Note that variables x, y and Q η are the functions of Q x , Q y , thus J can be derived from Eqs. ( 5), (6-1) and (6-2) as The subscript index in Eq. ( 9), of Q η denotes the specific point.
Since there are two skew lines from a same initial point (Fig. 1), both of them play the same role in wave propagation.The weight of field' amplitude at initial point should be divided equally to two rays.A factor of 1/2 should be multiplied at initial amplitude in Eqs.(8-1) and (8-2).Assuming geometrical optical field propagates a tiny length dl along a SLR of QS in Fig. 3, the phase changes d l ϕ .From Eq. ( 2), we can obtain that ( ) ( )

Phase propagation along the SLR
Combining with Eq. ( 10), we get arc tan .
Equation (12) shows that the phase shift along a ray can be divided into two parts.One is the Berry phase shift of b kl ϕ = .The other is the Gouy phase shift of ( )   which is related with rotation angle of l in terms of the initial point.A geometrical explanation of Gouy phase shift has been given by Landesman [20], as shown in Fig. 3 From Eqs. ( 7) to ( 13), the geometrical optical field of point ( , , ) S x y z should be Here, n Q A is the amplitude at initial points, and . According to Eq. ( 2), we get the initial amplitude ( ) ( ) . If we make the substitution of and Eq. ( 11) into Eq.( 14), it can be rewritten as This equation has the same form as Eq. ( 2).An intriguing phenomenon is that the parameter Q η in Eq. ( 5) must have a real root if the equation of skew line is to be determined.That means the initial points of skew line are restricted in a circle area with radius equal to D. In the subsequent of this paper, we call this circle area the "foci disk".This feature indicates that fields of nonparaxial Gaussian beam in other transverse planes ( 0 z ≠ ) are traced by rays from the confined foci disk.

Paraxial form of SLR model
then the beam propagates on the paraxial approximation.The condition of 1 kD  means that the size of beam waist is large enough compared with its wavelength.The condition of 1 η → makes those one-sheet hyperboloids tend to be cylinders which are very close to z-axis.
From the conditions of ( 17), we have divergence angle 0 , the parameter D tends to be where r z is the Rayleigh distance.On this approximation, we have ξ can be regard as a spherical surface at vicinity of the vertex 1 C .The radius of this spherical surface is ( ) R z at point 1 (0, 0, ) C z .Combining Eqs. ( 1) with ( 18), we can derive ( ) then the distance . Therefore, the optical path length l in Eq. ( 16) tends to be ( ) Meanwhile, the skew lines in Eqs.(6-1) and (6-2) have the paraxial form of : , , , Equations (21-1) and (21-2) are consistence with the skew line equation in [2,21].Then the determinant J of Jacobian is .In this case, the total phase shift l ϕ along a skew line is ( ) We choose the ( )  for the initial amplitude, the paraxial geometrical optical field at point ( , , ) S x y z is For paraxial Gaussian beam, ( ) ( ) x y ω ω = + + according to Eqs. (21-1) and (21-2).Then Eq. ( 24) can be rewritten as Obviously, Eq. ( 25) is the field of the well-known paraxial Gaussian beam [29].
In paraxial condition, the restriction of Eq. ( 5) vanishes for the radius of foci disk is much larger than beam waist.Although the paraxial SLR model can be enforced to trace the propagation of rays in nonparaxial regime, the initial points of rays might be out of the foci disk.Besides, we note that the optical path length calculated by Eq. ( 20) is an approximate treatment in paraxial condition.In practical ray tracing methodology, the optical path length is always rigorously equal to the distance between the initial and terminal points of a ray.In order to distinguish these two models, we define that the ray trajectories traced by Eqs.(6-1) and (6-2) is nonparaxial SLR model and by Eqs.(21-1) and (21-2) are paraxial SLR model, respectively.

Ray-based explanation on broken self-similar property
In previous works on nonparaxial Gaussian beam, the beam profiles are quite different from those of the conventional transverse pattern of paraxial beam [12,13].This phenomenon indicates that self-similar property of beam is broken in nonparaxial region.However, none explanation is given.Combining with ray tracing methodology, here we give an explanation in terms of rays.In the following simulations, we choose , the outer envelope is one-sheet hyperboloids as shown in Fig. 5(a).Trajectories of rays of paraxial SLR model for both kD = and 1 kD = (not shown) are consistent with those in Fig. 5(a).For 1 kD = , as shown in Fig. 5(b), the outer envelope is similar to that in Fig. 5(a) because of the nondimensionalization of coordinates.However, the distribution of rays' tips is more concentrated in central zone than that in Fig. 5(a).For detailed discussions, we set the rays with uniform distribution in the initial plane ( 0 z = ), and trace the rays' distribution in an observation plane as shown in Fig. 6.Gray images denote the ray's distribution pattern.The gray value at each pixel represents the ray's density of that area.Both x and y axes are nondimensionalized by  , the rays' distribution after propagation is similar to the initial distribution of rays.This property ensures the self-similar feature with beam propagation.Hence, the intensity distributions are modulated by rays' distributions.In this condition, the intensity patterns traced by the two models should be similar.6(e) has a bright flare at the center and gradually dim toward outside.There is no obvious edge of the distribution.Compared with the initial rays' distribution in Fig. 6(a), the distribution is obviously distorted after beam's propagation.The self-similar feature of beam is therefore broken.The intensity pattern traced by nonparaxial SLR model should be more concentrated than that of paraxial SLR model.Since self-similar property requires that the transverse distribution of rays should also propagate in a self-similar way [23].The change of distribution of rays' tips in Fig. 6(e) indicates that self-similar property is broken in nonparaxial case.
To compare intensity patterns traced by two models, we ought to choose the same distribution of amplitude in the initial plane.Note that Eq. ( 2) has amplitude singularity in 0 z = plane which can cause intensity distortion in other planes after ray tracing.Although the modified amplitude has been proposed in [25], it is a superposition of an outgoing wave and an incoming wave, and it has infinite energy.For practical purposes, we use the Gaussian distribution  I is narrower than 0,1 I due to the difference of the initial amplitudes between Eq. ( 2) and Gaussian beam.The curve 2 I is fitted with 0,2 I .The curve 1 I is obviously more concentrated in central zone than curve 2 I , and the peak of 1 I is also higher than the peak of 2 I .These appearances are consistent with discussions in paragraph above.While the paraxial condition is broken, the self-similar feature is vanished.The paraxial SLR model has a great deviation than nonparaxial SLR model, and makes the intensity distribution broaden and weaken.

Ellipsoid wavefront as the envelope of equidistance rays
Considering that the wavefront can be regarded as an envelope of rays' tips with equal length.We plot the rays' trajectories for 100 kD = and 1 kD = , where the lengths of rays equal to λ .
The trajectories of rays in each case are separately traced by paraxial and nonparaxial SLR models, as shown in Fig. 8.The transverse coordinate is nondimensionalized by

Fig. 2 .
Fig. 2. Sketch of the rays arriving at point S outside the waist plane.

1 Q l + and 2 Ql
through the observation point S, the equation of − can be derived by Eqs.(6-1) and (6-2).In the following text, the subscript "GO" denotes the optical field derived by geometrical optics.Then the geometrical optical field Ψ go (S) at point S is the sum of fields of the rays associated with point S[15,27].It can be expressed as , we choose the waist plane of 0 z = as the initial plane, so the ( )
and S η are the oblate spheroidal coordinates of point S. Note that the optical path length l (Eikonal) is equal to the distance between points Q and S. From Eqs. (6-1) and(6 . The edge ray (QS) from point Q has arrived at point S on a wave front S ξ .The azimuth angle from point Q to S is α.However, the central ray (OC) arrives at point C at the same time.If the edge ray passes through plane C, it should travel an additional path 0 Δl .The 0 Δl will lead to the azimuth angle difference Δα between S' and S. It can be find that arc tan .

Fig. 4 . 1 ,
Fig. 4. Sketch of the ray tracing along the skew line.As shown in Fig. 4, points ( ) , , S x y z and

1 S 1 Sξ . For the condition of 1 η → , the spheroidal surface 1 S
S can be replaced by the distance Δd between the transverse plane of S and spheroidal surface Eqs. (13), (21-1) and (21-2), the Gouy phase is given by ( ) g arc tan / r z z ϕ = a paraxial beam.In this case, the half divergence angle 8 of the nonparaxial SLR model are presented in Fig.5.Those rays are traced by Eqs.(6-1) and (6-2) and the initial points of rays are arranged in equidistant concentric circles.The coordinates are nondimensionalized by 0 1/ ω for x-/y-coordinates and by 1/ r z for z-coordinate.In the case of 100 kD =

Figure 6 (
a) shows the distribution of rays in the initial plane ( 0 z = ).The rays are uniformly distributed within the circle radius of 0 ω .Figures 6(b)-6(e) are distributions of rays in the observation plane the beam with 100 kD = in Figs.6(b) and 6(c), ray distributions are traced by paraxial and nonparaxial SLR model, respectively.In this case, the gray values are still uniformly distributed in a circle with radius approaches to 0 5ω .When 100 kD =

For 1
kD = , rays' distributions in Figs.6(d) and 6(e) are traced by paraxial and nonparaxial SLR models, respectively.Although 1 kD = breaks the paraxial condition, Fig. 6(d) has little difference with Fig. 6(b) which comes from the paraxial SLR model.

Figure
in the waist plane, where r is the distance from the initial point of ray to the origin.

Fig. 7 .
Fig. 7. Comparison of the intensity distribution in the plane of 5 r z z = between nonparaxial