Shaping caustics into propagation-invariant light

Structured light has revolutionized optical particle manipulation, nano-scaled material processing, and high-resolution imaging. In particular, propagation-invariant light fields such as Bessel, Airy, or Mathieu beams show high robustness and have a self-healing nature. To generalize such beneficial features, these light fields can be understood in terms of caustics. However, only simple caustics have found applications in material processing, optical trapping, or cell microscopy. Thus, these technologies would greatly benefit from methods to engineer arbitrary intensity shapes well beyond the standard families of caustics. We introduce a general approach to arbitrarily shape propagation-invariant beams by smart beam design based on caustics. We develop two complementary methods, and demonstrate various propagation-invariant beams experimentally, ranging from simple geometric shapes to complex image configurations such as words. Our approach generalizes caustic light from the currently known small subset to a complete set of tailored propagation-invariant caustics with intensities concentrated around any desired curve.

We assume that the beams with propagation-invariant caustics propagate in linear, homogeneous, and isotropic media, so that their rays are straight lines. Caustics are defined as the envelopes of these rays [1,2]. Each ray is described by its transverse coordinate Q = (Q x , Q y ) in the initial plane z = 0 and its transverse momentum P = (P x , P y ). The momentum is determined from the derivative of the 1D phase function Φ(φ) mapped onto the Fourier ring in cylindrical coordinates. The evolution in z of one ray is given by Q + zP [3].
For Bessel beams with topological charge l = 0 the envelope of the rays describes a circular caustic, as shown in Fig. 1(a) in the main part of this paper [3]. The rays of the astroid beam form an envelope that is characterized by four cusps, shown in Fig. 1 & 2 in the main manuscript. The left site of Supplementary Figure 1 shows the same subset of rays from several viewpoints, in order to show that the rays in the transverse plane emerge as a mapping from the three-dimensional rays to the plane. Recall that the full ray family is composed of a continuum of versions of the rays shown here displaced in the z-direction. This is demonstrated on the right site of Supplementary Figure 1, where we add for illustration purposes 5 z-shifted subsets of rays (blue), forming an invariant caustic (red). In each iteration, the additional subset is highlighted by bold lines. The 5 discrete sets should give an impression of the continuous ray picture. An animation of the rotation through the 3D subset ray volume can be downloaded here (Supplementary Movie 1), as well as the visualization of the assembly of the full ray picture and its caustic (Supplementary Movie 2). * a.zannotti@uni-muenster.de

