Engineering phase and polarization singularity sheets

Optical phase singularities are zeros of a scalar light field. The most systematically studied class of singular fields is vortices: beams with helical wavefronts and a linear (1D) singularity along the optical axis. Beyond these common and stable 1D topologies, we show that a broader family of zero-dimensional (point) and two-dimensional (sheet) singularities can be engineered. We realize sheet singularities by maximizing the field phase gradient at the desired positions. These sheets, owning to their precise alignment requirements, would otherwise only be observed in rare scenarios with high symmetry. Furthermore, by applying an analogous procedure to the full vectorial electric field, we can engineer paraxial transverse polarization singularity sheets. As validation, we experimentally realize phase and polarization singularity sheets with heart-shaped cross-sections using metasurfaces. Singularity engineering of the dark enables new degrees of freedom for light-matter interaction and can inspire similar field topologies beyond optics, from electron beams to acoustics.


Production of 0D point singularities
One way to produce 0D point singularities is to place the singularities on the axis of a cylindrically symmetric field. Specifically, cylindrical symmetry here means that the complex electromagnetic field components (Ex, Ey, Ez, Hx, Hy, Hz) are only dependent on the radial distance r from the optical axis and the longitudinal coordinate z. Such a cylindrically symmetric field can be generated by uniform illumination of a patterned aperture that imposes a cylindrically symmetric phase profile. By minimizing the field amplitude or maximizing the zdirected phase gradient at the desired axial positions, one can produce 0D singularities at these positions.
In the demonstration system, we parametrize the cylindrically symmetric patterned aperture (located at z=-1000 μm) with 251 radial phase pixels, each spaced 4 μm apart in the radial direction, for a total aperture radius of 1000 μm. The incident vacuum wavelength is 532 nm and is linearly polarized in the x direction. We propagate the complex wavefront into the domain z>-1 mm using the vectorial diffraction integral. The objective function used for optimization is: is the axial phase gradient located at z=zi. Minimizing F maximizes the z-directed phase gradient of the optical field at the three positions z1=0 μm, z2=1 μm, z3=2 μm. To improve convergence, we use a smooth approximation to the minimum function, which has the benefit of being analytic instead of piecewise continuous: In this smooth approximation, the sum inside the logarithm will be dominated by the term corresponding to smallest value of ! . is a scale factor chosen to normalize the input array values and avoid numerical underflow or overflow during the computation of the exponential.
We perform the optimization using gradient descent, where the step size and termination condition are chosen through the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm 1 .
Supplementary Fig. 2 plots the resultant field profile from the optimization. Supplementary Fig.  2a-b plot the field intensity (logarithmically-scaled) and phase as a function of (r,z). The zdirected propagation phase has been removed from the phase plot (by multiplication with exp(ikz)) to show the long-scale phase variations. The three 0D singularities are visible as the lowintensity features along r=0 and z=0, 1, 2 μm, at the intersection of the real and imaginary zeroisolines. The radial phase mask on the patterned aperture that produces this singularity pattern is plotted in Supplementary Fig. 2c.

Optimization details for the 2D heart-shaped phase singularity
The phase-controlled surface at z=0 is parametrized with a grid of 101 × 101 superpixels, each spaced 8 µm apart in the and directions. For each superpixel, we assign a value of ϕ for the propagation delay for light incident at that position. We then illuminate the phase-controlled surface with a uniform plane wave of unit intensity and vacuum wavelength 532 nm. This produces a complex wavefront at z=0, which we propagate into the domain z>0 using the vectorial diffraction integral.
At the focal distance of = = 10 mm, the heart-shaped singularity boundary is described by a parametric curve 2 : ( ) = 16 sin & , ( ) = (13 cos − 5 cos 2 − 2 cos 3 − cos 4 ), ∈ (0,2 ) The curve positions are scaled by the scale parameter = 1 µm and centered at the origin so that the heart centroid lies at the origin in the plane. A total of 50 values of linearly spaced between 0 and 2 (excluding 0 and 2 ) are used to parametrize the curve {( ! , ! )}, = 1, … , 50. For each point ( ! , ! ) on the curve, the inward-directed normal vector Y ! is also computed: The normal vector is used to compute the directional derivative of the field propagation phase The objective function is the minimum of the squares of the directional derivatives at each point on the heart: We minimize during the optimization so that the directional phase derivatives over the singularity boundary are maximized. To improve convergence, we use a smooth approximation to the minimum function.
We perform the optimization using gradient descent with the BFGS algorithm. To reduce the risk of getting stuck in a local minimum of the objective function, we implement an iterative refinement scheme. The optimization is performed at three spatial resolutions, beginning with the lowest spatial resolution mesh (largest spatial grid size). Upon convergence, the optimized solution 0 ( 1 , 1 ) at the lower spatial resolution is linearly interpolated onto a higher spatial resolution mesh as the starting condition, and the optimization is run again using the higher spatial resolution mesh. The spatial resolution in and is improved by a factor of 1.3 at each iteration and a total of three iterations is performed, including the final iteration with the desired highest spatial resolution. The starting point for the lowest spatial resolution mesh is selected by randomly sampling 100 configurations where each phase pixel 0 is drawn from a uniform distribution between [− , ], and then selecting the configuration with the best objective function value.
In order to realize the optimized phase profile in a metasurface that operates in transmission, we discretize each of the 101x101 superpixels (pitch 8 μm) into a 32x32 grid, where the smaller grid has a subwavelength pitch of 0.25 μm. Each grid square is associated with the phase of the larger 8 μm superpixel from which it was formed. Under the locally periodic assumption and the unit cell approach, we seek to place a nanostructure (meta-atom) at each grid square to enforce the required phase at that position. This nanostructure is chosen from a "library" of meta-atoms with pre-computed optical properties, where the numerical simulation was performed under the assumption that the meta-atoms are periodic in the transverse plane. Since the heart-shaped phase singularity can be realized in a scalar field, the polarization of light is not important, and we chose to work with polarization-insensitive cylindrical meta-atoms to build the meta-atom library. The nanopillars are made of 600 nm-tall amorphous Titanium dioxide (TiO2, n=2.40 at 532 nm), spaced 0.25 μm apart in a square lattice, and mounted on a substrate of fused silica (n=1.46 at 532 nm). These meta-atoms are exhibited schematically in Supplementary Fig. 13a.
For each nanopillar diameter, we simulate a uniform array of identical nanopillars using the RETICOLO package for Rigorous Coupled Wave Analysis 3 to obtain the zeroth-order transmitted phase. The circular cross-sections of the nanopillars are discretized using a staircase approximation with 800 points over the circumference, and 81 x 81 Fourier harmonics are used in the x and y directions on the plane. The dependence of the phase and transmission efficiency as a function of the nanopillar diameter is plotted in Supplementary Fig. 13b. Nanopillars of diameter between 60 nm and 215 nm provide 2π phase coverage and are thus used to construct the metasurface based on the required phase profile from optimization.

Optimization details for the 2D heart-shaped polarization singularity
The metasurface at = 0 is parametrized with a grid of 51 x 51 superpixels, each spaced 8.4 nm apart in the and directions. We treat the metasurface as a spatially varying wave plate which manipulates the local polarization state 4,5 . This allows us to assign three values for each metaelement position on the metasurface: 0 (for the propagation delay for light polarized along the fast axis of the meta-element), 2 (for the propagation delay for light polarized orthogonal to the fast axis) and 345 (for the overall rotation angle of the meta-element fast axis). The spatially varying Jones matrix of the meta-surface can then be written as: where ( 545 ) is the 2 × 2 rotation matrix.
We then illuminate the phase-controlled surface with a uniform plane wave of unit intensity and wavelength 532 nm polarized in the ( Y + Y) diagonal direction. This produces a wavefront of transverse electric field values at = 0, which we propagate into the domain > 0 using the vectorial diffraction integral. We consider all three Cartesian components of the electric field after propagation.
One can calculate the polarization azimuth by determining the direction of maximal oscillation exhibited by the transverse electric field. At a position in the transverse plane, the direction of the transverse electric field can be written as: The polarization azimuth is calculated from the Stokes parameters of the electric field: The normalized Stokes parameters are obtained by dividing through by / , the intensity component.
The complex σ field is constructed using the normalized Stokes parameters as σ = $ + " , noting that all the Stokes parameters are real numbers. The polarization azimuth is then half the complex angle of the σ field.
One should recognize the σ field as proportional to the time average of 0 + 2 : For linear polarization, the polarization azimuth is the direction of polarization relative to the transverse -axis.
At the focal distance of f = 10 mm, the heart-shaped singularity boundary is described by a parametric curve 2 : ( ) = 16 sin & , ( ) = (13 cos − 5 cos 2 − 2 cos 3 − cos 4 ), ∈ (0,2 ) (16) The curve positions are scaled by the scale parameter s = 1 μm and centred at the origin so that the heart centroid lies at the origin in the plane. A total of 50 values of linearly spaced between 0 and 2 (excluding 0 and 2 ) are used to parametrize the curve ( ! , ! ), = 1, … , 50. For each point ( ! , ! ) on the curve, the inward-directed normal vector Y ! is also computed: The normal vector is used to compute the directional derivative of the field polarization azimuth The objective function is the minimum of the squares of the directional derivatives at each point on the heart: We minimize during the optimization so that the directional phase derivatives over the singularity boundary are maximized. To improve convergence, we use a smooth approximation to the minimum function.
To produce a metasurface that is able to achieve the required spatially variant waveplate behaviour, we use the locally periodic assumption with the unit cell approach. Each of the 51 x 51 superpixels (pitch 8.4 μm) is partitioned into a 20 x 20 grid with subwavelength pitch of 420 nm. We associate the required (φx, φy, θrot) waveplate parameters of the larger pixel with the smaller grid positions. For each grid position, we pick the meta-atom from a pre-computed library that most closely satisfies the required (φx, φy, θrot) waveplate parameters. The meta-atom library was provided by Dr Noah Rubin and was described in a previous publication 6 . The meta-atoms are rectangular nanofins made of 600 nm tall amorphous Titanium dioxide with a pitch of 420 nm. These meta-atoms are exhibited schematically in Supplementary Fig. 13c and are parametrized by the width of the nanofin in the x and y transverse directions. The dependence of the phase delay along the x-direction of the meta-atom and the transmission efficiency on each of these x and y transverse widths are plotted in Supplementary Fig. 13d and e, respectively.
Characterization of the 2D heart-shaped phase singularity metasurface The phase of the optical field at each of the 41 z-positions was obtained by a modified version of the single-beam multiple-intensity reconstruction (SBMIR) technique 7 . This process is depicted schematically in Fig. 4c. This method involves repeated application of Gerchberg-Saxton forward and backward propagation between pairs of intensity captures at different z-positions. We downsample the images by a factor of 10 in both horizontal and vertical directions to speed up the calculation process. We use a cycle starting with the image at z=9.6 mm (image 1), forward propagating to z=9.62 mm (image 2), backward propagating back to image 1, forward propagating to image 3, back to image 1 and repeating this for each consecutive image until image 41. During each propagation step, we use the true intensity profile (captured experimentally) and the retrieved phase profile as the initial boundary conditions for the propagation. The resultant forward or backward propagated field phase at the target plane then becomes the updated retrieved phase at the target plane. The starting retrieved phase estimated for all planes is set to zero. We repeat this cycle 50 times. The forward propagation is performed using the x-polarization component of a full vectorial propagator. The backward propagation is performed by replacing the wavevector k in the forward propagator with -k. The accuracy of the phase retrieval is measured by the root-mean-squared (RMS) deviation between the estimated intensity (each normalized by their respective maximum intensities) after a propagation step and the true intensity map at that plane. We observe that by the second cycle, the maximum RMS deviation over all forward and backward propagation steps within the cycle remains between 6.2% to 6.5% for every cycle from cycle 4 onwards.
Characterization of the 2D heart-shaped polarization singularity metasurface Where, A = s0 + s1/2, B = s3, C = s1/2, D = s2/2. We thus extract the intensity variation of each pixel as the quarterwave plate is rotated and fit the intensity variation-angle relationship to that exhibited above, where A, B, C, and D are the fitting parameters. The polarization azimuth can then be computed using the four-quadrant arctangent 2Ψ = atan2(s2,s1) and the polarization ellipticity angle can be obtained using 2θ = sin -1 (-s3/s0). The fitted Stokes parameters at the plane z = 10 mm are displayed in Supplementary Fig. 12e-h and are compared to the numerically simulated Stokes parameters on the same plane in Supplementary Fig. 12a-d, exhibiting excellent agreement.

Mathematical relationship between phase gradients and singularities
Here, we elucidate the relationship between phase gradients and complex scalar field zeros (singularities). Consider a complex scalar field E(r), which can be written in polar coordinates as

E(r)=A(r)exp[iϕ(r)], A(r)∈ℝ. The phase gradient is ∇ϕ(r)=Im[∇E(r)/E(r)].
If a field can be written as a finite sum of plane waves, then the magnitude of the field and its derivatives are bounded, and therefore having an infinite phase gradient ∇ϕ(r0) at a point r0 implies that there is vanishing E(r0) at the same location -a singularity. However, the reverse implication is not true -vanishing field intensities do not imply infinite phase gradients: if the field gradient ∇E(r) converges to zero as r→r0 as a linear function of E(r), the phase gradient ∇ϕ(r) does not diverge when r→r0. An example of such a reverse implication can be found in systems with purely real E(r), such as a standing wave pattern of two counterpropagating plane waves, where the phase is a constant everywhere by definition, and thus zeros in the field all exhibit a zero phase gradient as well. Apart from this special case, however, in most complex wave systems, it is more common that there exists some direction of approach r→r0 such that ∇E(r) and E(r) have a complex phase difference. E(r) thus exhibits an infinite phase gradient in that direction. In these cases, an infinite phase gradient is synonymous with an optical phase singularity.

Why a 3D volumetric optical singularity is not possible
A 3D optical singularity such as a "ball" of darkness is not physically possible. To see why, let us assume that one such singularity exists in a non-trivial field distribution (i.e., the field is not zero everywhere). By the definition of a singularity, the field and all its directional derivatives must be identically zero inside its non-zero 3D volume. An analytic field, such as an electromagnetic field, can be expressed with a series function, such as a Taylor expansion, about a point in its occupying space. If that point is within the 3D singularity that the field is zero, the field expansion across the space must be equal to zero, which contradicts the initial assumption of 3D singularity in a non-trivial field distribution. Thus, a 3D singularity cannot exist mathematically. However, it is possible to create optical fields where the field value is very low and approaches zero in a finite volume, such as in the case of a "perfect" optical vortex with a hollow (dark) center 9 . Such cases are approximate but not true mathematical 3D singularities.

Singularity shaping with the Gerchberg-Saxton algorithm
We deploy a Gerchberg-Saxton algorithm (GS) to design dark patterns and compare the results against that obtained through phase gradient maximization. The GS algorithm follows these steps: 1. Initialize the starting pixelated plane (metasurface) with zero phase and unit field intensity. 2. Initialize the target plane with zero phase and a target pixelated intensity pattern. 3. Forward propagate the starting plane fields to the target plane (located a distance z>0) away using the x-directed component of the fully vectorial Green's function 27 . 4. At the target plane, compute the pixel-by-pixel standard deviation between the previous target plane phase pattern and the forward propagated field phase. 5. Replace the target plane phase pattern with that from the forward propagated field. 6. Reverse propagate (replacing k with -k) the target plane fields from the target plane back to the starting plane using the x-directed component of the fully vectorial Green's function.
7. At the starting plane, compute the pixel-by-pixel standard deviation between the previous starting plane phase pattern and reverse propagated field phase.
8. Replace the starting plane phase pattern with that from the reverse propagated field. 9. Check if the standard deviation values from step (4) and (7) are both smaller than an acceptable tolerance. If so, terminate. If not, return to step (3) and iterate. 10. The final phase profile at the starting plane is the desired phase profile that generates the target intensity pattern at the target plane.
To give the GS algorithm the best conditions to produce dark patterns, we consider two approaches. The first approach uses a binary target pattern in which the target field intensity is either zero or unity. The results of this study are exhibited in Supplementary Fig. 8. For all these GS designs, we keep the starting plane parameters fixed (aperture size 0.8 mm x 0.8 mm, phase pixel pitch 8 μm, incident wavelength 532 nm) and also fix the propagation distance to the target plane at 10 mm. This matched the geometry used to design the 2D phase singularity sheets. As a baseline, we consider using the GS algorithm to design a bright heart outline against a dark background, which is close to the standard use cases of the algorithm. We observe that when the heart shape and outline line thickness is large, the GS algorithm faithfully replicates the desired intensity pattern. However, as the target heart intensity pattern is reduced in size towards the scale used in the phase gradient maximization method, the fidelity of this replication becomes poorer, and the heart-shaped outline cannot be made as thin as that of the target pattern. In fact, the thickness of the GS-optimized heart-shaped outline closely matches the diameter of the first Airy disk minimum (1.22l0/NA, or 11.5 μm for the geometry examined) on the target plane, for an Airy disk generated by a perfect lens of diameter equal to the square aperture diagonal, focusing light onto the target plane. The physical mechanism for this is that the target plane intensity profile is generated by a superposition of waves that have a bandlimited transverse spatial frequency. To a good approximation, any field pattern on the target plane is generated by a complex superposition of plane waves emanating from the aperture. These plane waves have different transverse spatial frequencies based on the angle of incidence onto the target plane, with a transverse wavenumber k^ that is dependent on the angle of incidence from the normal q , k^=k0 sinq. The target plane field pattern is thus a bandlimited superposition of these plane waves with different k^, with the bandlimits corresponding to the extremal plane waves with the largest angles of incidence and hence the shortest periodicities on the transverse target plane. The Airy disk diameter is proportional to the shortest periodicity (specifically, 1.22 times) and thus sets an estimate for the minimum feature size on the target plane. The smallest linewidths on the target plane as obtained through the GS algorithm appear to be limited by this minimum periodicity. In order to achieve feature sizes that are smaller than the shortest periodicity in a bandlimited function, one has to construct superoscillatory functions 10 , which require specialised modifications to the GS algorithm, such as additional custom filters 11 , or gradient-free high dimensional optimization techniques 12 .
We perform the same GS optimization for a dark heart-shaped outline against a bright background as well, and these results are displayed in the lower two rows of Supplementary Fig.  8. Unlike the case of the bright heart-shaped outline, the dark intensity patterns produced by the GS algorithm exhibit much poorer replication fidelity, breaking apart into discontinuous dark patches as the target heart-shaped pattern scales down in size. Again, the minimum feature size of bright spots in the achieved intensity pattern is similar in scale to the Airy disk diameter. Since the bright background occupies the majority of the target intensity pattern, and the GS algorithm weights each pixel equally, and the GS algorithm acts to minimize the intensity deviation over a transverse area, it prioritises replicating this uniformly bright background at the expense of the dark heart-shaped features.
A major issue facing the GS design of a phase mask to attain a binary target pattern is that the binary pattern cannot be replicated perfectly since the binary pattern is not a solution to the electromagnetic wave equation. We can improve the behaviour of a GS-optimized phase mask by specifying a target intensity pattern that is a priori known to be part of a valid wave equation solution. This is the second GS design approach that we undertook. Instead of using a binary intensity pattern, we use the actual heart-shaped phase singularity sheet intensity pattern of Fig.  3d on the target plane, which is known to be a valid wave solution since it was designed using the phase gradient maximization technique. Performing the same GS algorithm with this nonuniform intensity mask, we obtain the target plane intensity and phase patterns in Fig. 3f and 3i, respectively. Although these target plane patterns more closely replicate the desired heart-shaped dark outline as compared to the binary target patterns in Supplementary Fig. 8, they still do not achieve the fidelity of the phase gradient maximization technique. Fig. 3f and 3i thus exhibit the best performing dark patterns attained through two GS approaches, which still pale in comparison to that attainable through phase gradient maximization.
Supplementary Figure 1. Closed loop to compute the topological charge of a 2D singularity with closed transverse cross-section. Black dots and solid lines represent the loci of 1D and 2D singularity cross-sections on the transverse plane, respectively. In order to compute the topological charge of the 2D singularity with a closed cross-section without including the influence of the 1D singularity it encircles, we define a closed loop C that comprises an anticlockwise inner loop and a clockwise outer loop connected by a cut ab. Although the cut traverses the singularity cross-section, since it is traversed in both directions, its contribution to the line integral is zero. The topological charge s=∮C(∇ϕ/2π)·dr computed along C is the conserved charge associated with the 2D singularity.
Supplementary Figure 2. Three 0D point singularities in a cylindrically symmetric field with radial coordinate r and axial coordinate z. a, Intensity and b, phase profiles of the field. The 0D singularities are located at z = 0, 1, and 2 μm along the optical axis (r = 0 μm). The blue contour is the isoline at which the real part of the field vanishes and the red contour is the isoline on which the imaginary part of the field vanishes. The complex scalar field has been multiplied by exp(-ikz) to remove the rapidly varying propagation phase. c Radial phase profile of the cylindrically symmetric phase mask placed at z = -1 mm that generates the three 0D singularities upon illumination with λ0=532 nm light.
Supplementary Figure 3. Real and imaginary zero-isosurfaces of the heart-shaped singularity sheets. a Real and imaginary zero-isosurfaces of the scalar field of the heart-shaped phase singularity sheet and b, the field σ=s1+is2 of the heart-shaped polarization singularity sheet.
Supplementary The area integration for the time-averaged OAM and time-averaged energy per unit length is performed over the region between the same two heart-shaped contours as in a to yield l = -0.0011.
Supplementary Figure 8. Results of Gerchberg-Saxton (GS) phase retrieval for designing bright and dark heart outlines using different line thicknesses. The first two rows exhibit the results of using the GS algorithm for designing a bright heart shape against a dark background. The lower two rows exhibit the results of the same for designing a dark heart shape against a bright background. The geometry is the same as that used for the heart-shaped phase singularity sheet (aperture size 0.8 mm x 0.8 mm, 8 μm phase pixel pitch, 532 nm wavelength, 10 mm propagation distance). Target pattern linewidths decrease from left to right. The smallest, rightmost target heart pattern is identical in size to the heart-shaped singularity patterns studied in this paper. First row: the target pattern of a bright heart shape. Second row: achieved intensity plots (logarithmically scaled) for each target pattern, normalized to the maximum intensity achieved on that plane. The circle on the lower right corner shows the 11.5 μm diameter of the Airy disk minimum (for an aperture of diameter √2×0.8 mm) to exhibit the characteristic size of features on that plane. As the line thickness becomes much smaller than this characteristic size, the resultant intensity patterns become poorer replications of the ideal heart shape. Third row: target patterns of a dark heart shape. Fourth row: achieved intensity plots for each dark heart pattern, with the same 11.5 μm diameter Airy disk in the lower right corner for comparison. Again, pattern replication and contrast become poor as the line thickness is decreased.
Supplementary Figure 9. Cross-sections of the 2D heart-shaped phase singularity sheet at various transverse planes. The columns are, from left to right, experimental intensity profile, simulated intensity profile, experimental retrieved phase profile, and the simulated phase profile.
The sheet singularity breaks into a collection of 1D singularities away from the target plane of z = 10.0 mm. Supplementary Figure 11. Optimized waveplate parameters for the heart-shaped 2D polarization singularity sheet. a Required phase delay along the fast axis of the local waveplate φx. b Required phase delay in the transverse axis orthogonal to the fast axis of the local waveplate φy. c Required rotation angle of the fast axis of the local waveplate θrot relative to the laboratory x-direction.