Vortex laser arrays with topological charge control and self-healing of defects

Geometric arrays of vortices found in various systems owe their regular structure to mutual interactions within a confined system. In optics, such vortex crystals may form spontaneously within a resonator. Their crystallization is relevant in many areas of physics, although their usefulness is limited by the lack of control over their topology. On the other hand, programmable devices like spatial light modulators allow the design of nearly arbitrary vortex distributions but without any intrinsic evolution. By combining non-Hermitian optics with on-demand topological transformations enabled by metasurfaces, we report a solid-state laser that generates 10 × 10 vortex laser arrays with actively tunable topologies and non-local coupling dictated by the array’s topology. The vortex arrays exhibit sharp Bragg diffraction peaks, witnessing their coherence and topological charge purity, which we spatially resolve over the whole lattice by introducing a parallelized analysis technique. By structuring light at the source, we enable complex transformations that allow to arbitrarily partition orbital angular momentum within the cavity and to heal topological charge defects, thus realizing robust and versatile resonators for applications in topological optics. A solid-state laser that generates 10 × 10 vortex laser arrays is demonstrated. The topologies are actively tunable.

I n many natural systems, ensembles of vortices can self-organize into regular patterns through their mutual interaction 1 . This is the case for clusters of cyclones encircling the poles of Jupiter 2 , knotted lattices in chiral liquid crystals 3 and motile haploid cells of sea urchins swimming in vortex trajectories 4 . In broad-area lasers, the crystal-like arrangements of optical vortices can also form spontaneously due to the nonlinear interaction of several competing transverse modes 5,6 . These crystals consist of arrays of intensity nulls, each corresponding to a phase singularity of the field carrying orbital angular momentum (OAM) with a topological charge 7,8 . The self-organizing nature of these vortex crystal lasers solves a difficult energy minimization problem, making them relevant in many areas of physics. Although fascinating as a nonlinear dynamic system, these optical vortex crystals are unsuitable for practical applications in topological optics, as the laser states can be non-stationary 9 and their topology cannot be tuned at will. Thus, the applications requiring vortex arrays 10,11 , including parallelized super-resolution imaging 12 , OAM-multiplexed communications 13 and multiparticle micromanipulation with optical tweezers 14,15 , have turned their attention to alternative generation schemes, ranging from spatial light modulators (SLMs) and Dammam gratings 16 to metamaterials [17][18][19][20][21] . However, all these extracavity implementations carve out the vortices from an input beam profile, resulting in reduced efficiency. Moreover, these devices only realize the photonic design programmed by the user without exhibiting an internal evolution, thus losing the appeal of an active system. A recent approach is to coherently combine multiple microring lasers based on supersymmetry, which is promising for power scaling but is currently limited to a topological charge of one 22 .
In this work, we present a structured light laser 23 capable of producing vortex laser arrays with a topology that is directly tailored at the source. Just as in a broad-area laser, our gain medium has a large transverse cross-section supporting many transverse modes. Instead of allowing these modes to structurally self-organize in the resonator, we force them to form a lattice of a hundred laser beams by inserting a mask that modulates the phase and amplitude of the field 24,25 . The phase profile is designed to impart arbitrary OAM values to each beam of the lattice; consequently, we artificially introduce an amount of OAM in the system that must be conserved. This effectively partitions the resonator into two sections-one on either side of the OAM-transforming mask-each emitting a laser array imbued with a different topological charge. We spatially couple the vortices using a non-Hermitian mechanism that introduces dissipative losses in the resonator, making the system stable and coherent, and resulting in a 'crystallization' of the array phase. The coupling network is not limited to the nearest neighbours as for arrays of Gaussian lasers with no topological charge, but instead can be tailored to mix vortices that are widely separated in the lattice. Although different schemes for generating a single vortex inside a resonator have been demonstrated 19,[26][27][28][29][30] , this work accomplishes the conditions for generating vortex laser arrays in a single cavity, realizing a platform to explore complex topological transformations and collective vortex effects at the source.

