Optomechanics of a stable diffractive axicon light sail

Beamed propulsion of a light sail based on radiation pressure benefits from a passively self-stabilizing “beam riding” diffractive film. We describe the optomechanics of a rigid non-spinning light sail that mitigates catastrophic sail walk-off and tumbling by use of a flat axicon diffraction grating. A linear stability analysis and numerical integration of the coupled translational and rotational equations of motion are examined. Stability is traded against longitudinal acceleration. The examined system achieves 90% of the theoretical longitudinal force limit and stability against a relative sail translation up to 30% of the sail radius when the payload is attached to a long boom.


Introduction
Optical momentum carried by light [1] may be redirected by means of reflection, diffraction, scattering, or absorption to achieve radiation pressure on a material body. The idea of propelling a spacecraft to high velocities using solar radiations was first proposed in the early Twentieth century by Tsiolkovsky [2]. With the advent of lasers in the 1960s, laserpropelled sails reaching relativistic velocities for interstellar travel were proposed [3][4][5][6][7][8][9][10][11][12][13]. One of the many challenges associated with a laser-propelled light sail is the achievement of "beam-riding," i.e., the autonomous ability to remain in the beam path without tumbling or sliding out. A decades-long approach considers shaped mirror structures [14][15][16][17][18][19][20][21][22][23][24]. The ponderomotive (or gradient) force as used in optical tweezers [25] is currently negligible for practical space systems. The proposed use of diffraction to impart optical momentum to a body [26] provides an opportunity to decouple the sail shape from the momentum transfer process, thereby affording a new degree of design latitude. For example, a one-dimensional bi-grating has been explored to demonstrate the principle of self-stability [27,28]. The first experimental verification of this principle, along with the measurements of parametric damping, was reported by Chu [29,30]. Furthermore, advancements in the design and fabrication of diffractive films using metamaterial principles provide opportunities to engineer desired optomechanical and other properties into the functionality of a sail [31][32][33][34][35][36][37][38][39][40][41][42][43][44]. In the near term we envision the integration of diffractive light sail components on future solar sailing missions to help resolve engineering challenges [45] such as attitude control (see for example, Near-Earth Asteroid Scout [46], Solar Polar Imager [47], and Solar Cruiser [48]). a e-mail: prs7786@g.rit.edu b e-mail: grover.swartzlander@gmail.com (corresponding author) This report extends our one-dimensional theoretical investigations of a bi-grating to a two-dimensional axicon grating sail. Section 2 describes the transfer of momentum from the light beam to the sail, making use of sail and observer reference frames moving at small relative velocities, so that the Doppler-shifted wavelength may be ignored. The equations of motion for linear and angular degrees of freedom owing to optomechanical force and torque are described in Sect. 3. A linear stability analysis is described in Sect. 4 where conditions for stable light propulsion are described. Numerical solutions of the equations of motion are presented in Sect. 5, including an analysis of motion in the stable regime. Important findings are summarized in Sect. 6.

