1. Introduction

Since the early days of astronautics, light-sails that harness optical forces have been in discussion [1]. Solar-sails using the light of the sun have been at the center of research for a long time and culminated with the launch of the IKAROS satellite in 2010 [2]. However, due to their limited speed, recent research rather revolves around laser-propelled spacecrafts [3–5]. Most notably, the Breakthrough Starshot initiative has set the goal to accelerate gram-scale spacecrafts to a fraction of the speed of light by an array of lasers, shining from the earth to the spacecrafts [6]. A concept art of such a spacecraft is shown in Fig. 1. The initiative has the goal to reach the next solar systems within a matter of decades [7–9].

 

Fig. 1. A concept art for a curved light-propelled metasurface being illuminated and hence driven by an earth-based laser array during the acceleration period. The concept art reflects the physical situation studied in this work. The sail is the curved surface that is shiny in the image. The spacecraft is the object in front of it (not considered in our analysis). The illuminating laser is indicated in faint red. Background image was taken with permission from Anvar Ghaderi. The Drawing of the image was inspired by figures published in [7,10].

Download Full Size | PDF

There are several challenges associated with this ambitious project [10,11]. Foremost, the light-sail has to be lightweight while being simultaneously highly reflecting at the operational wavelength. Moreover, it should provide a very low absorption. The combination of high reflection and low weight is important for efficient propulsion, and the low absorption is required to limit the heating during the acceleration to avoid possible thermal damage. Metamaterials, artificial materials that allow very precise control over the optical properties, are promising candidates that can cope with these restrictions [12]. Metasurfaces comprising particles made from low-absorbing and high-permittivity dielectric materials, e.g. silicon, meet simultaneously the above-mentioned requirements [13–18]. Therefore, metasurface-based light-sails were suggested in the past that provide extremely large propulsion [19].

But there is yet another challenge to be addressed when considering metasurfaces as a light-sail. During the acceleration period, the sail has to maintain a stable orientation relative to the laser. The light-sail has to be oriented normal to the principal direction of the illumination to transfer the linear momentum most efficiently. Moreover, any deviation from the optimal orientation should self-compensate, i.e. the orientation of the light-sail must be stable relative to an arbitrary orientation. If that would not be the case, the sail would automatically and nearly instantaneous get into a tailspin and an efficient acceleration would no longer be possible. While some control could be offered via the phase of the individual lasers forming the array [20], a much better solution is to design the metasurface such that an angular deviation tends to induce a torque of opposite rotational sense such that the sail automatically reorients suitable. That will be a much more appealing passive stabilization mechanism.

As optical metasurfaces promise to provide optical properties on demand, they can accept such a challenge and it is by no means surprising that various designs have been already presented to solve this issue. Most of these prior attempts considered planar bi-gratings that use space-variant metaatoms to achieve this effect of self-stabilization [21–24], or the response of the structure is assumed beforehand [25]. In contrast to these, we follow here an alternative approach and use a rigid curved metasurface with identical metaatoms as a light-sail. Other than in previous research performed on curved light-sails [26–28], which considered a perfectly reflecting bulk material, the parameters describing the position of the scattering elements that form the metasurface are a critical point of our research. Moreover, the interaction of light with the light-sail is studied using full-wave optical simulations of the full structure.

In Section 2, the theoretical background for the stability analysis will be provided, followed by the discussion of the setup in Section 3. Please note, we will consider a finite-sized metasurface based on identical scatterers arranged along a curved surface characterized by some finite radius. As a central part of our work, in Section 4, a fast method to determine rotational stability will be developed. Section 5 contains the results of our stability analysis and we propose an optimally stabilized curved metasurface based on that analysis. In Section 6 we will study the angular trajectory during the acceleration for the designed metasurface for some selected initial conditions. This will be done using full-wave optical simulations coupled to an evolution equation for the dynamics.

For the design of the optimized curved metasurface with realistic constituents, we have chosen c-Si spheres. Crystalline silicon has not only a high refractive index but also a very low absorption and relatively low mass density [10]. Other materials, as pointed out in Ref. [10] can be used as well. However, we deemed c-Si to be most suited, especially because of the very low absorption to refractive index ratio at the operating wavelength of 1000 nm.

Throughout the manuscript, the incident light is considered to be a plane wave. A time dependency of all fields according to $\mathrm {exp}(-2\pi \mathrm {i}\nu t)$ is considered.

2. Theoretical background

2.1 Optical force and torque

When light impinges on a particle, the light is scattered. The exchange of momentum between light and matter leads to the exertion of an optical force and torque on it [29,30]. In our case, the force and torque, and especially the latter one, will be calculated on the actual sail that we can consider in simulation. This optical force and torque can be calculated using the formalism of Maxwell’s stress tensor where it is expressed as a surface integral across the entire object of the time averaged Maxwell stress tensor or its cross-product with some spatial coordinate, respectively [31]. The Maxwell stress tensor depends only on the total field at the considered spatial location, i.e. the scattered plus the incident field in the language of scattering theory [32,33]. Assuming free space as the surrounding medium and assuming that the particle’s motion within one oscillation period of the electromagnetic field is negligible, the time averaged force and torque is equal to:

To compare the strength of the torque with some known fundamental physical limits, we always express it normalized to the maximal optical torque exerted by a plane wave on an isolated magnetic dipole in its resonance. It can be shown that this maximal optical torque is equal to [34]:

2.2 Euler’s equation of rotation

As the incident light is a plane wave, spatial changes will not influence the optical force or torque. Therefore, the rotation and the motion of the center of mass are decoupled. In this paper, we will only consider the rotation as it reflects the stability. For this, Euler’s equation of rotation will be used. The torque is defined as the time-derivative of angular momentum $\mathbf {L}$

The tensor of inertia $\overleftrightarrow {I}$ for an ensemble of $n$ point masses each with a mass $m_n$ is calculated as

3. Considered setup

We will restrict ourselves to rigid metasurfaces made from a finite number of particles. To construct the setup, the particles are initially arranged in a square array with a lattice constant of $\Lambda$ in the x-y-plane. It provides their coordinates $(x_i,y_i,0)$. Afterwards, only particles inside a circumscribing circle [shown in Fig. 2(a) in red] are retained. The radius of that circle dictates the finite size of the considered metasurface. To implement the curvature, the z-component of the position of each metaatom $z_i$ is set to

 

Fig. 2. (a) An exemplary simulation setup for $\Lambda /\lambda =0.75$ and $R=6\lambda$. The scatterers are represented as light blue points. We show them here only as point-like objects to preserve the technical nature of the drawing. In the final device to be studied these would be the actual objects that make the metasurface. The edge of the metasurface is marked with the red line. The radius of curvature is inserted in black. A part of the field sampling grid box is denoted in red, blue, and green small dots. (b) Sketch of the convention of Euler’s angles used. The order of rotations is $\alpha \rightarrow \beta \rightarrow \gamma$, which results in the transformation of the blue configuration to the red one. The sketch shows the axes in the body-fixed frame and the space-fixed system as well as the position of the rotated metasurface in the space-fixed system, here depicted with the blue disc.

Download Full Size | PDF

All simulations are performed with 9 scatterers along the diameter of the circumscribing circle, resulting in $N_{\mathrm {sca}}=49$ scatterers that form the metasurface. For the identification of a favorable curvature and lattice constant of the metasurface, we will initially only take into account hypothetical particles characterized by a magnetic dipole moment in resonance and vanishing contributions from all other multipole moments. At a later stage, actual particles are considered which offer, nevertheless, such a desired optical response.

To describe the rotation of the metasurface, Euler’s angles will be used. We will follow the intrinsic $zxz$-convention as shown in Fig. 2(b). To calculate the optical torque, the total scattered fields are required from the metasurface upon plane wave illumination. The scattered fields are calculated using the multipole expansion method [35,36]. This is done based on the addition theorem that allow describing the scattered fields of the surrounding scatterers as the incoming light for one specific scatterer. With this tool at hand, the field is sampled equidistantly on a cuboid surface enclosing the metasurface as seen in Fig. 2(a). From these points, the integral in Eq. (2) is approximated using Simpson’s method. The sides of this box are all planes aligned to the global $x-$, $y-$, and $z-$axis.

The position of each plane is determined such that the distance to the outermost scatterer is $\lambda$. This value is chosen together with the number of sampling points $n_{\mathrm {samp}}/\lambda =10$ such that the results are converged. At a distance of one wavelength, the near-fields sufficiently decayed and fewer points can be used for converging results.

4. Definition of stability

The key factor for deciding whether a metasurface is stable or not is the azimuthal angle $\beta$. A tilting angle $\beta =0$ is the ideal position of the metasurface (shown in Fig. 2). This is of paramount importance to transfer as much as possible linear momentum onto the structure. Therefore, the force exerted by the laser beam onto the metasurface is maximized, ensuring efficient propulsion. For a curved metasurface to stay stable during acceleration, its tilting from the main axis, quantified by $\beta$, should not take large values. For the stability measure, we only consider metasurfaces with a vanishing initial angular velocity. The reasons are two-fold. First, Euler’s equation of rotation can not be analytically solved for arbitrary torques that depend on the orientation of the object. Second, a numerical approach will be very time-consuming, as it has to be performed for an infinite set of initial conditions for a complete study. By performing the vanishing-initial-velocity approximation, we will arrive at a simple measure that expresses whether a surface is stable or not. In the following, we only consider metasurfaces for which the ideal equilibrium position is at $\beta =0$. We choose this as the equilibrium position because propulsion is most efficient when the surface normal of the metasurface points into the direction of the incident field.

Each curved metasurface, characterized by its curvature $R$ and normalized periodicity $\Lambda /\lambda$, is kept at different orientations characterized by Euler angles $\alpha ,\beta ,\gamma$ and is driven by the optical torque $M$. The goal is to calculate whether the metasurface rotation is towards the ideal equilibrium position or not; i.e. whether the azimuthal angle $\beta$ decreases or increases.

4.1 Stability condition

We assume a rigid metasurface with an initial orientation $\alpha ,\beta ,\gamma$, and zero initial angular velocities. Therefore, Euler’s equation of rotation in Eq. (6) in the body fixed-frame simplifies as

We will label the axes of the global coordinate system with small letters $x$, $y$, and $z$, and the body-fixed axes with capital letters $X$, $Y$, and $Z$ as seen in Fig. 2(b).

In the following, we will derive a condition that has to be satisfied for the azimuthal angle $\beta$ to decrease which is the requirement for stability. As we assume vanishing angular velocities, the rotation is only determined by the angular acceleration $\dot {\mathbf {\omega }}_B$. To determine the inward/outward orientation of the metasurface after the excitation, it is sufficient to only examine the change arising from an infinitesimal rotation around $\dot {\mathbf {\omega }}_B$. Instead of calculating the new value of $\beta$ after the rotation, it is easier to calculate the projection of the $Z$-unit vector to the $x$-$y$-plane $\rho =\sqrt {Z_x^2+Z_y^2}$.

In the body-fixed system the tensor of inertia is diagonal

Having the acceleration, we rotate the body-fixed $\mathbf {Z}$ infinitesimally around the axis $\dot {\mathbf {\omega }}_B$. For the rotation we use a first order rotation matrix for small rotation angles $\theta$ around a unit-vector $\hat {n}$, known as Rodrigues’ rotation formula, which is

With $\theta =\lvert \mathbf {\dot {\omega }}_B\rvert \epsilon$, where $\epsilon \rightarrow 0$ is used to move to an infinitesimal rotation, and $\hat {n}_B=\dot {\mathbf {\omega }}_B/\lvert \dot {\mathbf {\omega }}_B\rvert$, we can apply this rotation matrix to get the new body-fixed $Z$-unit-vector in the body-fixed system

Therefore, if $\varphi$ is smaller than zero, the torque resulting from the given angular position will rotate the surface inwards, and hence it is stable.

The stability condition measure depends on the three different Euler angles. However, we aim for a single measure for each metasurface regardless of the orientation. To solve that, we will define the cumulative stability measure.

4.2 Cumulative stability

To further compress the information, we will only consider those metasurfaces as stable, for which $\varphi$ is negative for all combinations of considered angles.

For these stable metasurfaces, we wish to furthermore define a single value to compare their stability. As the value of $\varphi$ describes how much the surface is rotated inwards, it can be treated as a gradient in $\beta$-direction of a potential in which the metasurface resides. Therefore, the stability can be measured by the depth of this potential wall along the $\beta$-axis. We will follow a weakest-link approach, where the value is given by the lowest absolute value for any possible path in the $\alpha$-$\gamma$-plane. Please note, that minimal stability translates to the maximal gradient $\varphi$. Therefore, we can define the Cumulative Stability measure as

To assign a single quantity to the metasurfaces for which $\varphi$ is positive for all angles, we followed an approach analogous to the cumulative stability. However, in this case, the minimum over $\alpha$ and $\gamma$ is used to give a lower bound for the instability. We will call those metasurfaces unstable.

Please note, while the interpretation as a potential is compelling, we use the gradient along one direction to find information about the three-dimensional potential. This can not be done fully and some information is lost. Additionally, the torque might induce an accelerated rotation around the z-axis that could not be described with a potential. Therefore, one has to be cautious about the limits of the method.

5. Results

Having all the tools at hand, we will begin our analysis by examining a metasurface made from metaatoms described only by a resonant magnetic dipole moment. It is well known that the magnetic dipole moment is the lowest order that can be driven into resonance in a high-index dielectric scatterer [37,38]. Exploiting such a lowest order resonance ensures the necessary small size and hence its low weight [39]. The operation of the particles at resonance is done to capitalize on the strong response. Hence, they are most suited for the application we have in mind.

The restriction to a magnetic dipolar scatterer is only preliminary. It brings insights to be exploited when designing realistic nanospheres made from crystalline Silicon (c-Si) that satisfy the conditions in a secondary step [40,41]. To ensure convergence, we analyze these realistic particles later up to quadrupole multipolar order within the framework of Mie theory [42–44].

5.1 Analytical Results: Magnetic dipoles

The simulations were performed with a constant periodicity of $\Lambda =100$ nm, and a $\Delta z$ that is constant for all $\Lambda /\lambda$. The parameter that we sweep is the wavelength. Thanks to the scalability of Maxwell’s equations, the results will be identical for any set of $\Lambda /\lambda$ and $\Delta z/\Lambda$ values, and all metasurfaces can be trivially scaled to operate at any desired wavelength. Therefore, sweeping the wavelength or the lattice constant are identical given that the normalized periodicity $\Lambda /\lambda$ is the same. Only when dispersive material properties are considered in the second step, the parameters are required to be fixed. What we calculate in the following is the cumulative stability $\Phi$ according to Eq. (16) for different ratios of lattice period over wavelength $\Lambda /\lambda$ and different ratios of $\Delta z/\Lambda$. This allows identifying the set of parameters that leads to metasurfaces that are rotationally stable. Moreover, from the absolute value of $\Phi$, we can conclude how stable the actual light-sail is, i.e. the smaller (more negative) this value is, the more stable the light-sail.

The parameter range for the Euler angles that were considered in the simulations are $\alpha =[0,\pi ]$, $\beta =[0, \pi /12]$ and $\gamma =[0,\pi /2]$. The reason for not performing a complete scan over $\alpha$ and $\gamma$ is that these values can be calculated by performing symmetry transformations. In the case of $\alpha$, the negative values can be calculated by performing a mirror transformation along the global $yz$-plane and keeping in mind that the sign of the $y$ and $z$ component of the torque changes. In the case of $\gamma$, the result can be periodically repeated, as the entire metasurface has a C4-symmetry. $\beta$ is only scanned for low values, as the metasurface is only examined to be stable for small perturbations.

The resulting value of $\Phi$ for metasurfaces with different $\Lambda /\lambda$ and $\Delta z/\Lambda$ is shown in Fig. 3. All stable metasurfaces are denoted in blue, and all unstable metasurfaces in red. There are several more complicated possibilities like, for example, different equilibrium positions, maximum angles to the stability, or mixed states where the behavior is difficult to judge. However, the examination of these less suited states is not of interest. Therefore, these cases are marked in grey in the figure and will not be discussed any further. By reading the data from the figure, it is immediately clear that metasurfaces with a positive $\Delta z$ are mostly stable, and metasurfaces with a negative $\Delta z$ are mostly unstable. This means that the surface should be curved away from the beam instead of being curved towards the beam like a sail catching the wind. Additionally, while some uncurved metasurfaces are shown to be weakly stable, the stability increases drastically by applying a small curvature. The sweet spot concerning the stability is around $\Lambda /\lambda =0.75$ and $\Delta z/\Lambda =0.5$ and coincides with an area with a large optical force in the z-direction. For our further investigation, we will use the most stable metasurface with $\Lambda /\lambda =0.77$ and $\Delta z/\Lambda =0.47$. Note that around the sweet spot, there is a relatively large region (dark blue), where the stability is still high. This provides design flexibility in the periodicity and curvature of the metasurface.

 

Fig. 3. The cumulative stability $\Phi$ for all stable metasurfaces in blue as well as the cumulative instability for all unstable metasurfaces in red over $\Delta z/\Lambda$ and $\Lambda /\lambda$. Additionally the corresponding value for all unstable metasurfaces is included in red. All other states are marked in grey and are discarded from further investigation.

Download Full Size | PDF

From this analysis, we know (a) what kind of particle is needed to form the metasurface (in terms of their optical response) and (b) how they need to be arranged in space to render a light-sail that is stable. This arrangement concerns both their lattice size and the desired curvature. At this point, the only leftover is the identification of an actual physical particle that offers these desired scattering properties. And of course, we would need a full-wave optical simulation to demonstrate the stability that is previously identified only in a perturbative sense. That will be done in the following.

5.2 Simulation results: Realistic particles

From now on we consider a rigid metasurface made from realistic particles. The rigidness is assumed to be enforced by having a connecting rod, made from a material low refractive index, between the nanospheres. As these connecting rods were shown to have minimal effect on the overall optical behavior of the structure [45], they are not explicitly considered.

The suggested wavelength for the Breakthrough Starshot initiative is $\lambda =1-2\mu \mathrm {m}$ as the atmospheric absorption is very low in this spectral region and the lasers can be shine on the metasurfaces without too much atmospheric power loss [10]. We will use the lower limit $\lambda =1000\,\mathrm {nm}$ as this will lead to larger optical forces that are beneficial for efficient propulsion [34].

The material used for the spherical metaatoms is c-Si. The particle is tuned to be at the magnetic dipolar resonance at $\lambda = 1000\,$nm. In the description of the scattering response from the spherical particles, we take into account up to the quadrupolar order to assure convergence. To tune the nanosphere into the resonance of the magnetic dipole moment, the radius is chosen as $r=132\,\mathrm {nm}$ assuming a permittivity for the material of $\epsilon =13.243$. Using this fixed wavelength, the lattice constant is $\Lambda =770\,\mathrm {nm}$ and the radius of curvature $R=13281\,\mathrm {nm}$ for the metasurface that was previously identified to be most stable.

To judge the stability of this metasurface, we calculated the stability $\varphi$ according to Eq. (15). This stability condition measure $\varphi$ as a function of the three Euler angles is shown in Fig. 4. In the entire parameter space that we study here, the value of $\varphi$ is negative. Therefore, it can be concluded this metasurface is stable. The cumulative stability $\Phi$ according to Eq. (16) that describes the stability amounts to $\Phi _{\mathrm {cSi}}=-0.073$ as opposed to $\Phi _{\mathrm {p}}=-0.088$ for the analytical magnetic dipoles. Therefore, the examination using just dipoles can serve as a first-order approximation for a quick survey to find a sweet spot. However, it is crucial to run a second simulation confirming the stable behavior for realistic particles. To check the stability, time evolutions of the orientation of the metasurfaces have to be performed with an initial inclination. That will be done in the following section.

 

Fig. 4. The stability condition measure $\varphi$ of the metasurface with $\Lambda /\lambda =0.77$, $R=13281\,\mathrm {nm}$ at $\lambda =1 \,\mathrm {\mu m}$ as a function of $\alpha$, $\beta$ and $\gamma$. The metasurface is made from c-Si spheres with a radius of $r=132\,\mathrm {nm}$. The initial conditions for the three trajectories considered in the following section and displayed in Fig. 4 a)-c) are denoted with green, yellow, and cyan dots, respectively in these figures.

Download Full Size | PDF

6. Simulation of the dynamics

In this section, we perform a numerical simulation of the trajectory for the c-Si-metasurface designed in the last section to observe the anticipated stability in its full dynamics.

As there is no appreciable drag in outer space, the metasurface will oscillate during the entire acceleration period. However, due to limitations concerning the length of the simulation times, it is not possible to show that the metasurface is not unstable, as the flip might just occur after the truncation of the calculated time evolution. We will, therefore, include a very small numerical drag, to prove that even the introduction of the smallest dampening will reduce the amplitude of the oscillation. This dampening can be introduced in Euler’s equation of motion as [21]

We will use the torque calculated in the last section as pre-calculated values and interpolate with splines between the sampled points if the torque needs to be known for some orientation at angles in between. Therefore, the simulation will terminate if the azimuthal angle $\beta$ is larger than $\pi /12$. The time evolution is calculated with a fourth-order Runge-Kutta method.

The time evolution, of course, sensitively depends on the initial conditions. Here, we will only consider initial conditions where the angular velocity is zero. But the metasurface has some finite potential energy as it is not placed in its equilibrium position. By assuming that no rotation is built up, the metasurface shall, in the temporal evolution, return to its equilibrium position without gaining any kinetic energy. In the presence of a small but finite coefficient of friction, however, kinetic energy is built up. The metasurface may, therefore, during the time evolution, exceed the initial azimuthal angle and therefore the maximum azimuthal angle. This happens, for example, if the initial condition of the metasurface is chosen to be an orientation where it has a rather large potential energy that can be converted into kinetic energy in the cause of its temporal evolution. Then, because of the missing data, the numerical evolution of the temporal behavior is halted. This implies that for sufficiently small coefficients of friction, we can not safely conclude that the metasurface is stable up to $\beta =\pi /12$ for all combinations of $\alpha$ and $\gamma$ as initial conditions. However, several systematic simulations have shown that the self-stabilizing behavior is retained for lower azimuthal angles independent of the choice of the other two angles.

Representative time-evolutions of the azimuthal angle $\beta$ for several different initial conditions are shown in Fig. 5. We have started all three simulations with an initial azimuthal angle of $\beta =\pi /18$ to see if the metasurface returns to its equilibrium position. The orientations of these initial conditions are indicated by colored dots in Fig. 4.

 

Fig. 5. Azimuthal angle of the metasurface as a result of the temporal evolution for some selected initial conditions. The temporal evolution has been made using a fourth-order Runge-Kutta method and $\Delta t=2.5\cdot 10^{-13}\,\mathrm {s}$. The initial conditions are marked in Fig. 4 in that order, in green, yellow and cyan. Time evolution simulations for b) and c) can be accessed in Visualization 1 (Ref. [46]) and Visualization 2 (Ref. [47])

Download Full Size | PDF

In the first case, shown in Fig. 5(a), both $\alpha$ and $\gamma$ are set to zero. In this case, the behavior is exceptionally simple. Due to symmetry reasons, only the torque along the $y$-axis is nonzero. Therefore, the trajectory is merely a harmonic oscillator with a finite dampening. Note that $\beta$ is positive definite, the deflection to the other side is realized by a jump of $\alpha$ and $\gamma$.

In the second case, shown in Fig. 5(b), some different initial conditions for $\alpha$ and $\gamma$ are chosen. This leads to a more complicated behavior. However, the overall trend is that the amplitude of the azimuthal angle $\beta$ reduces over time without major blowups.

The third case, shown in Fig. 5(c), is included to further elaborate on the discussion of the conversion of potential to kinetic energy. In this case, the initial conditions were chosen to coincide with an orientation that is characterized by a relatively high potential at a low azimuthal angle. This can be seen from the cyan dot in Fig. 4 that marks the initial conditions for this case. When looking into the emerging temporal evolution, we notice that at a later moment in time, the azimuthal orientation is larger when compared to the initial condition, i.e. there is no steady return to the equilibrium position. These peaks appear in some extended duration and repeat periodically. The peaks come up whenever the current position of $\alpha$ and $\gamma$ has a lower potential gradient (here $\alpha$ and $\gamma$ close to zero) and therefore reaches the initial potential only at a higher azimuthal angle. Nevertheless, the overall magnitude of these local peaks gets lower and ultimately the metasurface can be judged as stable.

For completeness, in the Supplemental Material, we have included two video simulations for the time evolution of orientation of the sail. The initial conditions are the ones in Fig. 5(b) (Visualization 1 in Ref. [46]) and (c) (Visualization 2 in Ref. [47]), respectively.

7. Conclusion

Identifying the composition of a light-sail that can be used to propel gram-scale spacecrafts to a fraction of the speed of light by an array of lasers is a major challenge. In particular, the light-sail has to be designed such that it always orients towards the light source. This has to be done passive and automatically, i.e. the metasurface must have a stable equilibrium position at the desired orientation.

In this paper, we have examined the stability of curved metasurfaces made from identical particles with magnetic dipolar resonance. We have shown that a small curvature applied away from the beam increases the stability of the metasurface drastically. A lattice constant around $\Lambda /\lambda =0.75$ is shown to be optimal. Full-wave optical simulations of actual metasurfaces combined with simulations on their dynamics have validated the results beyond the dipolar approximation.

To classify the behavior of many different curved metasurfaces in a systematic approach, a method based on the torque with vanishing initial angular velocities is developed. By comparing analytical results and simulation of realistic c-Si-based curved metasurfaces, we have shown, that the proposed method can be very helpful to design stable metasurfaces. Finally, we have calculated the full-time evolution dynamic of the proposed metasurface for some initial conditions and have shown their stability.

To apply this research to the Breakthrough Starshot initiative, several challenges have yet to be overcome. First, the performed time-average assumes that the particle’s motion within one wave period is negligible, which is violated for relativistic velocities. Additionally, the resulting Doppler-shift calls for a broad stability region along the $\Lambda$-axis in Fig. 3. For the metasurface to carry any payload, the size of the metasurface has to be increased drastically. This calls for a way to calculate the optical response of curved metasurfaces with a continuous description. However due to the reason that even the smallest curvature increased the stability we are confident that the results can be reproduced for the envisioned light-sails as well. Lastly, we explicitly only considered rotational stability. This is justified, as the beam diameter of the laser is many orders of magnitude larger than the spatial extent of the sail. The positional stability can therefore be considered to be less important. Of course, the position has to be considered in a complete analysis of a lightsail, which then also should account for a changing beam profile. However, to keep the degrees of freedom and the dynamics of the system as understandable and intuitive as possible, the examination of the position was neglected.

We note that the self-stabilizing lightsail considered here is entirely passive. Adding active elements to the lightsail [48], e.g. based on liquid crystal elements [49] or using phase changed material [50], might further improve the stability of the device. However, this adds equipment to the device that has to pay off eventually because it comes with the disadvantage of additional weight. A clear answer on that aspect of a possible net benefit still needs to be investigated.

We hope that our methods and designs proposed can make a small step for the big optics-enabled space journeys in the future.