SUPPLEMENTARY NOTE 2: ESTIMATION OF THE INVARIANT PROPAGATION LENGTH
Theoretically, beams that obey Eq. (1) of the main manuscript do not change their transverse intensity distribution, independently of any propagation distance z. However, experimental limitations like finite transverse sizes of the beams corresponding to the SLM dimensions and their imaging lenses naturally lead to finite but still comparatively long longitudinal distances (several Rayleigh lengths) in which the beams can be considered as being propagation invariant, by keeping, in good approximation, the original transverse intensity distribution I 0 (r ⊥ ) at z = 0 mm [4].
This concept is visualized schematically in Supplementary Figure 2(a) at the example of two finite plane waves with wave vectors k 1 = k ⊥ and k 2 = −k ⊥ , located on the Fourier space ring, that interfere in the far field (the real space) and form intensity fringes I(y) ∝ cos 2 (k ⊥ y). The volume in which invariant fringes occur is finite, and as shown here it has a rhombus cross-section. Its dimensions depend on the transverse size g of the finite wave contributions (due to apertures in the experimental setup) and the angle θ = sin −1 (k ⊥ /k 0 ) = sin −1 (λ 0 /a) ≈ 2 • (cf. Sec. Caustics in propagation-invariant beams), thus in turn on the light's wavelength λ 0 = 532 nm and the characteristic real space structure size a = 2π/k ⊥ = 15 µm.
In order to estimate the quality of our beams in terms of the length in which such beams can be considered to have an invariant transverse intensity profile, we investigate the propagation of a Bessel-lattice beam as an example. Note that the quality depends in parts on the individual experimental realization and the particular setup. Supplementary Figure 2(b) shows a yz crosssection through the intensity volume, obtained by measuring 397 transverse intensity distributions over a range of 40 mm, symmetrically around the origin at z = 0 mm. As one Rayleigh length of this beam is z e = 2k 0 a 2 = 5.31 mm, the length of 40 mm corresponds to 7.5 Rayleigh lengths. The five transverse intensity distributions below are obtained at the longitudinal positions indicated in the yz cross-section by red lines. They show the transition from the best experimentally realizable state at z = 0 mm, to which we refer as I 0 (r ⊥ ), to states that propagated so far that the pattern is beyond the invariant area (|z| > 11.5 mm). In order to highlight the darker areas in the yz cross-sections, for visual convenience we show the square root of the (normalized) intensity. This visualization allows to recognize the rhombus shape of the diffracting beam.
The rhombus shape is estimated and indicated with red dashed lines in Supplementary Figure 2(c). The aim of this evaluation is to get a qualitative explanation for the finite length of the area with invariant intensity. From this sketch, we estimate that the lengths of the two diagonals of the rhombus are A = (23 ± 5) mm and B = (820 ± 25) µm, hence its height can be calculated to be g = AB/ √ A 2 + B 2 = (819.5 ± 24.9) µm. The size g corresponds to the dimension of the contributing (theoretically infinitely) extended plane waves that form the Bessel-lattice beam. From our estimation, we calculate an opening angle of θ = sin −1 (2AB/(A 2 + B 2 ))/2 ≈ (2.0 ± 0.4) • , which is in agreement with the previously stated opening angle.
As a consequence of a finite wave expansion g, the Fourier ring is not infinitesimal thin, but has a size of 2π/g. Supplementary Figure 2(d) compares, true to scale, the radius of the Fourier ring 2π/a, where a = 15 µm, to the thickness of the ring. The image shows the calculated Fourier space intensity distribution. From this, we see that the ring still is comparatively thin and thus ensures several Rayleigh lengths in which the beam can be considered as propagation invariant. Further, taking into account the physical size of the SLM L SLM = 15.36 mm, imaged by the telescope L 1 -L 2 with a demagnification of M = 9.4, gives a corresponding feature size of g = L SLM /2/M ≈ 820 µm.
A good measure for the similarity of the original transverse intensity distribution I 0 (r ⊥ ) at z = 0 mm with transverse intensity distributions I(r ⊥ , z) at any longitudinal position z is the normalized cross-correlation function (denoted by ) where σ I0 and σ Iz are standard deviations of the respective intensities [5].
For the experimentally obtained 400 transverse intensity distributions, we calculate the z-dependent normalized cross-correlation γ (r ⊥ ) and show its respective maximum value in Supplementary Figure 2  intensity (sim.) origin at z = 0 mm, the cross-correlation achieves its maximal value of unity. Around the origin, a plateau is apparent, with correlation values of more than 80%, symmetrically extended over more than 20 mm. In this area, the realized Bessel-lattice beam has, in very good approximation, an invariant transverse intensity pattern. Beyond these distances, diffraction influences the beam center and the degree of the cross-correlation decreases. Similar results are existing for all presented propagation-invariant beams, as they obey the same construction rules and are realized in the same experimental system. The beams of Fig. 2 are shown as simulated intensity and phase distributions in Supplementary  Figure 3(with a very good agreement). They have comparable, experimentally obtained cross-correlation signatures. The extended plateaus around z = 0 mm prove that these beams propagate invariantly for several Rayleigh lengths. Further examples of propagationinvariant beams and the evaluation of their invariance are shown also in Supplementary Figure 6.
The analysis presented in this section shows that the length of the invariant volume can be further enhanced using a larger SLM, a smaller wavelength λ 0 , or a larger characteristic structure size a. The implementation of these enhancements is important for advanced imaging applications with a range of distances typically used in microscopy with ultra-long focal depths.

SUPPLEMENTARY NOTE 3: GENERALIZED MOMENTUM TRANSFER IN PROPAGATION INVARIANT CAUSTIC BEAMS
We discuss how a propagation-invariant beam can be made to morph into another by making it pass through a transparent mask with the appropriate azimuthallydependent phase.
Consider first the example of a propagation-invariant Bessel-lattice beam ψ BL ,q (x, y) shown in Fig. 2. As shown there, by modifying its phase with a vortex mask of charge m = 2, the field acquires an intensity pattern that changes under propagation until it settles (for z > z e ) into the propagation-invariant field (at least in the center) of a corresponding Bessel-lattice beam with changed topological charge ψ BL +m,q (x, y). This process is represented symbolically as where θ is the azimuthal angle. We generalize this behaviour to any initial and final propagation-invariant beams, ψ i and ψ f , whose phase functions are given by Φ i (θ) and Φ f (θ), respectively. The transition can be then implemented by using a transparent mask whose phase is We demonstrate this effect by transferring the three fundamental caustic shapes from Fig. 2,

SUPPLEMENTARY NOTE 4: BESSEL PENCIL FOR COMPLEX HIGH-INTENSITY FEATURES IN PROPAGATION-INVARIANT LIGHT
This section provides more experimental results and shows sophisticated high-intensity distributions realized with the Bessel pencil method. We show the spectra of the beams in Fig. 3, present several further examples whose high intensity curves describe both fundamental and rather complex shapes, and demonstrate their invariant propagation.
In order to create desired propagation-invariant highintensity curves in the transverse plane, we integrate 0 thorder Bessel beams, whose caustics are points [3], coherently along these paths. We modulate both the amplitude and phase of the spectrum, which is confined on a ring. Supplementary Figure 5 shows the calculated intensity and phase distributions of the spectra and the corresponding propagation-invariant experimentally obtained transverse intensity and phase distributions.
All images in Supplementary Figure 6 are obtained experimentally. The initial transverse intensity distribu-  tions can be seen in column 1, while column 2 shows the corresponding phases. The transverse intensity distributions after propagating one Rayleigh length z e = 5.32 mm are shown in column 3. Column 4 shows 3D intensity volumes created from taking 50 transverse intensity distributions in equidistant steps from z = 0 mm to z = z e = 5.32 mm.
We demonstrate a deltoid (row 1), and astroid (row 2), which are hypocycloids with 3 or 4 cusps, respectively. They are followed by a hypocycloid with 5 cusps (row 3). From the class of epicycloids, we realized the cardioid (row 4) and the nephroid (row 5). Further, a parabola (row 6) and a cusp (row 7) may act as building blocks for rather complex shapes, such as words. For example, we create the letters LIGHT in a propagation invariant light field, demonstrating the potential of our approach (row 8).

SUPPLEMENTARY NOTE 5: ACCURACY OF THE (PHASE) MEASUREMENTS
In this section we provide some details regarding the uncertainties in the experimental results. The dynamic range of the camera corresponds to 12 bits. Experimentally, we verified the linearity of the detected intensity (gray values) of the camera. The goodness of the linear fit is characterized by a value of R 2 = 0.99962. The discretization of the generated intensity values is limited by the 256 gray values that can be imprinted by the SLM.
Using a digital holographic method [6], we perform a phase measurement by superimposing the signal beam with a slightly tilted plane wave. The Fourier transform of the recorded fringe pattern shows the 0th, 1st and -1st diffraction order. Both higher orders contain the full information of the complex field, so that a bandpass filter is applied followed by an inverse Fourier transform of one of these higher orders. We correct all measured transverse spatial phase distributions with a reference phase, obtained by interfering a plane wave as signal beam with the tilted plane wave.
In order to characterize the accuracy of the phase measurement in our system, we measure the phase difference of a reference area on the SLM with respect to a second area where we increase the phase value gradually and compare the result to the desired phase, as shown in Supplementary Figure 7. This is done for a constant laser power as a function of the exposure time of the camera, as the error increases for low intensities, which is related to the decrease of the detectability of the interference fringes and their declining contrast. Therefore, we obtain the mean (blue values) and the corresponding standard deviation (error bars) of the measured (spatial) phase value distribution over the 256 steps. The ratio of the desired phase to the measured phase should ideally be 1 and is shown in purple. Each set of measured values is fitted with a line. The goodness of the fit is stated in each sub-image and summarized in image (b), which shows the phase uncertainty ∆φ as a function of the exposure time. We indicate the trend of the distribution of these uncertainties with the dashed red line, which follows a hyperbolic trajectory. The deviation of the measured values from a line reduces with increasing exposure time and converges asymptotically to the minimum uncertainty of our measurement system given by ∆φ ≈ 0.024 rad. For the cases presented in our work, where we use typical exposure times of 5-10 ms at the same constant laser power as for these measurements shown here, we can assume an error of the phase measurement in the order of ∆φ ≤ 0.05 rad.