Arbitrary engineering of spatial caustics with 3D-printed metasurfaces

Caustics occur in diverse physical systems, spanning the nano-scale in electron microscopy to astronomical-scale in gravitational lensing. As envelopes of rays, optical caustics result in sharp edges or extended networks. Caustics in structured light, characterized by complex-amplitude distributions, have innovated numerous applications including particle manipulation, high-resolution imaging techniques, and optical communication. However, these applications have encountered limitations due to a major challenge in engineering caustic fields with customizable propagation trajectories and in-plane intensity profiles. Here, we introduce the “compensation phase” via 3D-printed metasurfaces to shape caustic fields with curved trajectories in free space. The in-plane caustic patterns can be preserved or morphed from one structure to another during propagation. Large-scale fabrication of these metasurfaces is enabled by the fast-prototyping and cost-effective two-photon polymerization lithography. Our optical elements with the ultra-thin profile and sub-millimeter extension offer a compact solution to generating caustic structured light for beam shaping, high-resolution microscopy, and light-matter-interaction studies.


Constructing caustics into a spatial focal curve
We consider a single caustic point that propagates along a curved trajectory in 3D space (Fig. 2a-c).This trajectory produces a spatial focal curve that is fundamental for arbitrary caustic engineering.
To realize this spatial focal curve, we design the initial amplitude and phase distributions in Fourier space and then calculate the light field by using the method of angular spectrum 26 .The curve is decomposed into two analytical functions given by () Xz and () Yz describing the trajectories of the beams in the x-z and y-z planes respectively, where z is the coordinate on the optical axis.The angular spectrum integral is then asymptotically calculated using the method of stationary phase [3][4][5][6] .Thus, the solutions ( , ) xy kk form a continuous locus in the spatial frequency domain, which is expressed explicitly as a circle Yz respectively; () z  represents the radius.That is to say, the wavevectors of the light rays that form a caustic point on this focus curve constitute a circle in the Fourier space.If any wavevector ( , ) xy kk from the locus is mapped to the distance z , a two-variable function ( , ) xy z k k is obtained.From this perspective, the movement of the circle center ( ( ), ( )) kX kY zz  serves as the fundamental physical mechanism in driving the curved propagation trajectory.To achieve a proper design of spatial focal curve, the displacement of the circle center in Fourier space should be no more than the variation of radius () z  .Finally, we obtain the phase in spatial frequency domain with normalized unity amplitude: , respectively (Fig. 2a,b).The red, orange, and blue solid dots in Fig. 2a denote the positions of the centers of these circles.As the propagation distance increases from 1 z to 3 z , the center of the circle moves along the y k axis, which is correctly predicted by the previous physical analysis.The modulated optical field then passes through a converging lens (at the bottom of Fig. 2b).Thus, the conical ray bundles emitted from the circles intersect at corresponding points on the curve, which form the caustic points that propagate along the parabolic trajectory.Figure 2b illustrates the formation of this parabolic focal curve.To investigate the caustic characteristics of this curve, we simulated the transverse projections of optical ray (blue) and caustics (red) at  for calculation methods and analysis).The caustic point, which is the singularity of the gradient mapping, shifts in the positive y direction as the propagation distance increases.Moreover, the simulation results of the transverse intensity and phase distributions are provided in Supplementary Fig. 1.From the transverse profiles at different distances, the intensity maxima of the field, characterized by their caustics, exhibit an accelerated movement along the desired parabolic trajectory under wavefront control.d Schematic of one meta-atom.The parameters W, L, H and θ denote the width, length, height and in-plane rotation angle of the polymer nanofin, respectively.LCP and RCP refer to left-and right-handed circular polarized light.The height (H) and in-plane rotation (θ) enable independent control of both the amplitude and phase characteristics of the transmitted light.e, f SEM images of the metasurfaces for cases in Fig. 1c, d, together with their magnified images in the top-and oblique-views.g Schematic of the experimental setup for optical caustic characterization.LP, linear polarizer; QWP, quarter wave plate; MO, microscope objective; CMOS, complementary metal oxide semiconductor.