Photon momentum transfer to a diffractive sail
Let us consider a laser beam of characteristic radial width w incident upon a sail of radius a. Radiation pressure applies a local force at all sail points, resulting in longitudinal acceleration along the optical axis, lateral force, and torque. The minimum beam size, w 0 (the waist), is positioned at the origin of the observer coordinate system (X, Y, Z ), as illustrated in Fig. 1, and the beam propagates in the Z -direction (the optical axis). The electric field profile of a monochromatic beam of wavelength λ and constant power P may be expressed [49] where k Z = 2π/λ, ω = ck is the angular frequency, c is the speed of light, I 0 (Z ) = 2P/πw 2 (Z ) is the irradiance on the optical axis, w(Z ) = w 0 1 + (Z /Z 0 ) 2 1/2 is the radial beam size, Z 0 = πw 2 0 /λ is the diffraction length, (r, z) = k z (X 2 + Y 2 )/2R(z) − arctan(Z /Z 0 ), and R(Z ) = Z 1 + (Z 0 /Z ) 2 for a TEM 00 Gaussian beam. Assuming the beam is much larger than the wavelength (w 0 >> λ), we ignore the transverse component of the wave vector, k X = ∂ /∂ X and k Y = ∂ /∂Y , which are much smaller than k Z . That is, the paraxial approximation is made such the incident wave vector may be expressed We consider a sail comprised of a reflection grating that diffracts light toward the sail axis when illuminated at normal incidence. That is, the sail functions as an optical axicon (see inset of Fig. 1), having a periodic phase profile, axicon (ρ + ) = axicon (ρ ) = −2π(ρ / ) where ρ = (x 2 + y 2 ) 1/2 . For analytical convenience we assume a single diffraction order, noting that this analysis may be readily extended to include multiple reflection and transmission orders. The axicon grating vector K lies in the plane of the sail and points radially toward the sail axis (see inset of Fig. 1).
The grating vector (see inset of Fig. 1) of the sail is directed radially inward from the center of the sail and is expressed where is the grating period and ψ is the polar angle measured counterclockwise fromx . At normal incidence, the angle betweenẐ andẑ is zero, i.e., the sail normal and incident wave vector are perfectly aligned and the grating functions as a reflective axicon.
For an arbitrary attitude, the momentum imparted to the sail may be determined from the difference of linear photon momenta before and after diffraction. This difference is quantified where k i ( k i ) is the incident wave vector in the sail frame (stationary frame) and k d ( k d ) is the diffracted wave vector in the sail frame (stationary frame). For example, if k d = − k i = −(2π/λ)Ẑ then η = 2Ẑ . We note that for a Doppler-free elastic process | η | = | η|. For an arbitrary sail attitude the method of Euler angles is used to relate the wave vectors in the two reference frames (see "Appendix 1"). However, it is instructive to first consider a sail that is tipped in a single direction as depicted in Fig. 1. Let us therefore set ζ Y = ζ Z = 0 and consider a rotation angle ζ X about theX axis. The angle ζ X represents the attitude of the sail normal (ẑ ) with respect to the beam axis (Ẑ ) and is measured counterclockwise fromẐ , i.e., ζ X < 0 for the attitude of sail shown in Fig. 1. The angle of incidence β i,x is measured counterclockwise from the sail normal such that β i,x = −ζ X and β i,x > 0 for the orientation shown in Fig. 1.
In the sail reference frame the incident wave vector may be expressed The diffracted wave vector k d is determined from the phase matching condition, whereby the phase of the electric field tangential to the sail surface is continuous at the interface: where m is the integer-valued diffraction order. For a normally incident beam where k i · (x +ŷ ) = 0 and k d = −m K , the beam is diffracted toward the sail axis as desired and discussed below when m = −1.
Let us express the components of the diffracted wave vector by use of a unit vectorÂ: where phase matching and elastic scattering ( where the − (+) sign corresponds to a reflection (transmission) grating. To achieve efficient acceleration along the beam axis we assume a reflection grating in this report.
Let us now describe diffraction in the stationary reference frame where k i = (2π/λ)Ẑ and where the unit vectorB is the rotated version ofÂ: General expressions relating rotated vectorsÂ andB are described in "Appendix 2". We therefore find the components of the efficiency vectors:

Optomechanics of a diffractive sail
The force and torque imparted to the sail produce both linear and angular displacements that depend on initial conditions and other factors such as the beam power, sail shape, and the spatial distribution of the grating vector. As depicted in Fig. 1 we assume a rigid circular sail of radius a whose distribution in the sail reference frame F may be expressed: where the function Circ(s) has a value of unity (zero) if |s| < 1(|s| > 1). A payload of mass M p is attached to the sail of mass M s by means of a rigid boom of mass M b and length D b and negligible thickness. A positive (negative) value of D b corresponds to a non-exposed (exposed) payload. For convenience we assume M s = M p such that the center of mass coincides with the mid-point of the boom. For this configuration the principal moment of An observer standing next to a stationary laser system will observe the sail moving through space in the F = (X, Y, Z ) coordinate system, where the reference frame F is described by a right-handed set of unit vectors {X ,Ŷ ,Ẑ } and origin O. We wish to predict the position, velocity, and attitude of the sail in that inertial reference frame. However, radiation pressure exerted on the sail is more readily described in the non-inertial reference frame of the sail, F , with right-handed coordinate system (x , y , z ) and origin O (see Fig. 1). In a homogeneous coordinate system (see "Appendix 3"), an arbitrary point in F F is expressed as a column vector [X, Y, Z , 1] T [x , y , z , 1] T , where the 4th component is a scaling factor set to unity.
Radiation pressure on a sail gives rise to forces and torques that may translate and rotate the sail. The translation of the sail in the frame F may be described by the displacement vector We represent the attitude of the sail in this frame in terms of ZYX sequence of Euler angles {ζ Z , ζ Y , ζ X } (see "Appendix 1"). For an arbitrary rotation and translation the relationship between the two frames of reference may be expressed ⎡ where H is the Homogeneous transformation matrix described in "Appendix 3", and the elements containing factors of c X,Y,Z = cos ζ X,Y,Z and s X,Y,Z = sin ζ X,Y,Z belong to the rotation matrix described in "Appendix 1". The net radiation pressure force in the stationary reference frame is found by integrating over the local force elements: where P F is transformed into the reference frame F by the expression P F = H −1 P F , φ is the angle between the sail normal and the incident wave vector (i.e., cos φ =Ẑ ·ẑ ), I = |E(X, Y, Z )| 2 is the beam irradiance described in Eq. (1), c is the speed of light, and we have applied Newton's second law to the right-hand side where M = M s + M p + M b is the total light sail mass. Unlike the net force, the net torque N net measured about the center of mass of the sail is calculated in the sail reference frame F and may be found by integration: where r = x x + y ŷ − (D b /2)ẑ is the moment arm. Euler's equations for rotational degrees of freedom may be expressed where the angular velocity of the sail measured in the reference frame F is related to the time rate of change of Euler angles (see "Appendix 1").
and where the dot symbol represents the time derivative. The displacement, velocity, attitude, and angular velocity of the sail may be found by simultaneously solving coupled equations, Eqs. (13)- (16).

Linear stability analysis of a diffractive sail
From a practical point of view we desire the sail to accelerate in theẐ direction, while otherwise at an equilibrium position centered on the beam and an equilibrium attitude with the sail axis parallel to the optical axis. To determine whether a given set of system parameters satisfies this requirement, linear stability analysis is applied [50]. Let us define a state vector: The linearized equations of motion for translation and rotation may be expressed: where 0 is calculated at the equilibrium state q 0 = 0: By calculating the eigenvalues of the Jacobian of 0 , we determine complex frequencies that correspond to state solutions having the time-dependent form exp(γ a,b t), where real values of γ a,b provide exponential damping or gain, and imaginary values provide oscillations. Four complex frequencies are found which satisfy: where γ a,r , γ b,r , ω a , ω b are real values. The conditions for linear stability are γ a,r ≤ 0 and γ b,r ≤ 0, i.e., exponential growth is prohibited. For γ a,r = γ b,r = 0 as found below, the sail oscillates about the equilibrium point with four characteristic periods that depend on system parameters such as the grating period, the size of the sail, the beam size and power, and the moment of inertia of the light sail. What is more, for the symmetric system considered in this report 1 = 5 < 0, 4 = 8 = 0, and 2 3 = 6 7 < 0, 2 1 > 4| 2 3 | and we therefore find two degenerate frequencies: a high frequency ω h and a low frequency ω l satisfying where ω 2 0 = 1 + 4 = 5 + 8 and 2 = (( 1 − 4 ) 2 + 4 2 3 ) 1/2 = (( 5 − 8 ) 2 + 4 6 7 ) 1/2 . Therefore we expect the system to display two oscillation modes when excited close to equilibrium.

Numerical solutions
Closed-form solutions of the system equations of motion generally do not exist, and therefore, numerical integration methods must be applied. For a representative non-optimized case we examined a laser-sail system having parameters listed in Table 1. We assumed a beam power We numerically computed Eqs. (13)-(16) for different initial values of linear and angular displacement, plotting the results in Fig. 2. The linear nature of the force and torque nearequilibrium is clearly evident in Fig. 2 for the range |δ X,Y /a| < 0.5 and |ζ X,Y | < 2.5 • . We also observe that the force along the beam axis reaches roughly 90% of the maximum theoretical value of 2P 0 /c. Furthermore, the value of the roll torque N z is zero, and thus the system does not acquire angular momentum about the sail axis. Changing only the angle ζ X (ζ Y ) at equilibrium we also find that the torque N Y (N X ) is zero valued.
A perspective of the net force exerted on the sail at equilibrium is depicted in Fig. 3a where local transverse components of force are displayed, resulting in no net transverse force. Similarly, the net torque exerted on the sail at equilibrium is depicted in Fig. 3b where local transverse components of torque, are displayed, resulting in no net transverse torque. If the sail is displaced from equilibrium to the right, as in Fig. 3c the net force drives the sail to the left. In Fig. 3d the net torque is in the −Ŷ direction.
Values of the slopes at the equilibrium points are obtained from Fig. 2, which along with the mass and moments of inertia in Table 1 allow us to determine the values of j (see km] = 25Z 0 , assuming the beam size is controlled, so that it does not overfill the sail. As expected from our linear stability analysis, the system remains stable under this condition. The acceleration a Z = 0.51[m/s 2 ] may be increased in proportion to the laser power, thereby providing values of v Z that are relevant for orbit-changing maneuvers, although the high oscillation frequencies (see above) may become mechanically intolerable if not damped.
An examination of Fig. 5 indicates that force and torque are nonlinearly related to linear and angular displacements for |δ X,Y /a| 0.5 and |ζ X,Y | 2.5 • . Below these bounds the system may be characterized by linear and torsional spring models with stiffness values equal to the slopes in Fig. 2. Close to the nonlinear bounds the springs become soft and less able Eur. Phys. J. Plus (2020) 135:570 to provide a restoring force or torque. Beyond these bounds the system is driven away from equilibrium. To explore how the departure from linear behavior affects the range of stable motion for the system described in Table 1 we varied the initial conditions across the range δ X,Y ∈ [−a, a], or ζ X,Y ∈ [−10 • , 10 • ], withδ X,Y =˙ X,Y = 0. We then numerically integrated the coupled equations of motion, categorized the observed motion as stable or unstable, and summarized the results in the stability maps shown in Fig. 5. The stable range of linear displacement (assuming ζ X,Y =δ X,Y =˙ X,Y = 0 at t = 0 indicates a stability zone defined by δ 2 2 where the radius 0.3a is significantly smaller than the bound δ X,Y = 0.5a. We attribute this smaller zone to the weak force stiffness at 0.3a and coupling to motion in other degrees of freedom that do not provide an attraction to equilibrium. A linear zone boundary was found when varying both δ X and ζ Y (with other state parameters equal to zero), and is shown in Fig. 5b. An examination of Fig. 5b indicates that the force at ζ Y = 6 • is equal and opposite to the force at δ X = 0.3a, suggesting both a reason and an equivalence for the stability boundaries at δ X = 0.3a and ζ Y = 6 • . The same zone boundary relation was found when varying δ Y and ζ X . According to Fig. 5c the system stability is more robust to simultaneous displacements along and rotations about a common axis. Finally, we explore an example where variations of the boom length and beam size affect stability. In this example, we selected the initial condition: δ X = −δ Y = 0.1a and ζ X = −ζ Y = 1 • . As shown in Fig. 5d the system is generally more stable for long boom lengths, but for a given Fig. 5 Regions of optomechanical stability for a relative linear displacement δ X /a versus δ Y /a, b orthogonal linear displacement and attitude axes, δ X /a versus ζ Y (or δ Y /a versus ζ X ), c parallel linear displacement and attitude axes, δ X /a versus ζ X (or δ Y versus ζ Y ) and d relative laser beam width versus relative boom length w 0 /a versus D b /a beam size there is a minimum boom length below which the system is unstable. For example, if the beam radius equals half the sail radius, w 0 = a/2, as listed in Table 1, we predict a minimum boom length of D b = 10a. In comparison we made our numerical studies in Section 3 for a boom length of D b = 15a, well into the stable regime. We also predict that stable motion may be achieved when the beam overfills the sail (i.e., w 0 > a), but only if the boom length is made significantly larger than the sail radius. For example if w 0 = a, stability requires D b > 28a.

Summary
Diffraction-based light sails provide a design flexibility that is not afforded by reflective sails. This is attributed to the controlled redirection of light by an engineered diffractive surface rather than a deformed reflective surface. We have described the optomechanics of a rigid non-spinning laser-driven sail comprised of a reflective axicon diffraction grating and a payload attached to a boom. A single diffraction order is assumed, producing diffraction toward the optical axis. Such diffraction affords passive stability while also providing longitudinal acceleration along the beam axis. Our example exhibited 90% of the maximum theoretical value of force along the optical axis, with the 10% deficit sacrificed to achieve beam-riding stability against transverse perturbations as large as 30% of the sail radius and attitude perturbations as large as 6 • . Numerical methods were used to integrate the coupled equations of motion, allowing us to illustrate stable oscillations and to map regions of stable and unstable motion. A linear stability analysis predicted as many as four modes of oscillation, reducing to two degenerate frequencies for our symmetric structure. Both the attitude and transverse motion exhibited the two frequencies. The squared frequency was found to increase linearly with beam power. Our optomechanical model may be readily extended to include complex diffractive structures, complex beam shapes, modulated beam power, and a spinning sail. The model requires further work to include the Doppler effect [51], mechanical compliance, and center-of-mass center-of-pressure offsets. Advanced features that may be integrated into the diffractive sail include active attitude control [31,45,52].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix 1: Euler angles
Let F = {X ,Ŷ ,Ẑ } represent the reference frame of a stationary coordinate system and F = {x,ŷ,ẑ} represent a rotated reference frame. The two frames of reference with orthonormal coordinate systems can be aligned by three consecutive rotations given by Euler angles ζ ζ ζ ∈ {ζ Z , ζ Y , ζ X }. As per Euler theorem, two consecutive rotations have to be made about different axis. We make use of the ZY X Euler angle sequence with rotation matrices: where ξ j is the rotation about j-th axis by angle ζ j and c j = cos ζ j and s j = sin ζ j , and j = X , Y , or Z , with the singularity restriction ζ Y ∈ (−π/2, π/2) and ζ X,Z ∈ (0, 2π). The rotated frame may be expressed F = ξ ξ ξ F where Conversely we may write F = ξ ξ ξ −1 F where The angular velocity depends on the time rate of change of Euler anglesζζ ζ = [ζ X ,ζ Y ,ζ Z ] T whereζ j = dζ j /dt. For the ZY X sequence of Euler angles, the angular velocity in the two frames may be expressed [53]: and the angular velocity vectors are expressed

Appendix 2: Rotation of diffraction efficiency vectors
For an arbitrarily orientated sail with respect to a stationary frame, the incident beam, grating vector, and diffracted beam are expressed: where In the stationary frame, the diffracted beam components are expressed where ⎡ The photon momentum transfer efficiency imparted to the sail may be expressed η = k i − k d /(2π/λ) and η = k i − k d /(2π/λ) (31) or and Similarly, a second frame F with origin O and right-handed basis {x ,ŷ ,ẑ } is expressed A displacement δ δ δ = [δ X , δ Y , δ Z ] T and rotation ξ ξ ξ of F with respect to F may be written as a matrix operation: where H is the homogeneous transformation matrix [55] and ZYX sequence of Euler angles is assumed to calculate the rotation matrix ξ ξ ξ (see Eq. (23)). The inverse mapping from F to F may be expressed: