{'fraction_non_alphanumeric': 0.027000691519456842, 'fraction_numerical': 0.011771547117621173, 'mean_word_length': 4.796046638732009, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structural Stability of a Lightsail for Laser-Driven Interstellar Flight

The structural stability of a lightsail under the laser flux necessary for interstellar flight is studied analytically and numerically. A sinusoidal perturbation is introduced into a two-dimensional thin-film sail to determine if the sail remains stable or if the perturbations grow in amplitude. A reflective material that gives specular reflection of the laser illumination is assumed in determining the resulting loading on the sail, although other reflection models can be incorporated as well. The quasi-static solution of the critical point between shape stability and instability is found by equating the bending moments induced on the sail due to radiation pressure with the restoring moments caused by the strength of the sail and the tension applied at the edges of the sail, permitting analytical expressions for the elastic modulus and boundary tension magnitude to be found as a function of sail properties and the amplitude and wave number of the initial perturbation. These expressions are also derived from a formal variational energy approach. A numerical model of the sail dynamics is developed by discretizing the sail into finite elements. By introducing torsional and rectilinear springs between the elements into the numerical model, a hierarchy of models is produced that can incorporate the effects of bending and applied tension. The numerical models permit the transient dynamics of a perturbed sail to be compared to the analytic results of the quasi-static analysis. The analytic expressions can correctly predict the stability boundary found in the numerical simulations. The stiffness required to make a sail stable against uncontrolled perturbation growth is shown to be unfeasible for known materials. A modest tensioning of the sail (e.g., via an inflatable structure or spinning of the sail) can maintain a stable sail shape under all wavelengths and amplitudes of perturbations.


Introduction
Concentrating laser-light energy onto a reflective foil to permit fast transportation within the solar system and beyond has been actively considered since the 1980s [1,2] with serious academic discussions dating as far back as the 1960s [3,4]. Because of the exponential rate of development of fiber-opticbased lasers within the telecommunication and laser machining industries, the laser-driven spacecraft is now steadily turning from concept to a present-day reality. Arbitrarily large laser beams can now be made by constructing phased arrays of lasers using inexpensive optical components [5,6]. With modern initiatives such as the Breakthrough Starshot project, which intends to use a ∼1 gram lightsail laser-accelerated to 20% the speed of light to reach the nearest star Proxima Centauri within 20 human years, interstellar flight in the 21st century may become an actuality, but for this to happen various interdisciplinary scientific and engineering challenges still need to be overcome [5][6][7]. One such challenge that needs to be solved, given the large laser intensities involved (∼10 GW/m 2 ) and the ideally low sail inertia, is the problem of the dynamic and structural stability of the lightsail.
The directional or beam-riding stability of laser-driven lightsails has already been studied for a variety of rigid sail shapes. The beam-riding stability of conical and spherically curved lightsails attached to the spacecraft via a rigid boom have been investigated [8,9], and Srinivansan et al. have shown that a hyperboloid shaped lightsail impinged upon its convex surface by the laser beam has passive directional stability [10]. Manchester and Loeb also demonstrated that a spherically shaped lightsaila lightsail akin to a balloon-is inherently directionally stable provided that the laser beam intensity profile was Gaussian multi-modal or donut-shaped [11]. The influence of dynamic dampening upon directional lightsail stability has also been considered recently [12,13]. The inclusion of dissipation effects appears to effectively dampen the lightsail motions lateral to the beam axis thereby further increasing directional stability. Shirin et al. have numerically shown, for example, that dynamically dampened conical-shaped lightsails become exponentially stable-as opposed to their marginally stable undampened counterparts [12]. Directional stability has also been shown to be possible without recourse to curved sail shapes by means of engineered optical lightsail surfaces and materials. This has been achieved through the addition of nanoscale (dielectric) structures to the surface of a flat lightsail, a procedure enabled through the recent advances in optical design and nanofabrication. These so-called photonic metasurfaces allow for control over the magnitude and direction of the reflected and transmitted incident laser light while adding relatively little mass to the lightsail because of their nanoscopic nature. Swartzlander et al., in a series of theoretical, computational, and experimental studies have shown that a flat lightsail whose reflecting surface is equipped with diffractive gratings is directionally stable [14][15][16][17][18]. Myilswamy et al. have also shown that a nonlinear photonic crystal lightsail can help minimize the dynamic asymmetry caused by the atmospherically distorted laser beam [19]. Lightsails employing  Bloch-wave type scatterers and other engineered optical metasurfaces have also been reported to be directionally stable [20][21][22]. Santi et al. have also suggested the use of thin-film multilayered optical structures (Bragg mirrors) for the actualization of curved lightsails with enhanced passive beam-riding stability [23].
The studies mentioned so far have assumed the lightsail to be rigid and either ideally flat or perfectly smooth, that is, absent of deformation of the ideal sail shape. Deformation of the lightsail has been built into several other models in an attempt to quantify their influence upon the orientation and trajectory of solar-driven sails [24][25][26]. Huang et al., for example, inquired into the deviation of the resultant solar radiation pressure force due to sail deformation effects caused by wrinkling and billowing via the use of point cloud and triangular mesh methods [26]. Structural analysis of the sail, considering beams/booms and membranes, has also been undertaken [27][28][29][30][31]. Liu et al. studied the attitude dynamics and the vibrations of a square solar sail supported by four beams, the presence of which allowed them to neglect the detailed sail membrane vibrations and wrinkle effects [28,29]. Wong and Pellegrino inquired theoretically, numerically, and experimentally into the visible membrane wrinkling amplitude and wavelength growth when tension is gradually applied to the corners of an initially flat, square solar sail membrane [27]. Other studies concerning the problem of sail structural response have also been undertaken for particular sail shapes [32,33]. Cassenti and Cassenti have also proposed the tensioning of the lightsail via the use of a boundary ring-much akin to a drumskin-as a potential solution to the structural vibration problem caused by the presence of a non-uniform laser beam [34]. A more recent study inquired into the thermal and mechanical stresses that the lightsail experiences during laser-driven acceleration and has concluded that spherically curved lightsails of appreciable curvature are better suited to sustain the large laser loads [35]. Of note are also the experimental efforts of Myrabo et al. who conducted investigations on the problem of sail stability by subjecting lightsail prototypes to laser loads in vacuum [36].
Recently, the application of the radiation pressure regime for laser-driven acceleration of thin foils under extremely intense fluxes (exceeding 10 21 W/cm 2 ) has been considered as a technique for heavy ion acceleration [37][38][39]. Such radiation pressure acceleration technology could overcome the velocity limitations of more traditional pulsed-laser acceleration technologies, such as ablation, and has potential application to fast ignition in inertial confinement fusion. Under such intense fluxes, the structural stability of the thin-foil lightsail under radiation pressure can be treated using the Rayleigh-Taylor formalism of interface instability [40]. The regime of laser flux in this application is fifteen orders of magnitude greater than that proposed for the laser-driven interstellar lightsail and thus is not likely of direct relevance to the problem under consideration in the present study.
Altogether, while some studies have inquired into the structural stresses and strains supported by solar sails, it should be noted that solar-driven sails are generally not designed to tolerate high photon pressures and resultingly most solar sails envisioned so far-like the IKAROS, NanoSail-D2, and LightSail 2 solar sails [41][42][43][44]-sustain accelerations orders of magnitude less than the accelerations that, for example, the Breakthrough Starshot laser-driven lightsail would need to withstand [6,7]. Consequently, the structural analysis of solar-driven sails has been chiefly concerned with the influence of lightsail deformations upon the spacecraft trajectory whereas, by contrast, the comparatively large accelerations sustained by laser-driven lightsails could cause a deformed lightsail to crumple beyond functional use. Thus far, no inquiry into the detailed vibrations and deformations of a lightsail under high photon radiation loading has been completed. In particular, no study inquired into whether an ideally thin lightsail is capable of sus-taining large laser-driven accelerations despite the inevitable presence of perturbations (see Fig. 1). These perturbationsarising from multiple possible sources such as atmospheric disturbances, beam non-uniformities, shape distortion of the sail, etc.-prevent perfectly uniform loading of the sail, and this complication of the laser-lightsail dynamics may cause the distortions in sail shape to grow in amplitude. Even if the perturbations are, through careful lightsail deployment, kept small in magnitude, the following question remains: Under the large photonic pressure loads required for feasible interstellar flight, will the structural perturbations of the lightsail grow out of bound or will the lightsail remain flat?
The present study addresses the question of lightsail shape stability with perturbations by considering a first-principles approach to the problem. First, a continuous model of a sinusoidally perturbed lightsail under radiation load allows an analytical, quasi-static analysis of the critical point between lightsail shape stability and instability in Section 2. A Lagrangianbased finite element (FE) numerical model of the perturbed lightsail is then constructed using rigid rectangular slices to simulate the full dynamics of the lightsail while undergoing acceleration in Section 3. The rigid-element numerical model was further generalized to include torsion and rectilinear springs to investigate the influence upon lightsail structural stability of material bending stiffness and applied tension, respectively. In Section 4, the quasi-static derived analytical expressions are compared to the numerical results, which consisted of multiple lightsail dynamics simulations with each simulation varying lightsail geometry, material modulus, and/or applied tension magnitude. The engineering implications of the analysis are then explored.

Theoretical Considerations
As a preliminary consideration of the problem of a lightsail under radiation loading, an L × W × h continuous elastic plate model was first constructed that would allow for a static analysis of the criticality between a structurally stable and unstable lightsail (see Fig. 1). The 3-dimensional plate model was then simplified to a 2-dimensional beam model by setting the lightsail width, W, equal to unit length (in meters). The continuous elastic beam model is shown in Fig. 2. The following analysis considers introducing a sinusoidal perturbation into the lightsail shape. The restoring bending moment caused by this deformation is then compared to the moment induced by the incident radiation interacting with the curve surface of the lightsail. If the radiation-induced moment exceeds the bending moment associated with the imposed perturbation and acts in the same direction as the perturbation, then presumably the lightsail continues to further deform. If the radiation-induced moment is less than the restoring bending moment of the perturbation, then the sail would be expected to return toward its original, flat configuration. By equating the bending moment of the sinusoidal perturbation with the radiation-induced moment, a critical condition for lightsail stability can be defined in terms of the radiation intensity, the elastic modulus and dimensions of the lightsail, and the amplitude and wave number of the perturbation. The formalism of this approach follows here.
The beam model analysis is here conducted in a non-inertial reference frame accelerating at g 0 , the lightsail's vertical acceleration, with a body force term ρ h W g 0 acting in the ydirection on the mass elements of the lightsail in accordance with D'Alembert's principle, thereby reducing the dynamic problem to a quasi-static problem. 3 The presence of radiation pressure is modeled by the distributed loading, p r . The expression of p r can be derived by first considering the kinematics and dynamics of an infinitesimal flat sail element of length ds as shown in Fig. 3 for reference. For a lightsail element made of a material with reflection coefficient r, absorption coefficient a, and transmission coefficientτ, the total force imparted by radiation of oblique incidence angle θ can be resolved into the normal-tangential components where s stands for the fraction of light that is specularly reflected; ε f and ε b stand for the emissivity of the front (reflecting) and back (non-reflecting) surfaces of the sail, respectively; and B f and B b are the coefficients accounting for the potentially non-Lambertian nature of the sail element surfaces. Together, these two forces generate a resultant whose direction, m, is skewed away from the element's normal: This optical lightsail model was first proposed by Forward [45] and then was further discussed by Wright and McInnes [46,47] now termed the Forward-Wright-McInnes or FWM model. 4 For the purposes of this paper, the lightsail element will be assumed to be perfectly reflective (r = 1 =⇒ a =τ = 0) with all beam reflections being specular (s = 1). The resulting force per element then becomes a force that is entirely normal to the sail element. Normalizing (4) with respect to the element area generates the radiation pressure: 3 Given the large magnitude of the vertical accelerations, it is here assumed that g 0 =ÿ ẍ,θ(s), and thus the horizontal and rotational inertia of the lightsail were deemed negligible and their d'Alembert equivalent was not included into the lightsail beam model. 4 The authors are aware that more accurate lightsail optical models using vector theories such as the Rayleigh-Rice theory have been proposed [48,49], but the FWM model was used here to ensure a straightforward simplification to the case of an ideal reflector. The implementation of physically more accurate optical models can be incorporated into future studies. To model the presence of surface defects, the lightsail was initially sinusoidally deformed with a deformation profile given by where ν is the mode number of the perturbation, which was taken as whole and halved positive integers, i.e., ν = 1 2 , 1, 3 2 , 2, . . . The acceleration of the sinusoidal lightsail can be found by integrating the forces acting on the sail in the y-direction: where θ is the local angle of the lightsail. Using the small angle approximations sin(θ) ≈ θ ≈ dw dx and cos(θ) ≈ 1 − θ 2 2 , 5 the acceleration of the lightsail can be solved for The term in the brackets represents the second order correction or decrement of lightsail acceleration due to the momentum of photons being scattered off the perturbed reflective surface. Of note here is that as the amplitude of the perturbations, a 0 , goes to zero, that is, as the lightsail flattens, the vertical acceleration becomes the acceleration of a flat lightsail (g 0 ≈ 2I0 c ρ h ≈ 2I0 c m ) [5]. This same result for g 0 can be obtained if only a first order small angle approximation is employed by letting cos(θ) ≈ 1 instead, making equation (8) then yield g 0 = 2I0 c ρ h . To keep all subsequent analytical expressions tractable and conceptually insightful, for the rest of this section, only a first-order approximation will be employed.
To investigate the material strength and/or boundary tension that is required to keep the deformed lightsail in equilibrium, a standard method of sections was performed by conceptually cutting the lightsail at a given point O (see Fig. 2b). Summing the moments acting at that point then gives: The tension T is tangent to the L end of the lightsail, however, given the small amplitude of the sinusoidal deformation, the dynamic contribution of its vertical component is neglected. The 5 The second order small angle approximation for the secant function here used is sec(θ) ≈ 1 + θ 2 2 .
first two integral terms in equation (10) represent the moment contributions of the vertical and horizontal components of the radiation pressure loads, respectively. The third integral term containing ρ h g 0 in the integrand represents the moment arising from the d'Alembert body force, M * . Again using the small angle approximations and retaining only first order terms, equation (10) yields Substituting the first order result for the lightsail acceleration, g 0 = 2I0 c ρ h into equation (11) causes the moments generated by the vertical component of the radiation pressure loading and by the d'Alembert body force to cancel each other to first order. Equation (11) thus reduces to To carry out the remaining integral, recall that the lightsail deformation is defined by equation (6). Substituting the appropriate expressions for w, w 0 , and dw dx , the integration can be carried out to yield −M(x 0 ) − I 0 c W a 2 0 sin 2 2 π ν x 0 L + T a 0 sin 2 π ν x 0 L = 0. (13) Equation (13) defines a state of static equilibrium of the lightsail: The net moment applied to the sail by the acceleration body-force term and the radiation pressure is balanced by a material moment resulting from the bending of the lightsail. The elastic modulus necessary to provide the moment M 0 (x 0 ) can be found from the well known moment-curvature relation where E is the material elastic modulus (Young's modulus), I is the second moment of area about the z-axis of the sail crosssection, 6 and κ is the curvature of the deformed lightsail. Using the small angle approximation for the curvature, κ ≈ d 2 w dx 2 , the material strength required for equilibrium can be solved for: By comparing the applied moments due to radiation and acceleration (in the non-inertial frame) and the moments due to the bending of the stiff sail and due to the applied tension, criteria for the structural stability of the lightsail can be obtained. For example, if the moments due to the lightsail's bending stiffness and tension exceed the moments applied by radiation and acceleration, it is plausible that the sail remains stable. If the moments caused by radiation pressure exceed the bending moment of the lightsail, however, the lightsail is likely to undergo further-potentially catastrophic-deformation. From these considerations, it is possible to define a critical value of elastic modulus in the absence of tension, that is, by letting T = 0: where the extrema values, E crmax , occur at the peaks and troughs of the sinusoid. For the case of tension (in the absence of bending stiffness), the critical value of tension (per unit width of the lightsail) can be found The results of this quasi-static analysis provide a candidate for the functional form of the relation between material properties (i.e., elastic modulus), perturbation amplitude and mode number, and the radiation intensity. The results of equation (18) also suggest a critical value of the tension necessary to prevent the development of instability. Appendix A contains a derivation of the functional expression of the critical material modulus based on a different, energy-based approach. These findings will help guide the computational simulations of the full lightsail dynamics found in the next section of this paper.

Rigid-Element Lightsail Model
With the intent of achieving the rapid transit of a smallscale spacecraft equipped with a 1 meter sail to Mars (1 AU) in 30 minutes or to 0.3 c for interstellar flight within three minutes, the mission at launch would require the use of a 0.1 − 100 GW laser load [5,6]. The dynamic stability of the lightsail at launch was analyzed by first considering a rectangular sail with dimensions L × W × h where L = W = 1 m and h is the sail's thickness. The sail was then divided into n equal slices of length l = L/n and mass m, each attached together via frictionless hinges (see Figs. 4 and 5). The system was further subjected to a laser beam of uniform intensity distribution. Given the idealized flat surface of each slice, specular reflection was assumed to occur between an individual photon and a sail element. Because the reflection was specular and complete, the force resulting from the interaction between the laser beam and the lightsail was normal to the sail element at the point of photon incidence. This uniform pressure thus results in the presence of concentrated forces along the center of mass (CoM) of each individual slice (recall Fig. 3b). In the limit where the number of rigid slice elements is very large (i.e., when n → ∞) this assumption holds true even for non-uniform laser intensities because of the infinitesimal length of each element. The force imparted to the ith element by the photon pressure is thus where A = W × l is the slice element area and θ i is the ith element's inclination with respect to the horizontal.
To investigate the dynamical behavior of the system, the Lagrangian formalism of mechanics was used. First, the Lagrangian of the system was constructed alongside its constraints (note that I G = 1 12 ml 2 is the moment of inertia of the sail element about its center of mass).

Torsion Lightsail Model
To include the material bending stiffness of the lightsail within the study, the finite-element model was further generalized by attaching torsional spring elements with stiffness constant k t = EI/l to the hinges that connect the rigid slices together (see Fig. 6). To see a full derivation of this discrete torsional spring constant from the moment-curvature relation, the reader is referred to [50]. The quantity EI stands for the flexural rigidity of the individual slices with E being the sail material elastic modulus and I being the second moment of area of each rigid element: Employing the same variational approach as for the rigid lightsail model, the Lagrangian of the system is modified with the appearance of a potential energy term arising from the presence of the linear torsion springs: Coupling equation ( θ 1 equation of motion for the torsion model −6(n − 1)l sin θ 1ẍ1 + 6 (n − 1) l cos θ 1ÿ1 θ k equations of motion for the torsion model where k = 2, 3, 4, . . . , n − 1 sin θ k sin θ j cos 2 θ j + cos θ k cos 3 θ j ; θ n equation of motion for the torsion model with the same auxiliary coefficients a xy j , a θ1 j , and a θk k j as in the rigid case. Note that the only difference between the equations of motion of the torsion and rigid models are the addition of the torsional spring restoring moments boxed in the RHS of equations (35) to (37). For the sake of completion and clarity, however, the entire set of ODEs is given. The Supplementary Material accompanying this paper contains a section describing how the torsion lightsail model was validated against classical Euler-Bernoulli beam theory.  Figure 6: Finite-element model generalization using torsional spring elements to simulate material bending stiffness.

Tension and Torsion (TnT) Lightsail Model
To include both bending stiffness and tensile, linear stiffness of the material, the lightsail FE model was further generalized by the addition of rectilinear springs with stiffness constant k s = (n − 1) E h W L and initial/rest length l s . Further, to simulate a pre-tensioned sail, boundary tension was applied at both ends of the discretized lightsail. This new model is depicted in Fig. 7. The reader is referred to the Supplementary Material for a validation of this model against linear wave theory.
The formulation of this new model is now considerably increased in complexity since, to keep track of the rectilinear spring deformation, an additional set of generalized coordinates z = [z 1 , . . . , z n−1 ] T needs to be included along with the previous set of coordinates x 1 , y 1 , and θ = [θ 1 , . . . , θ n ] T . This makes for a full vector of generalized coordinates q TnT = [x 1 , y 1 , θ, z] T , which contains 2n + 1 equations instead of the previous n + 2. Given the presence of rectilinear springs, the Lagrangian of the system was also modified with another potential energy term: The constraint equations were also revised; to relate the coordinates of each element's center of mass to the first element's coordinate (x 1 , y 1 ), the length of the springs, (l s + z k ), now needed to be accounted for in addition to the length of each rigid element: Given the inclusion of boundary tension, the generalized force components were also modified: where r Tn = x n + l 2 cos θ n , y n + l 2 sin θ n , 0 , -y 1 equation of motion for the TnT model -z k equations of motion for the TnT model where k = 1, 2, 3, . . . , n − 1 θ 1 equation of motion for the TnT model sin θ 1 sin θ j cos 2 θ j + cos θ 1 cos 3 θ j θ k equations of motion for the TnT model where k = 2, 3, 4, . . . , n − 1 sin θ k sin θ j cos 2 θ j + cos θ k cos 3 θ j θ n equation of motion for the TnT model with the new auxiliary coefficients

Numerical Methods
Before pressing on with the numerical implementation of each model, note that each set of ODEs can be rewritten in the more compact matrix form where the subscript "m" attached to the generalised coordinate vectors in the above denotes the respective model being studied be it the rigid, torsion, or TnT lightsail model. As per numerical fashion, the second-order system described by equation (48) can readily be turned into a first-order system through the introduction of an auxiliary vector of unknowns u =q m . The above system can thus be rewritten in the more numerically friendly form wherê Given the dependency of the modified mass matrixM upon the numerical vector of unknowns x, the derivative of the unknown cannot be isolated in equation (49). This renders the use of more conventional ODE solvers inadequate, and thus, to solve the numerical system of equations at hand, the set of firstorder ODEs was instead treated as a set of differential-algebraic equations (DAEs). Opting for solvers better adapted to treat such DAEs, Mathematica's NDSolve package was employed with the enforced residual method of simplification whereby the NDSolve function first rewrites (49) into the fully implicit form As its name implies, this MEBDFV solver employs a modified version of the backward differentiation formulae whereby, if x i is the solution of the unknown at the current time step, a superfuture value, x i+1 is computed to improve computational stability by refining the initial guess of the derivative of the unknown, x i , thereby allowing for a more accurate Newton-Raphson solution to the nonlinear system (50). For a concise description of the overall MEBDFV scheme, the reader is directed to [52]. The MEBDFV solver also includes additional intricacies such as adaptive time stepping and Newton-Raphson scheme convergence rate estimates, but their detailed description goes beyond the scope of this present paper. For a thorough treatise of the inner workings of the MEBDFV code, please refer to the papers written by Cash et al. on the subject [53][54][55][56]. A Mathematica vs FORTRAN lightsail simulation output comparison test case can be found in the Supplementary Material.
While the authors recognize that a commercial FE software would have been able to numerically address this problem, the main goal of this study was to develop an analytic theory for lightsail stability under large radiation pressure. As such, it is important that the numerical method used to validate the theory should be derived from the same underlying physical model, and this would only be possible in a finite element model developed from the ground up. Further, the development of a FE model from first-principles allowed the authors full knowledge of the assumptions and limitations of each lightsail numerical model, which in turn granted the ability to isolate the influence of various effects upon lightsail shape stability such as bending stiffness, boundary tension, etc. The fact that bending stiffness was incorporated in this study is of note since conventional thin film theory usually dispenses with the resistance to bending de- formation altogether. Finally, this first two-dimensional, rectangular model here studied both analytically and numerically was also chosen for inquiry as it could potentially act in future studies as an elementary building block for the construction of more general (i.e., non-rectangular) perturbed lightsails by coupling together, for example, multiple torsion or TnT lightsail models.

Torsion Model Results
Given the rather large number of parameters required to perform a torsion model simulation, the theoretical expression (17) for the critical lightsail material modulus E crmax was used as a guideline to reduce the parameter space to a more manageable size. As a result, only the material modulus, E, and the lightsail thickness, h, were varied during the torsion model simulations. In particular, h was varied from 0.01 µm to 10 µm in jumps of powers of ten while the material modulus was varied over a few magnitudes from its theoretical critical value. Each simulation was run for a total of 1 s and failure was considered whenever the perturbation amplitude doubled from its initial value. That is, at each numerical time increment, the maximum amplitude of the sail, amp max was computed by considering the vertical position of the CoM of all sail elements: Whenever amp max (t) ≥ 2a 0 , the simulation was stopped and the time to failure, τ, was set to τ = t. If the amplitude of perturbations never doubled during the simulation, then the time to failure was set to the total runtime, τ = t final . In total, for the torsion model, 900 lightsail simulations per mode number for each of modes ν = 1 and ν = 3/2 were run using 50 elements by varying the thickness, h, from 0.01 µm to 10 µm and the elastic modulus, E, from 4.5 × 10 −3 GPa to 4.5 × 10 7 GPa with increments equally spaced on a logarithmic scale. Each data point was a 4-tuple recording the lightsail simulation thickness, elastic modulus, initial perturbation mode number, and time to failure (h, E, ν, τ). Using 50 elements per simulation appeared justifiable because, despite the fact that asymptotic numerical convergence starts at about 100 sail elements (see Supplementary Material), the qualitative structural behavior of the sail was shown through convergence tests to remain the same independent of the number of elements used (i.e., a stable sail remained stable and an unstable sail remained unstable independent of the number of sail elements used; see Supplementary Material). All other simulation parameters were kept constant from simulation to simulation and their values are given below:  sail torsion model run where E crmax = 45 GPa. Given the quasiperiodic nature of the maximum amplitude value, the lightsail remained stable, that is, its initial perturbation amplitude never doubled during the runtime. Figures 8a-d shows four snapshots of the lightsail's shape each recorded at different simulation times, which visually further reinforce the fact that the lightsail shape dynamics for this particular set of torsion model pa-rameters is quasi-periodic. 7 Such lightsail configurations were termed stable oscillatory. Figure 9e, on the other hand, depicts the time evolution of the maximum amplitude for a torsion model run where h = 10 µm, E = 0.108 GPa, ν = 3/2, τ = 0.541 s . Notice that the sail perturbation amplitude grew without bound 7 In order to help better appreciate the lightsail shape dynamics, all numerical lightsail snapshot figures found in this paper were visualized in a reference frame fixed to the body of an ideally flat lightsail accelerating at g 0 = 2I 0 c ρ h . as indicated by the exponential behavior of the maximum amplitude, which effectively blew up as we approached the end of the simulation. Snapshots of the lightsail shape were taken at different times for this simulation as well and of note are both the absence of periodic behavior and the fact that the lightsail appears to have failed at a trough (see Figs. 9a-d). Simulations such as these, which did not exhibit stable oscillatory behavior but whose initial perturbation amplitude doubled some time before the runtime ended, were termed potentially unstable.
A summary of 1800 simulations (900 per perturbation mode) in the form of stability plots is found in Fig. 10 where Fig. 10a and Fig. 10b contain the data points for perturbation mode 1 and 3/2, respectively. Onto each plot was added the theoretical stability threshold predicted by the analytical equation (17). Of most remarkable note is that, for both perturbation modes, nearly all stable oscillatory sample points (i.e., all points with a τ value of 1 s) fall above the theoretical curve and vice versa for potentially unstable sample points where, the further away from the theoretical curve the unstable data point was, the more unstable the simulation was (i.e., the faster the doubling of initial perturbation amplitude occurred). This distribution of stable/unstable data points indicates agreement between theory and numerical simulations: if the elastic modulus of a given sample run stood a few orders of magnitude above the critical value predicted by theory, that particular lightsail configuration behaved in a quasi-periodic, stable manner, otherwise the configuration displayed unstable behavior. Of further note is the fact that qualitative and quantitative differences between the two initial perturbation modes appear relatively minor, indicating that the agreement between theory and computation is independent of initial perturbation mode.

Tensile Model Results
As with the torsion model, the expression (19) for T crmax as derived from theory was used to make the parameter space of the TnT model more tractable. Consequently, only the boundary tension, T , and the lightsail thickness, h, were varied. A total of 900 lightsail runs for each of modes ν = 1 and ν = 3/2 were performed using the TnT model by varying the lightsail thickness, h, from 0.01 µm to 10 µm and the boundary tension magnitude, T , from 3.34 × 10 −7 N/m to 3.34 N/m. 8 Each simulation or data point once more recorded the lightsail thickness, boundary tension magnitude, initial perturbation mode, and time to failure (h, T, ν, τ). Note also that, in order to solely investigate the influence of membrane action or tensile stiffness upon the lightsail structural stability, the bending stiffness of the lightsail model was manually tuned down to 0 by setting k t = 0 (although, because of the dependency of bending stiffness upon 1/h 3 , the presence or absence of torsion springs had no significant influence given the range of the thickness values studied). Once more, each simulation was run for a total of 1 second using 50 elements and failure was considered to have occurred when the light-sail perturbation amplitude doubled. The simulation properties were as listed below: As in the bending stiffness case, the simulation of a lightsail whose boundary tension magnitude was at least an order of magnitude above that of its critical tension value displayed stable oscillatory motion. Conversely, decreasing the boundary tension value by an order of magnitude or more below the critical value always made the lightsail model display potentially unstable behavior. See, for example, the h = 10 µm, T = 3.09 × 10 −6 N/m, ν = 3/2, τ = 0.494 s lightsail TnT configuration whose uncontrolled perturbation growth is shown in Fig. 12. Figure 13 depicts the stability maps for perturbation modes 1 and 3/2 for the TnT model simulations. As with the torsion model, plotting the theoretical stability curve onto the stability maps indicates agreement between theory and simulations: nearly all stable data points lie above the theoretical stability curve and almost all unstable data points lie below the theoretical curve, with the data points farthest below the theoretical curve being more unstable than the data points closest to the theoretical curve. Noteworthy is also the fact that, once more, there appears to be little qualitative and quantitative difference between the different perturbation mode stability maps.

Engineering Implications
The results encapsulated in the stability maps (Figs. 10 and 13) showed good agreement between theory and numerical computations. It should be noted, however, that stability is not strictly guaranteed when either the material modulus, E, or the boundary tension magnitude, T , are greater than their theoretical critical maximum value; typically, it is necessary to exceed the critical values by an order of magnitude to ensure that the lightsail shape remains stable. The inclusion of a safety factor in equations (17) and (19) should, in principle allow for the construction of a structurally stable lightsail configuration.
Upon closer inspection of the torsion model stability maps depicted in Fig. 10, it should be remarked that the material moduli required to generate a stable sail configuration are impractically large (E crmax = 45 000 GPa for a 1-µm-thick sail with a mode 3/2 and 0.01 mm amplitude perturbation under a laser flux intensity of 10 GW/m 2 , compared to the elastic modulus of diamond, E ≈ 1000 GPa). This is believed to be in large part due to the dependence of the required material strength on the inverse of the thickness cubed, E ∝ 1/h 3 , a dependence emerging from moment-curvature relation (14), where the second moment of area of the lightsail is I = Wh 3 /12. Consequently, if the lightsail is to be kept thin for practical accelerations to be reached, the inclusion of boundary tension into the design of the lightsail would be necessary in order to achieve structural stability since the boundary tension thresholds for stability appear feasible. For example, for the same 1-µm-thick sail with a 0.01 mm amplitude perturbation, T crmax = 0.334 mN per unit width for a laser flux intensity of 10 GW/m 2 , which can be achieved by inflating the lightsail (much akin to a balloon) using a gas pressure of 1.4 mPa or by spinning the lightsail about the laser beam axis at a rate of 1.6 rad/s-although, in the spinning case, the critical tension would be achieved only at the periphery of the lightsail. Ultimately, how the boundary tension is added into the lightsail concept is left at the discretion of the design engineer, and other ideas like the Cassenti and Cassenti proposed boundary ring [34] and the inverse cat eye lightsail optical metasurface concept of Siegel et al. [21], etc. have been proposed in the literature. Not only do the boundary tension magnitudes required for structural stability appear readily achievable, but of important note is the absence of material strength and lightsail thickness in the functional expression of T crmax . The use of boundary tension would then effectively leave the design space of material strength and lightsail thickness entirely open and thus allows the conceptual engineer to, for example, base their choice of material on optical properties alone without having to worry about structural stability implications.
Further, a closer look at equations (17) and (19) reveals that both the critical material modulus and boundary tension magnitude are directly proportional to the perturbation amplitude, a 0 .
Given this direct correspondence between the critical stability values and the amplitude of the perturbations, it is an imperative that good quality control be undertaken during the manufacturing and packaging/handling processes of the lightsail. As already discussed, the previously quoted T crmax = 0.334 mN per unit width was taken for a perturbation amplitude of 0.01 mm. If the amplitude of perturbation is increased such that a 0 = 10 cm instead while maintaining all other parameters constant, the critical maximum value becomes T crmax = 3.34 N per unit width. We see here a substantial increase in the threshold value and the reader is reminded that, as per numerical simulations, stability is only guaranteed when the boundary tension or material modulus are at least an order of magnitude above their respective thresholds.

Lightsail Model Limitations
It remains of note that the analytical and numerical models and results studied thus far hold certain limitations. The laser beam intensity was here taken as uniform, and the nonuniformity of a more realistic laser beam intensity distribution like that of a Gaussian beam may drive further instability. Further, and perhaps of greater importance, is the fact that no lightsail model included the presence of dissipation. This is especially important given the stable oscillatory behavior of the torsion and TnT models when above their respective E crmax and T crmax values. In particular, it was noted during the generation of the TnT model stability maps that further increasing the boundary tension beyond its critical value increased the frequency of oscillations of the lightsail. Increasing frequency of oscillation in the presence of dissipation effects can potentially increase  the internal heat generation to a point where the lightsail may melt, even in the absence of any radiation absorption. Finally, each sail element was modeled as perfectly reflective. As noted in the introductory section, the material optical properties can potentially be used to help enforce lightsail stability. In order to improve the physical accuracy of the lightsail model it will be thus necessary in the future to include: a realistic, non-uniform laser beam intensity distribution, heat dissipation effects, and material optical properties.

Conclusion
The large directed-energy loads required for laser-driven propulsion pose the problem of the stability of the lightsail shape. Given the unavoidable presence of perturbations along the reflective sail surface, it is important to determine whether the sail will retain its shape or crumple under large photon pressures. An analytical theory of stability was first proposed wherein critical values for the material modulus, E cr , and for boundary tension, T cr , were derived using a moment balance, involving key lightsail parameters such as thickness and laser intensity. The analytical model was compared against various numerical simulations of the dynamics of an accelerating lightsail, and agreement was found between theory and numerical computations. In the absence of boundary tension, both theory and numerical simulations predict the necessity of unrealistically large elastic moduli to guarantee the lightsail's stability. Increasing the sail thickness can potentially generate stable sail configurations, but to the detriment of spacecraft acceleration. With the addition of boundary tension, it was found from both theory and numerical simulations that the lightsail can be kept under a stable oscillatory state through the application of relatively small tensile loads, even for very thin (sub-micron) lightsails. The tension required to maintain a stable shape could be readily provided via a sail inflation, spinning of the sail, or other tensioning mechanics. The parameter maps provided by the analytical theory and verified through numerical simulation could potentially provide the mission designer with a convenient method to engineer a lightsail capable of sustaining the large directed-laser pressures needed for high-acceleration spacecraft applications.
Integrating the above result while employing the moment-curvature relation and a first order small angle approximation, the total internal strain energy becomes 3) The first variation of the lightsail's internal bending energy can thus be found by applying the functional derivative: After performing integration by parts twice, δU can be expressed only in terms of the first variation of the vertical deformation 9 In the absence of tension, the virtual work applied to the lightsail is composed of the vertical and horizontal components of the radiation pressure loads and of the d'Alembert body force, Once again employing small angle approximations, the vertical radiation pressure and d'Alembert body force virtual work contributions cancel out and the total virtual work simplifies to The virtual displacement in the horizontal direction, δu, present in equation (A.7) can be related to the virtual displacement in the vertical direction, δw, by assuming the continuous beam lightsail model to be inextensible. Inextensibility can rightly be assumed in the case where the tension running through the sail before the application of loads is of such magnitude that 9 The boundary term E I d 2 w dx 2 δ dw dx L 0 is neglected because the curvature of the lightsail is zero at the endpoints by construction. any additional deformation will not induce the change of the body's total length. The inextensibility assumption can also be applicable in the case where the tensile stiffness of the lightsail is altogether ignored with the investigation efforts being solely focused on the bending stiffness of the lightsail as is presently done. Assuming an inextensible lightsail model, it can be shown that a standard result found in the literature [57][58][59]. From the above kinematic constraint, finding the first variation of u(x) becomes a straightforward application of the functional derivative followed by an integration-by-parts: where the partial derivative symbols have been replaced by a total derivative because the vertical displacement here spatially and temporally depends by construction only on ξ, which is here a dummy variable of integration for x. The boundary term in equation (A.9) can be more explicitly written to yield with the second boundary term, dw dx δw x=0 , not being dependent upon x. Substituting equation (A.10) into equation (A.7) and expanding the virtual work expression gives The second integrals can be solved for to yield 0 given that the term in brackets, dw dx δw x=0 , effectively acts as a constant in the integrand. The order of integration of the double spatial integral in equation (A.11) can be interchanged by properly accounting for the limits of integration and by noting that the dummy variables of integration x and ξ can themselves also be interchanged: Having found an expression for the virtual work depending only on the vertical virtual displacement, δw, the principle of virtual work may finally be used to investigate the question of structural stability. Substituting equations (A.3) and (A.13) into equation (A.1) results in the following: (A.14) Given the arbitrary nature of δw, it is noted that for the second term in equation (A.14) to be equal to zero, the entire quantity found inside the braces in the integrand must be equal to zero. Thus, 15) or, 16) is noted to look the same as the familiar Euler-Bernoulli (EB) beam theory result, E I d 4 w dx 4 = q(x) where q(x) stands for the distributed loading applied onto the EB beam. The radiation pressure loading and d'Alembert body forces applied onto the lightsail are also noted to generate a RHS in equation (A.16) whose dependency on the lightsail slope and curvature renders the differential equation nonlinear. Recalling that the initial deformation of the lightsail is, by construction, w = a 0 sin 2 π ν x L , the critical material elastic modulus can be solved for in terms of the lightsail parameters by substituting in the appropriate expressions for w, dw dx , d 2 w dx 2 , and d 4 w dx 4 into equation (A.16): Comparing equations (17) and (A.18) shows that the direct approach and the energy approach predict near identical results as to the critical material strength required for equilibrium: These two equations only differ by a factor of 1/2. 1 and mode 3/2 stability maps previously encountered in Section 4. After laying out both of the theoretical curves onto the numerical maps, there appears to be little difference between the direct method and the energy method in terms of qualitatively predicting the lightsail stability boundary: no matter the analytical expression chosen, if the material strength of the lightsail is an order of magnitude or so above that of the theoretical maximum critical value, then the lightsail configuration ends up being stable and vice versa. The factor of 1/2 difference between equations (17) and (A.18) is ultimately inconsequential since, as previously remarked in the Engineering Implications Section 4.3, both equations require the addition of a safety factor to ensure that the lightsail configuration is stable as stability is guaranteed only when the material modulus of the sail is an order of magnitude greater than either E crmax or E crmax δW . This strong similarity between the two analytical approaches helps further verify the agreement between theory and numerical computations found in this paper with respect to the lightsail structural stability problem. This Supplementary Material documents the methods used to validate the various lightsail numerical models. In particular, this material contains the Euler-Bernoulli beam validation case for the torsion lightsail model (Section 1) and the harmonic string validation case for the tension and torsion (TnT) model (Section 2). A study of the convergence of both lightsail models is described in Section 3. Section 4 contains a Mathematica vs FORTRAN lightsail simulation output comparison test case.

Euler-Bernoulli Beam Validation
Given that the torsion lightsail model was constructed specifically to account for material bending stiffness in the elastic regime, a relevant validation case here is that of the classic Euler-Bernoulli cantilever beam. The beam, clamped at its left end, sustains an applied concentrated load P at its right end (see Fig. 1). The Euler-Bernoulli theory predicts, for small deformations, the deflection The setup for the numerical simulation is similar to the torsion lightsail model with only a few minor exceptions to account for the applied forces and boundary conditions. Since the left end is clamped, the coordinates of the center of mass of the first element is set to (x 1 = l 2 , y 1 = 0) and the angular position of the first element is held constant as θ 1 (t) = 0. This, in turn, reduces the size of the system of equations from n + 2 to n − 1. Finally, the forces imparted to each sail element by the photon pressure f i were set to 0 and replaced by a static concentrated force, P = (0, −P 0 , 0) applied at the center of mass of the last element (x n , y n ). This new FE model can be seen in Fig. 2. δ(x) P Figure 1: The analytical Euler-Bernoulli model.  To properly simulate static loading while keeping the computational demands of the validation to a minimum, it was decided not to make the transverse load a function of time. Instead, the applied load was kept as a simple step input while the numerical model was modified to include (linear) angular dashpots at its hinges to help dampen the response of the system following the application of the load. The addition of these angular dashpots further modifies the numerical model's EOMs ever so slightly due to the presence of dissipative moments. These dampening moments can be modeled using the Rayleigh dissipation function The generalized dissipative forces that need to be accounted for are whose final form is In the absence of such dissipative forces the numerical model would output periodic to quasi-periodic behavior given the step load input and would never settle to a static deformation curve within computationally reasonable time. To find the magnitude of the dissipative forces required to ensure a static numerical response, the material and geometrical simulation parameters first had to be set. The simulation parameters for the beam validation are: the static value as seen in Fig. 3. Ideal here means that reducing the dampening constant value by more than one order of magnitude will cause the numerical system to overshoot its static curve before settling, that is, reducing the value of the ideal c d = 10 7 N m s by more than an order of magnitude will cause the numerical system to respond like an underdamped system.  To validate the torsion model, the Euler-Bernoulli beam simulations were run for different numbers of elements and the maximum error between the analytical deflection, δ(x), and the numerical beam deflection at end of runtime was computed for each run by taking the absolute difference between a given element's vertical coordinate and its theoretical deflection: e EBmax = max δ(x j (t final )) − y j (t final ) for i = 1, 2 , 3, ..., n. Figure 4 plots the maximum absolute error vs the number of elements. The number of elements was increased from n = 5 to n = 215 in increments of 5. Note that, after about n = 150, the absolute maximum error seems to start plateauing with the final maximum error at n = 215 being e EBmax = 0.002 57 mm. Figure 5 shows the deflection curves of both the numerical and analytical model for the n = 215 simulation run where it is noted that the two curves are nearly identical. This relatively small error between theory and numerical computations instills confidence into the ability of the torsion model to correctly describe material bending stiffness.

Harmonic String Validation
Given its inclusion of rectilinear springs and boundary tension, a pertinent scenario to validate the TnT model for transient problems against is that of a vibrating string's harmonics. From Newtonian first principles, it is known that the fundamental frequency, f 0 , of a string under constant tension, T , and kept fixed at both ends can be expressed as where L is the working length, that is, the stretched length of the string and M is its total mass. To help make the validation as general as can be, as set of geometric and material simulation parameters were first set. From these starting parameters, a theoretical fundamental frequency for the string could be computed using equation (6). This theoretical value was then compared to the fundamental frequency of the simulated spring measured via the use of a discrete Fourier transform (DFT). The simulation parameters chosen for the validation test were the following 3,4 : In adapting the TnT model to the taut string scenario, the following issue is encountered: in order to keep the right end of the TnT model fixed, it is not enough to simply set y n = constant. This holonomic constraint would, in turn, eliminate the connection between the last TnT element and the remainder of the string. It would be as though the hinge connecting the last extensible element to the other elements was deleted at the start of the simulation. Ideally, one would want the vertical acceleration component of the last element to be zero throughout the simulation, that is, ideally,ÿ n (t) = 0. In reality, however,ÿ n = f (q,q), and, in fact,ÿ n is non-linear inq and thus Lagrangian multipliers cannot be used to account for the generalized force(s) needed to maintain theÿ n (t) = 0 constraint. Fortunately, maintaining the first string element fixed, that is, setting (x 1 = constant, y 1 = constant) and ensuring that the two applied boundary tension forces remain co-linear to their respective element, that is, enforcing T 1 = T (− cos θ 1 , − sin θ 1 , 0) and T n = T (cos θ n , sin θ n , 0) appears to keep the end element from straying too much from its starting position.
Because the objective of this procedure was to validate the fundamental frequency of the simulated string against that of a physical string, the numerical string was initially deformed in a half sine wave fashion where y(x) = a 0 sin π x L with a 0 = 10 mm. An initial amplitude of 10 mm was used instead of the 0.01 mm used throughout the main paper for the lightsail stability maps in order to ensure that the initial mode shape of the simulated spring was of significant scale as opposed to a perturbation. Further, to simulate pre-tension, the initial elongation of each rectilinear spring was set such that they each exerted a total force equal in magnitude to the boundary tension, that is, initially z i = T/k s for i = 1, 2, 3 ..., n − 1. Figures 6 and 7 display the validation results which were gathered by increasing the number of elements used to simulate the string from n = 5 to n = 135 in increments of 5. To properly measure the initially deformed string's fundamental frequency, the frequency of vertical oscillations of each string element with the exception of the ends was measured via a discrete Fourier transform. The mean fundamental frequency of the total elements is plotted in Fig. 6. Note that, as per theory, each element ended up oscillating with the same frequency and thus the standard deviation from the mean is always zero, hence the absence of error bars in Fig. 6. The mean fundamental frequency of the string appears to have settled, from the very start of the convergence test at n = 5, to that of 199.96 Hz which is less than 3 Hz away from the target fundamental frequency of 196.00 Hz. The discrepancy is suspected to be due at least in part to the DFT scheme and perhaps in part also due to the "wobbling" of the nth element (i.e., the final element, which should remain fixed). Figure 7 portrays the maximum amplitude of the vertical oscillations of the nth element. Note that the wobbling of the nth element steadily tends to disappear as the number of elements is increased. Figure 8 plots the convergence data for the torsion model. The torsion model appears to converge to the τ ≈ 0.7 s value following the n = 70 simulation case. The convergence study results for the TnT model are shown in Fig. 9 where it can be seen that, for the chosen lightsail configuration, the time to doubling of initial perturbation amplitude steadily converges to τ ≈ 0.05 s with the asymptotic behavior clearly starting at around n = 100 elements. The lightsail simulation runs performed for the construction of both stability maps found in the Results and Discussion section are noted to have used a lesser number of elements with n being set to 50 elements for the sake of numerical efficiency. Given the convergence results, it may be suspected that the number of elements used in the stability maps -especially for the stability maps employing the TnT model -was too low. However, while increasing the number of elements does affect the results quantitatively, the number of elements appears not to affect the qualitative stability results. That is, a stable lightsail configuration remains stable no matter the number of elements employed and likewise for an unstable lightsail configuration.

Mathematica NDSolve vs FORTRAN MEBDFV Comparison
All results shown in the above two validation sections, in the numerical convergence section, and in the main paper were a consequence of numerically integrating a discrete lightsail model or slight modification thereof using the NDSolve solver package of Mathematica. To further gain confidence into the validity of the numerical models, a number of lightsail configurations where also simulated using the FORTRAN modified backward differentiation formulae (MEBDFV) solver of Abdulla and Cash from Imperial College, London (Department of Mathematics), a FORTRAN solver written to deal with the kind of numerical systems encountered in this study (see Section 3.4 of the main paper for a more complete discussion of the inner workings of the MEBDFV solver). This section of the Supplementary Material will focus on the comparison between the Mathematica NDSolve and FORTRAN MEBDFV outputs for a representative lightsail configuration of the mathematically and numerically more complex TnT lightsail model. The TnT lightsail model configuration in question is the h = 10 µm, T = 1.29 × 10 −2 N/m, ν = 3/2, τ = 1 s simulation run whose behavior was visually displayed in a set of temporal snapshots in the main paper. For the sake of completion and convenience, all other parameters relevant to the lightsail configuration studied here are restated below from Section 4.2 of the main paper: 1. Sail length, L = 1 m 2. Sail width, W = 1 m 3. Sail density, ρ = 1000 kg/m 3 4. Material Elastic Modulus, E = 5 GPa 5. Incident laser intensity, I 0 = 10 GW/m 2 6. Initial perturbation amplitude, a 0 = 0.01 mm 7. Perturbation mode number, ν = 3/2 8. Number of elements, n = 50 9. Total runtime, t final = 1 s or to failure Figures 10a-c show the Mathematica lightsail output (red) and the FORTRAN lightsail output (blue) at three different snapshots in time as seen from a reference frame attached to a flat lightsail-the flat lightsail reference frame was used to help better appreciate the perturbed shape of the lightsail. At each instant in time there was no visually appreciable difference in lightsail shape between the Mathematica and FORTRAN outputs with the greatest percentage difference between the two solvers in lightsail shape maximum amplitude being of 2.4% and occurring at the end of the simulation time, t = 1 s. While there was some appreciable difference between the two solver outputs in absolute vertical position of the sail with the Mathematica output being consistently ahead of the FORTRAN output, the greatest difference in vertical position was of 0.025 mm and also occurred at the end of the simulation, that is, at time t = 1 s. Considering the fact that, given the 10 GW laser intensity, both lightsails traveled a total of 3335 m during the 1 s runtime, this difference in absolute position between the two solvers is here considered acceptable. Of note is that both solvers displayed stable oscillatory lightsail behavior, and it is important to emphasize that the Mathematica vs FORTRAN comparison test case shown in this section is a representative of all the comparison tests made between the two solvers: Barring the presence of minor differences in absolute positioning of the lightsail, the evolution of the lightsail shape throughout simulation runtime was identical between the Mathematica NDSolve and the FORTRAN MEBDFV outputs, that is, whenever a Mathematica lightsail simulation displayed stable oscillatory behavior, the FORTRAN lightsail simulation also displayed stable oscillatory behavior, and whenever a Mathematica lightsail simulation displayed potentially unstable behavior so did the FORTRAN lightsail simulation.