Design and characterization of the complex-amplitude 3D-printed metasurface
To achieve compact and lensless caustic engineering, we introduce the design and optimization of the 3D-printed metasurface with complete and independent control of amplitude and phase.The complex-amplitude information encoded in the metasurface is derived from the FT of the angular spectrum.We design the meta-atom to be a rectangular nanofin made of polymerized IP-L photoresist (refractive index ~1.52 in the visible band) (Fig. 2d).This nanofin works as a truncated waveguide with high aspect ratio up to 10.To find the amplitude and initial phase of cross-polarized light transmitted from the nanofin, we use Lumerical finite-difference time domain (FDTD) to simulate the 3D-printed meta-atom with varying height (H) and length (L), whilst the width-to-length ratio is fixed at 0.5 (see Methods and Supplementary Fig. 2a, b).Following the principles of the PB phase, we use the in-plane rotation angle ( ) of the nanofin as an additional degree of freedom to tune the phase of cross-polarized light transmitted from the nanofin.The PB phase is given by  = 2 .Thus, the complex-amplitude of the beam is fully controlled by nanofins with various geometries.
For experimental demonstrations, we choose nanofins with uniform transverse dimensions ( 400 nm W = and 800 nm L = ) and varying heights ranging from 3.0 μm to 3.7 μm (represented by the black line in Supplementary Fig. 2a, b).The characterization results are consistent with the simulation (Supplementary Fig. 2c).A meta-atom library utilizing both the height and in-plane rotation angle of meta-atoms is created (Supplementary Fig. 2d).We then fabricate two 3D-printed metasurfaces for the cases in Fig. 1c, d.Each metasurface comprises an array of 300 × 300 nanofins with periodicity 1.25 μm, yielding an area of 375 μm × 375 μm (see Methods and Supplementary Video 1 for fabrication details).The scanning electron microscopy (SEM) images of the entire metasurfaces are depicted in Fig. 2e, f, respectively.High-magnification and tilted SEM images in the insets show the details of complex-amplitude 3D-printed metasurfaces and the constituent nanofins.
To accurately construct and analyze the caustic structured light field using the 3D-printed metasurfaces, we employ a home-made optical setup as presented in Fig. 2g.First, the emitted light ( 532 nm

 =
) from the NKT supercontinuum laser source is converted to a left-circularly polarized (LCP) beam through a linear polarizer and a quarter-wave plate.The beam is then transmitted through the 3D-printed metasurface to obtain the desired complex-amplitude information of the caustic beam.Afterwards, the output is magnified by a microscope objective (Nikon, NA 0.45 = , 20×) and a tube lens ( ), then analyzed by another circular polarizer (composed of a quarter waveplate and a linear polarizer with 45° angle).Finally, these caustic beams are captured by a complementary metal oxide semiconductor (CMOS) camera.