Results
Operating principle of the metasurface laser. The resonator comprises a diode-pumped neodymium-doped yttrium aluminium garnet (Nd:YAG) laser with a high Fresnel number in a self-imaging (4f) configuration (Fig. 1a), supporting a very large number of transverse modes (~10 5 ) (ref. 31 ). Two intracavity lenses, which form a telescope system, image the plane of the metasurface array onto the opposite end of the cavity, where an output coupler (OC) is present (Materials). The other OC is displaced from the metasurface array to couple the laser beams via the Talbot effect. The metasurface array provides the necessary intracavity topological charge conversion. It consists of 10 × 10 metasurfaces in a square-lattice geometry with amorphous Si nanopillars on a fused SiO 2 substrate (Fig. 1b). Different vortex-plate designs are possible, either based on the Pancharatnam-Berry phase (such as q plates 32,33 that are used in this work) or on the propagation phase such as (J plates 34 ; Supplementary Information). The general requirement is that the circulating beam does not acquire any net topological charge on a double pass through the metasurface array. This condition cannot be satisfied using a polarization-independent spiral phase plate in a plane-parallel resonator 35 , but we achieve this by exploiting the engineered birefringence of a metasurface, in combination with the appropriate intracavity polarization optics, for complete spin-orbit control of the beam (Supplementary Information provides details on how the round-trip condition is satisfied). The spin-orbit conversion is provided by the geometric phase elements of the metasurfaces, which convert circularly polarized light carrying spin angular momentum (SAM) into a beam with opposite handedness and an OAM quantity dictated by the spatial orientation of the nanopillars 34 .
The telescope configuration gives direct experimental access to the far-field (FF) plane, which is physically located between the two intracavity lenses and contains the optical Fourier transform (FT) of the near field (NF) of the metasurface array. By placing a pinhole in the FF, it is possible to filter high-spatial frequency components of the beam without reducing the gain effective area or affecting the laser efficiency 31,36,37 . Since the metasurface array separates the cavity into two sections with different topological charges, an FF pinhole allows us to force the existence of the Gaussian beams-having the smallest spatial frequency content in the FF-on the telescope segment of the cavity, leaving the vortex beams in the Talbot section.
The interstices between the metasurfaces of the array are filled with a metallic mask that spatially filters the beams' transmission. The beams diffract at every metasurface aperture, resulting in an array of interfering vortices (Supplementary Videos 1 and 2). Due to the periodic arrangement of the vortices, the Talbot effect ensures that self-images of the array are created at specific planes along the propagation direction (Fig. 1d, Talbot carpet), whose positions, in general, depend on the phase relation among the vortices. At the characteristic Talbot distance z T = 2d 2 /λ of a square array, where d = 0.3 mm is the period and λ = 1,064 nm is the laser wavelength, an array of both in-phase and out-of-phase vorticeswhere the nearest neighbours have a phase difference of 0 and π, respectively-will form a self-image 37 . This degeneracy of the phase solutions is lifted at z T /2, where only an out-of-phase array can produce a self-image; therefore, we use this distance in our system. Since the Talbot propagation path is folded in the resonator, this requires us to set a distance of z T /4 = 42 mm between the metasurface array and the nearest OC (Fig. 1a). The amplitude mask combined with Talbot propagation allows to select the lowest-loss phase solution for the array over several round trips. Since the lasers of the array operate at several frequencies corresponding to different longitudinal modes of the cavity, a phase-locking mechanism is needed to synchronize them 24 . In our cavity, the linear coupling together with the saturable nonlinearity of the gain medium provide mode competition and favour the minimal-loss solution ( Supplementary Fig. 13).
The generation of a vortex laser array is demonstrated in Fig. 1c. Different NF intensity distributions are observed at each end of the cavity: one corresponding to an array of Gaussian beams and the other to an array of donut beams with dark cores, signalling the occurrence of a topological charge conversion inside the cavity. The phase coherence of the arrays is witnessed by the Bragg diffraction peaks present in the FF distribution, obtained by the FT of the NF distributions. The Bragg pattern of the vortex array can be easily understood using the convolution theorem. Since a vortex array corresponds to the convolution of a Dirac comb with a vortex beam (Fig. 1e, first column), its FT is simply given by the point-wise product of the FTs of the Dirac comb and the vortex, resulting in a single, large pixelated vortex (Fig. 1e, second column). The experimental results are in very good agreement with the numerical Fox-Li simulations of the laser, which include the saturable gain nonlinearity, proving that the measured lasing mode is the one experiencing the lowest loss in the cavity (Supplementary Fig. 6 and Methods).
Parallelized topological charge analysis. The Bragg diffraction peaks are a signature of phase locking in the vortex laser array, but they do not directly show the value of the topological charge carried by the vortices. A standard approach to analyse the topological charge of a vortex is to perform a modal decomposition by measuring the inner product between the field and appropriate correlation filters. In general, a single beam could be selected by applying a digital aperture on an SLM in the NF of the array and its charge revealed by displaying a phase and amplitude hologram of the opposite charge, leading to an intensity peak at the beam's optical axis in the FF 38 . However, this is an inefficient technique for an array, with a long sequential processing time proportional to the number of vortices. By exploiting the fact that each vortex in the array is imaged to the same FF profile regardless of its position in the array (Fig. 2a) and by using the linearity of the FT, we introduce a parallelized technique to analyse the topological charge of each beam in the array all at once. Using an SLM, we carry out a modal decomposition of all the vortices in the FF and measure their respective on-axis intensity in the NF (Materials). The modulated beams will now have a charge given by the sum of the decomposition charge displayed on the SLM and the intrinsic one of the array. This corresponds to a lattice of Gaussians when the SLM and array charges have equal values but opposite signs, and to an array of donuts of variable radii in all the other cases (Fig. 2f). In Fig. 2b-e, we present the results for a metasurface array with charge ℓ a = 1. The measured charge purity statistically averaged over the array is near 100% for right-and left-circularly polarized light. We also note that the total summed OAM of the beam is conserved in the NF-to-FF propagation: in the NF, we have 100 beamlets separated in space, each with a charge of ℓ a , whereas in the FF, we have 100 superimposed beams-still with a charge of ℓ a -that interfere and produce a pixelated donut beam. By changing the intracavity polarization optics, it is also possible to generate vector vortex arrays with azimuthally and radially polarized light, corresponding to a superposition of opposite charges (ℓ = ±1; Fig. 2d,e). In this case, a small contribution at ℓ = 0 is also present (4%) due to a residue of the beam that is not converted by the metasurface array and not filtered outside the cavity (Supplementary Information). In the rest of this work, we will concentrate on scalar vortex beams, keeping the focus on the topological optics aspects of the laser.
Scaling up the topological charge. Our metasurface laser scheme allows to scale up the charge of the vortex array by increasing the azimuthal gradient in the phase design of the metasurfaces. We characterized an array with a charge of 5 ( Supplementary Fig. 7). Due to a charge-dependent increase in the vortex core, the typical diameter of the donuts observed in the NF is 150 ± 10 μm (Fig. 3a), larger than that of the charge-1 donuts shown above with 91 ± 5 μm (Fig. 1c)-a size change in agreement with the numerical simulations (138 ± 8 and 90 ± 9 μm, respectively). The larger topological charge content of each laser in the vortex array is also manifested Fig. 2 | Parallelized topological charge characterization of a vortex laser array. a, Each vortex in the array is projected to the same location in the FF plane using a lens. In the FF, a hologram encoded to modulate both phase and amplitude is displayed on an SLM to carry out the charge analysis. A second lens re-images the entire array in the NF, where the new beam profiles carry information about the charge of the input vortices. The SLM operates in reflection with a first-order diffraction grating, but here it is shown in transmission and at zero order for simplicity. b-e, Topological charge spectra statistically averaged over the whole array for a metasurface array with a charge of one integrated in the laser. The different cases correspond to left-circular (b), right-circular (c), radial (d) and azimuthal (e) polarization outputs. The dots indicate the mean and the error bars indicate the standard deviation of the values in the vortex array. f, Topological charge decomposition of each vortex in the array for the laser with left-circular polarization (top). Images of the modulated array obtained with different SLM holograms are also shown (bottom).
in the FF distribution (Fig. 3a), corresponding again to a pixelated vortex but with a dark core larger than the charge-1 case.
The topological charge purity remains high, with a peak at ℓ = 5 around 90% (Fig. 3b). The small spurious contributions at the other charges are due to imperfect Talbot reconstruction in the peripheral regions of the array. This can be understood by examining the coupling network of the array, which is more complicated than for arrays without topological charge, where coupling is dominated by the nearest neighbours (ℓ = 0; Fig. 3c). The coupling matrix κ i,j of the vortex array is non-Hermitian 39,40 , being complex, with an imaginary component due to transmission losses at the mask, and symmetric, like the field propagation from one vortex to another; thus, κ * i,j ̸ = κ j,i = κ i,j . We calculate the magnitude of κ i,j for any position j in the array and i corresponding to a central position 37 . This quantity is normalized to the self-coupling, which is the amount of energy that laser i transfers back to itself. Considering a vortex array with ℓ = 5 (Fig. 3c), the coupling is strong, as the main energy comes from relatively distant vortices lying on a large circle. On the other hand, the self-coupling tends to zero due to the intensity null of the vortex that-after expansion in a Talbot round trip-transfers negligible energy to the initial site. When i corresponds to the position of a vortex far from the centre of the array, the power coupling will be anisotropic due to the finite size of the array, forming a vortex with an asymmetric donut intensity profile and fractional OAM producing small sidebands in the topological charge spectrum. These effects are less pronounced for small charges (Fig. 3c; ℓ = 1) and can be mitigated by increasing the array size.
Partitioning OAM inside the cavity. The spin-orbit transformations implemented in our cavity, in principle, allow for several topological solutions, each comprising a pair of vortex arrays with different topological charges occupying the two sides of the metasurface array. We denote a generic topological solution as (m|n), where m and n are the charges of the arrays on the telescope and Talbot segments of the cavity, respectively. OAM conservation imposes the condition m + ℓ a = n, where ℓ a is the metasurface array charge. Until here, the use of the FF pinhole forced the appearance of beams with zero charge in the telescope section, imposing the solution (0|ℓ a ).
After removing the FF pinhole, there could be many degenerate solutions. However, this degeneracy is lifted by two main elements in the cavity, giving topological charge-dependent losses and favouring a specific solution. Both types of loss are dissipative and originate from the spatial filtering of the mask, after the beam has made a round trip in each segment of the cavity (Fig. 4a). The first contribution relates to a round trip in the telescope section. We find that a small deviation Δz in cavity length from the ideal 4f value produces a topological charge-dependent transmission efficiency through the mask, due to the divergence of the beams (Fig. 4c). The second contribution is connected with a round trip in the Talbot section. Due to the finite size of the array, part of the Talbot self-image does not overlap with the original array but rather spills out around its edges 41 , getting blocked by the mask. Since the spillover increases with the divergence of the diffracting beams, the transmission efficiency diminishes with the charge (Fig. 4d). The two described effects are in competition and the lowest-loss solution chosen by the laser depends on Δz.
By combining the telescope and Talbot efficiencies, we calculate the total transmission efficiency for a full round trip in the resonator. As shown in Fig. 4e, in the case of ℓ a = 1, there are two possible topological solutions, namely, (0|1) and (−1|0), dominating in different ranges of the tuning parameter with a transition point around Δz/4f = 0.8%. We experimentally verify the existence of these two solutions by displacing the OC on the telescope side by a few millimetres and measuring the change in topological charge spectrum of the array emitted from the Talbot side of the laser (Fig. 4a). Remarkably, the beams are donut shaped even when they do not carry any charge. This feature arises from the phase-only transformations operated by the metasurfaces, which unwind the phase but leave an amplitude term in the NF of the beams, as also observed in the laser simulations (Supplementary Information). These findings can be generalized for higher array charges, resulting in a set of ℓ a + 1 different topological solutions, ranging from (0|ℓ a ) to (−ℓ a |0) (Fig. 4b).
Healing of defects. So far, we have presented vortex laser arrays with topological charges determined by arrays of identical metasurfaces. Next, we break the symmetry of the array by introducing an intentional topological charge defect and study the properties of the resulting vortex array. We fabricate a 10 × 10 array with metasurfaces carrying a charge ℓ a = 1, except for a defect device with charge ℓ d = 2 located in the central region of the array (Fig. 5a). We remove the FF pinhole from the cavity and tune its length to sustain the (0|ℓ a ) solution. At the output of the metasurface laser, we observe a vortex array with donut beams of the same size (Fig. 5b), as in the defect-free laser (Fig. 1c). Further examination of the spatially resolved topological charge spectrum reveals that in correspondence to the defect metasurface, the dominant component is ℓ = 1, as for the rest of the array, indicating the ability of the metasurface laser to heal the topological charge defect (Fig. 5c). The defect recovery engine built into our laser is based on the topological version of the Talbot self-healing effect 42 : multiple passes through the Talbot section of the cavity remove non-periodic topological faults in the array, only strengthening the periodic components at the self-imaging plane. Following a full round trip through the cavity, we can see how the topological charge is repeatedly converted at the defect metasurface, resulting in a self-consistent solution ( Fig. 5d; for the defect healing transient, see Supplementary Video 3). The additional charge introduced by the defect device does not disappear from the system, but, instead, it manifests itself in the telescope section as a beam with a new charge ±(ℓ a − ℓ d ). As long as this beam can propagate through the telescope segment, any value of the topological charge of the defect will be recovered in the vortex array with high efficiency (Fig. 5e). Since this beam carries a charge larger than the surrounding array of the Gaussian beams (that is, in our experiment, |ℓ a − ℓ d | = 1 > 0), it suffers more losses by self-imaging in the telescope (Fig. 4c at Δz/4f = 1%), resulting in a donut with lower intensity than its neighbours (Fig. 5b). The fact that in our non-Hermitian system 39,43 , the presence of the defect increases the losses in the cavity, resulting in a change in the gain condition, is in sharp contrast with topological insulator lasers 44 , which are also robust to defects but feature a gain that is not affected by their presence owing to topological protection-a concept that does not easily lend itself to solid-state lasers. Finally, we note that multiple defects can be healed in the laser, provided that each defect is surrounded by a suitable neighbourhood as determined by the vortex coupling network (Fig. 3c). Supporting numerical studies investigating the dependence of self-healing on different properties of the system, such as coupling, defect location, charge and density, are described in the Supplementary Information.