Sculpting caustics into desired patterns propagating along arbitrary trajectories
Here, we assume the caustics of the beam at plane z is a transverse path, which is parametrically expressed as with  being the length of the path.In this case, the propagation trajectory of any caustic point on the path in the z direction is represented as . According to Eq. ( 1), the required phase in Fourier space is The key to designing the light field in real space is to construct the angular spectrum by coherent integration along the desired path  .The amplitude is chosen as to achieve uniform intensity along the path.To ensure that the interference is constructive along the desired path only, the compensation phase would be accumulated through the optical path length  by rays forming the desired caustics.The compensation phase (see Supplementary Note 3 for derivation) is given by ) where indicates the derivative of the function describing the propagation trajectory; For simplification, we can define length trajectory is governed by the length of the transverse path, and trajectory ( is jointly determined by the propagation trajectory and the transverse path of caustics at the corresponding plane z .Therefore, the spectrum distribution can be calculated by the following integral: As an illustrative example, we present a deltoid-shaped structured light beam with a parabolic propagation trajectory.The trajectory equation and the maximum propagation distance are designed the same as above.For charity, we still choose three specific planes as above, denoted as Initially, the focal points on the cross-section are arranged along a deltoid-shaped path, as depicted by the red lines in Fig. 3a.We sequentially superpose these focal points, designated as ➊, ➋, ➌, and so forth, along the transverse path indicated by the arrows in Fig. 3a.To ensure that the interference is constructive during the superposition of fields, it is essential to account for the corresponding compensation phase  at each superposed focal point, as illustrated in Fig. 3b. The compensation phase is calculated by Eq. ( 3), which comprises two subcomponents, For arbitrary plane z, this construction method is the same as described above.Consequently, based on Eq. ( 4), we can compute the amplitude and phase distribution in Fourier space, as shown in Fig. 3c.The amplitude and phase information indicated by circles ①, ②, ③ correspond to the propagation distances 1 z , 2 z , and 3 z , respectively.By using angular spectrum transformation, the desired structured light beam in physical space is obtained (see Supplementary Fig. 3 for the complex-amplitude information in real space at the initial plane).Figure 3d displays the 3D intensity distribution of the deltoid-shaped structured light beam.Figure 3e shows the caustic surface of the light beams, along with the ray pictures (blue) and caustics (red) at 0 z = , max /2 = zz , and max = zz .Simulated intensity at the detector plane can be seen in Fig. 3f.These results illustrate that the beam maintains a deltoid intensity profile while bending during propagation.Furthermore, the experimental results (Fig. 3g) agree well with the simulation, validating the correctness of our complex-amplitude 3D-printed metasurface design and the high accuracy of fabrication.To demonstrate the importance of the compensation phase, we consider four scenarios in Supplementary Fig. 4

Customizing caustics into morphed patterns with flexible trajectories
Since one caustic shape at plane z is determined by the complex-amplitude information on a ring in Fourier space, multiple caustic shapes can be created at different planes by tailoring the information on corresponding rings.Hence, the intensity profile of the beam dynamically changes along the propagation trajectory.
To demonstrate this concept, we design a morphed light field that not only follows a parabolic trajectory, but also exhibits different caustic shapes in three distinct regions: deltoid caustic when z , respectively.The roman numerals I, II, III correspond to three propagation regions.Fig. 4d presents the simulated evolving intensity volume throughout the propagation, demonstrating morphed optical caustics.Moreover, Fig. 4e shows the caustic surface, ray pictures (blue) and caustics (red) corresponding to the three propagation regions.The simulation and experimental results of the intensity profiles at different propagation distances are shown in Fig. 4f, g, respectively, with good agreement.Along the parabolic trajectory, the shape of the beam evolves from a deltoid to an astroid, and then to a hypocycloid with 5 cusps.The importance of the compensation phase in the construction of this morphed caustic field is shown in Supplementary Fig. 7.Note that the desired caustic patterns can only be shaped when the entire compensation phase  is considered.
Thus far, our approach has realized multi-dimensional customization of caustics.The engineered caustics can trace arbitrary contours, creating shapes such as deltoids, astroids, and hypocycloid with 5 cusps.These caustic patterns will propagate along anticipated trajectories.Furthermore, they can morph from one shape to another during propagation.

Discussion
To summarize, we have demonstrated arbitrary engineering of spatial caustics with morphed intensity patterns along curved trajectories.The engineered caustics are achieved by introducing a compensation phase to ensure constructive interference at every caustic point.The complex-amplitude information of the designed caustic beam in Fourier space is obtained by the method of angular spectrum and the method of stationary phase.This information is transformed into physical space and then encoded into the 3D-printed metasurfaces, which are fabricated by cost-effective, rapid-prototyping, and high-resolution TPL.By varying the height and in-plane rotation angle of each meta-atom (a polymer-based nanofin), we achieve independent and complete control over both the amplitude and phase of the transmitted cross-polarized light.Experimental lensless constructions of caustic structured fields with arbitrary spatial distributions show good agreement with the simulation results.Our approach extends the limited scope of caustic light to arbitrarily customized spatial caustics, with intensities precisely focused along predetermined curves.We anticipate that our results will inspire new ideas for designing more complex propagation scenarios of incoherent caustic fields in future works.Finally, arbitrary caustic engineering can be used in many potential applications, such as optical high-resolution imaging, optical micromanipulation techniques, and nanofabrication and processing.

Optical simulation of meta-atoms
The amplitude and initial phase retardation of cross-polarized light beams transmitted through a meta-atom are simulated using the commercial software Lumerical FDTD.The boundary conditions are set to periodic in the x and y directions, and set to perfectly matched layers in the z direction.In the simulation of nanofin, a constant pitch distance of 1.25μm and a fixed width-to-length ratio of 0.5 are used.The nanofin's height and length are set as sweeping parameters.According to the analytical results based on Jones matrix calculus, the amplitude and initial phase response of a meta-atom are calculated as axis respectively (a detailed derivation is provided in Supplementary Note 4).To design the metasurface with a precise phase distribution, the initial phase introduced by nanofins with different heights is compensated by the PB phase.The initial phase retardation is calculated with respect to a plane above the highest nanofin.

Fabrication of complex-amplitude 3D-printed metasurfaces
The metasurfaces are fabricated by a Nanoscribe Photonic Professional GT machine, based on TPL within a tightly focused femtosecond laser beam.The polymer metasurface samples are fabricated on indium tin oxide-coated glass substrates using a ×63 Plan-Apochromat objective ( NA 1.40

=
) in the dip-in configuration and IP-L 780 resin (refractive index~1.52 in the visible band).In the fabrication process, ContinuousMode in GalvoScanMode is used, with a scan speed of 7000μm/s, LaserPower of 90 mW, Slicing (laser movement step in the longitudinal direction) distances of 20 nm, GalvoSettlingTime of 2 ms, and PiezoSettlingTime of 20 ms.Each scanning field is confined to a square unit cell with a size of 100 µm×100 µm.After the exposure process, the samples are sequentially immersed in propylene glycol monomethyl ether acetate for 15 min, isopropyl alcohol for 5 min (simultaneously exposed by a Dymax BlueWave MX-150 UV LED curing system set at 70% maximum power), and methoxynonafluorobutane for 7 min.The samples are then allowed to dry in ambient air through evaporation.

Fig. 1
Fig. 1 Principle of arbitrary engineering of spatial caustics based on 3D-printed metasurfaces.a Schematic of customized optical caustics, generated by the metasurface.b(i)-d(i) Amplitude and phase information of the desired caustic beams in Fourier space, together with a FT lens.b(ii)-d(ii) Propagation dynamics of realized caustic structured light in real space.
limits the maximum propagation distance max z of the caustic point.More details of the derivation are provided in Supplementary Note 1.We present an example of the spatial focal curve with a parabolic trajectory denoted as field.The calculated amplitude and phase distributions in Fourier space are shown in Fig.2a.We use red (①), orange (②), and blue (③) dashed circles to represent the complex-amplitude distributions corresponding to the propagation distances

Fig. 2
Fig. 2 The design principle of arbitrary spatial caustic engineering, the fabrication of complex-amplitude 3D-printed metasurfaces, and the characterization of optical caustics.a Initial complex-amplitude distributions of the spatial focal curve in Fourier space.The dots are the circle centers corresponding to different propagation distances.b Rays emitted from the circles intersect on a specified focal curve.c Visualization of the caustic curve and transverse projections of rays (blue) aligning with the caustics (red) at 0 z = , max /2 zz = , and of the path with l being the total length of the path.

3 z
, showcase the process of constructing the previous caustic point into a predefined pattern in the transverse plane.
Fig. 3b1.The details of the phases are presented in Fig. 3b2.Notably, length  exhibits high-frequency oscillations, while trajectory  resembles a square wave distribution.
: without compensation phase  , with length with  .We analyze the complex-amplitude distributions in Fourier space and intensity distributions at plane max /2 zz = .Evidently, the desired complete caustic pattern is sculptured only when length  and trajectory  are combined.Using this design method, we customize a Z-shaped caustic that follows a parabolic propagation trajectory (Supplementary Fig.5).Further results of the measured intensity distributions of geometric-shaped and letter-shaped caustics at different wavelengths are shown in Supplementary Fig.6a, b, respectively.

Fig. 3
Fig. 3 Design and characterization of deltoid-shaped caustic beams with a parabolic trajectory.a Schematic of the optical caustic paths at transverse planes.b Compensation phase with components

. 1 z , 2 z , and 3 z
We use the same trajectory equation and maximum propagation distance as described in the previous sections.The transverse paths of caustic points (red lines) at three propagation distances are shown in Fig.4a.Similar to the previous approach, we continuously superpose caustic points such as ➊, ➋, ➌, and so forth, following the direction of the arrows.At each point, we evaluate the corresponding compensation phase  (Fig.4b) to ensure constructive interference.Schematics of the compensation phases ( functions are shown in Fig.4b1,The computed complex-amplitude distributions of the angular spectrum are illustrated in Fig.4c, where the circles labeled ①, ②, ③ correspond to the propagation distances 1

Fig. 4
Fig. 4 Design and characterization of caustic beams with shapes varying from deltoid to astroid, and then to hypocycloid with 5 cusps, along a parabolic trajectory.a Schematic of the optical caustic paths at different transverse planes.b Compensation phase, composed of length  t are the S-parameters for the polarization along the long axis and short