Discussion
We have demonstrated a metasurface laser capable of controlling topological transformations in optical vortex arrays based on arbitrary spin-orbit conversions. The underlying non-Hermitian mechanism provides stability and coherence to the system, as well as the ability to partition OAM within the cavity. The concepts presented in this work open up several new possibilities. The spatial coupling geometry of a laser array could be engineered by judiciously tailoring the topological charge of the beams. Anisotropic configurations could also be realized by employing arrays with fractional OAM. This could serve to individually control mutually coupled channels in a network of communicating nodes. We have focused here on the steady-state generation of optical vortex laser arrays, but in the future, their charge can be dynamically modulated by adjusting the cavity length, either mechanically or optically. By combining the two laser arrays emitted from the opposite sides of the cavity, one can utilize the interference of the vortex array with the complementary zero-topological-charge Gaussian array to generate topologically protected lattices known as skyrmions 45 , endowed with the tuneable feature of our system. By exploiting the design freedom allowed by metasurfaces, many other topological coupling schemes can be explored, such as non-reciprocal OAM conversions, where the imparted charge depends on which side of the  metasurface is seen by the beam propagating through the cavity. The healing of defects more generally shows that the system can respond and adapt to aperiodic features of the phase mask. This opens up the possibility of iteratively manipulating information encoded as topological charges and could be used to develop photonic simulators 46 using OAM as a synthetic dimension, including an optical simulator of synthetic gauge fields 47 , and implementing quantum search algorithms in a laser solver with a marked array 48 . For instance, our laser arrays could be used to model the population dynamics of other physical systems, such as interacting galaxies or hydrodynamic vortices, mapping their angular momenta to different topological charges of the laser. In particular, it will be interesting to explore what happens when the system is tuned away from the regime of strong coupling and whether distinct domains of different topological charges can form in the laser array. Finally, our scheme is scalable and allows for the generation of large optical vortex arrays that could exploit spatial-division multiplexing 49 for optical communications using multidimensional structured light 50 with OAM, frequency and polarization degrees of freedom.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41566-022-00986-0.

Methods
Metasurface arrays. The metasurfaces consist of 600-nm-tall Si nanopillars arranged on a hexagonal close-packed lattice and lying on a 500-μm-thick fused SiO 2 substrate. The phase library of nanopillars was simulated for the Nd:YAG laser wavelength of 1,064 nm using a finite-difference time-domain software (Lumerical) and a complex-value refractive index measured by ellipsometry. Each metasurface has a diameter of 200 μm and is arranged in a 10 × 10 square array with 300 μm period. For these array parameters, the Talbot distance is z T = 169 mm. The metasurfaces employ either geometric or propagation phase to impart a topological charge to the beam. The interstices between the metasurfaces are filled with a thin Au layer, which serves as an amplitude mask. The fabrication of the metasurface arrays includes the following main steps. A 600-nm-thick Si film was deposited onto the fused SiO 2 substrates by plasma-enhanced chemical vapour deposition (STS LPX PECVD system). Before electron-beam exposure, a 10-nm-thick Au film was deposited on top of the poly(methyl methacrylate) layer by physical vapour deposition, thus ensuring that there were no charging effects. Electron-beam direct writing of the metasurface arrays was carried out using an ultrahigh-resolution Raith150 Two electron-beam lithography system. During exposure, a beam energy of 20 keV was employed, setting the beam current to 140 pA. The Cr metasurface pattern obtained through lift-off was transferred to the Si film by means of a plasma-enhanced reactive ion etching system (Sentech SI500). The etched depth and sidewalls verticality were optimized by carefully tuning the process parameters (radio-frequency power, 18 W; pressure, 0.32 Pa; inductively coupled plasma power, 200 W; C 4 F 8 , 32 s.c.c.m.; SF 6 , 30 s.c.c.m.; Ar, 10 s.c.c.m.; temperature, 0 °C). The amplitude mask was finally realized using standard ultraviolet lithography (Süss MicroTec MA6/BA6 mask aligner system) and chemical etching.
Metasurface laser. The active medium of the laser is an Nd:YAG module with 40 laser diode bars (Cutting Edge Optronics, Northrop Grumman). The Nd:YAG rod is low doped (<1%), 146 mm long and 7 mm in diameter. The module is water cooled at 26 °C, and operated in the quasi-continuous-wave mode with 250 μs current pulses at 1 Hz repetition rate to avoid thermal lensing effects 25 . Different intracavity optics schemes can be employed depending on the type of phase used in the metasurface design. We first align the resonator using a simple array of apertures (not generating a vortex array), and then we replace this by the array of metasurfaces. The intracavity lenses forming the telescope have an anti-reflection coating designed for the lasing wavelength, and a focal length of 150 mm. The OCs have a reflectance above 97% (EKSMA Optics). The FF pinhole spatially filters the Fourier spectrum of the Gaussians, making their intensity distribution more uniform and improving the quality of the vortices generated from the metasurface. At the same time, the pinhole also blocks possible reflections from the metallic mask of the metasurface array. The FF pinhole can be removed if a non-reflective (non-metallic) amplitude mask is employed or when using J-plate metasurfaces in the vector-vortex-beam cavity configuration.
Numerical simulations. We carry out numerical simulations of the metasurface laser cavity based on the iterative Fox-Li method 51 , including the nonlinearity of the active medium via a gain matrix 37 . We start every simulation from a numerical noise seed to generate the initial complex field, which is then run through the cavity using the Fourier beam propagation method with an absorbing boundary layer 52 . In particular, the propagation through the telescope section is not realized by means of FTs converting the beam from the NF to the FF and vice versa; we rather implement phase matrices for the lenses and propagate the beam over the actual distances corresponding to our experimental setup. This allows calculating FF patterns in scale with the experimental ones, as well as to account for the finite aperture of the intracavity optical elements. Topological charge transformations are accounted for by implementing two different phase profiles on the opposite sides of the metasurface array, as those seen by the propagating beam.
Topological charge analysis. The setup for parallelized topological charge analysis consists of a 4f telescope with a phase-only SLM (HOLOEYE GAEA-2) lying at the Fourier plane in the middle of the two lenses. The SLM operates at the first diffraction order in the reflection mode, although in Fig. 2a, it is shown in the transmission mode for simplicity. As an initial charge analysis, we decompose the beam using a pitchfork hologram, corresponding to the interference of a spiral phase with a diffraction grating (for the operation at the first order). However, to correctly extract the relative weights among the different components of the topological spectrum, complex amplitude modulation is needed 38 . This translates into an additional phase mask showing up as an annular element in the SLM holograms. To define the hologram used in the complex amplitude modulation, we first fit the radius of the vortices imaged in the NF of the metasurface array and extract the corresponding beam waist ω 0 of the embedded Gaussian 53 . Then, we simulate a single donut with the measured ω 0 and the charge deduced from the pitchfork analysis, propagate it to the FF and extract the waist of the corresponding FF beam, which is used to program the SLM hologram. As the SLM only modulates horizontally polarized light, the circularly polarized light of the vortex laser array generated from q-plate arrays is converted into horizontally polarized light using a combination of quarter-and half-waveplates (not shown in Fig. 2a). A pinhole is used to block unwanted diffraction orders before the camera (Spiricon LT665) acquisition. The modal decomposition of the array is carried out over a basis of Laguerre-Gaussian modes with topological charge ranging from −2ℓ to 2ℓ, where ℓ is the array charge and p = 0 radial mode. The optical axis of the setup was aligned by sending an array of Gaussian beams generated using a calibration mask of apertures instead of the metasurface array. The modal decomposition and normalization of the topological charge spectra was carried out according to the method described elsewhere 38 . For each individual beam, the topological charge spectrum was extracted by measuring the intensity at their optical axis for every detection mode and then normalized by their sum. We employed background subtraction to remove noise from the acquired intensity images and set a noise threshold to exclude contributions of vortices that are too weak (<2%) from the statistical average of the array.

Data availability
All the relevant data that support the findings of this study are available from the corresponding authors upon reasonable request.