Metalens formed by structured arrays of atomic emitters

Francesco Andreoli; Charlie-Ray Mann; Alexander A. High; Darrick E. Chang

doi:10.1515/nanoph-2024-0603

Article Open Access

Metalens formed by structured arrays of atomic emitters

Francesco Andreoli , Charlie-Ray Mann , Alexander A. High and Darrick E. Chang

Published/Copyright: January 31, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Nanophotonics Volume 14 Issue 3

Abstract

Arrays of atomic emitters have proven to be a promising platform to manipulate and engineer optical properties, due to their efficient cooperative response to near-resonant light. Here, we theoretically investigate their use as an efficient metalens. We show that, by spatially tailoring the (subwavelength) lattice constants of three consecutive two-dimensional arrays of identical atomic emitters, one can realize a large transmission coefficient with arbitrary position-dependent phase shift, whose robustness against losses is enhanced by the collective response. To characterize the efficiency of this atomic metalens, we perform large-scale numerical simulations involving a substantial number of atoms (N ∼ 5 × 10⁵) that is considerably larger than comparable works. Our results suggest that low-loss, robust optical devices with complex functionalities, ranging from metasurfaces to computer-generated holograms, could be potentially assembled from properly engineered arrays of atomic emitters.

Keywords: metalens; metasurface; quantum optics; atomic emitters; color centers; atomic arrays

1 Introduction

Light-mediated dipole–dipole interactions in dense ensembles of atom-like emitters, and the wave interference encoded in them, can lead to a cooperative response that is markedly different from that of an isolated emitter [1], [2]. This resource is most effectively harnessed in ordered arrays of emitters with subwavelength lattice constants, where the collective behavior leads to nontrivial phenomena, including an efficient, directional coupling to light. Capitalizing on these properties, many works have explored classical and quantum optical applications of atomic arrays [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], such as the realization of an atomically thin mirror [23], [24], [25]. Perhaps most relevant to the theme of this paper, these arrays have been proposed to implement various classical optical functionalities, including nonreciprocity [26], optical magnetism [27], [28], [29], wavefront engineering [28], [29], [30], polarization control [31], [32], and chiral sensing [33]. Here, we explore a distinct route toward their application as an optical metalens, which only requires the ability to design the positions of identical emitters.

Metalenses have recently emerged as a promising alternative to traditional bulk optics, enabling complex optical operations while retaining subwavelength thicknesses [34], [35]. Their functionality demands simultaneous control over both transmission intensity and phase pattern. In conventional metasurfaces, this is achieved by spatially varying the size, shape, and orientation of individual nanoscatterers, which generally support both electric and magnetic modes. In contrast, the optical response of atom-like quantum emitters is usually dominated by electric dipole transitions, and it offers limited control over their radiative properties. On the other hand, atomic emitters represent an excellent playground to engineer collective effects, as their electronic transition can provide a low-loss, near-resonant optical resonance, with a large scattering cross section ∼ λ 0 2 , compared to their point-like, physical size [36]. Inspired by the paradigms of conventional metasurfaces, previous works have proposed to engineer an optical metalens out of a bi-layer atomic array, by locally shifting the resonance frequencies of the individual emitters with additional dressing lasers, whose intensities should vary on a subwavelength scale [28], [29], [30]. A similar approach was also proposed in Ref. [37], involving a disordered sheet of atoms.

With one eye on integrated photonic devices, here we propose a different mechanism to realize an efficient metalens, which only requires a suitable choice of the positions of solid-state, atom-like emitters. Specifically, we demonstrate that one can achieve full control of the transmission phase in a bi-layer, rectangular array, while maintaining unit transmittance, by simply varying lattice constants and layer spacing. Moreover, by adding a third layer, we show that these transmission properties can be robustly maintained even in the presence of nonradiative losses or other imperfections, owing to the enhanced collective response. Finally, we demonstrate that these structures can be used as building blocks of an efficient metalens, which we verify through large-scale numerical simulations involving a substantial number of emitters (up to N ∼ 5 × 10⁵), which is considerably higher than comparable works [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49]. The corresponding code is available for public use at Ref. [50], provided with a broader, user-friendly toolbox to simulate the linear optical response of an arbitrary set of two-level, quantum emitters.

The rest of the paper is structured as follows. First, in Section 2, we review the concept of metalenses, and we introduce the physical system under analysis and its theoretical model. Then, in Section 3, we show how arrays of atomic emitters can be engineered to guarantee unit transmission and tunable phase shift. In Section 4, we use these elements to design an illustrative metalens composed of atomic arrays, and in Section 4.1, we test its behavior through extensive numerics, while optimizing its free parameters via a global particle-swarm algorithm [50]. Finally, in Section 4.2, we probe the resistance of that design against different sources of losses or imperfections.

2 Overview of metalens concept and presentation of our system

Conventional refractive lenses rely on local variations of the optical path inside the lens (where light experiences a higher, positive refractive index) to induce a spatially dependent phase shift. Thereby, the wavefront is shaped in such a way that the output beam focuses at a designed distance, as pictorially represented in Figure 1a. In the past couple of decades, however, the novel idea of developing flat metalenses with much smaller footprints has emerged [51], [52], [53], [54], [55]. These metalenses rely on the electromagnetic response of tailored nanostructures to locally impress abrupt phase shifts on the transmitted light [35], [56], [57], [58], while maintaining a thickness on the order of the wavelength or less [34], [35].

$Figure 1: Pictorial comparison between a textbook bulk lens and an atomic metalens. (a) Bulk lens of refractive index n, whose spatially variable optical path d(R) induces a phase delay ϕ lens(R), which curves the incident wavefront, to make it focus at the target distance. (b) Schematic structure of an atomic metalens. Its building blocks consist of at least two atomic arrays in series, whose subwavelength lattice constants d x,y,z < λ 0 can be engineered to ideally ensure a fully directional transmission, with an arbitrary phase shift. For a realistic, lossy system, three atomic layers are required to enhance the robustness to losses.$

Figure 1:

Pictorial comparison between a textbook bulk lens and an atomic metalens. (a) Bulk lens of refractive index n, whose spatially variable optical path d(R) induces a phase delay ϕ _lens(R), which curves the incident wavefront, to make it focus at the target distance. (b) Schematic structure of an atomic metalens. Its building blocks consist of at least two atomic arrays in series, whose subwavelength lattice constants d _x,y,z < λ ₀ can be engineered to ideally ensure a fully directional transmission, with an arbitrary phase shift. For a realistic, lossy system, three atomic layers are required to enhance the robustness to losses.

Regardless of physical implementation, the function of a simple ideal lens of focal length f on a monochromatic input beam of light with wavevector k = ( 2 π / λ 0 ) z ̂ = k 0 z ̂ is to impart the position-dependent phase profile

(1) ϕ lens R = k 0 f − R 2 + f 2 + ϕ 0 ,

upon transmission. This phase is defined modulo 2π, and here we adopt the convention −π ≤ ϕ _lens ≤ π. Moreover, we define the transverse coordinate R = x 2 + y 2 , while the parameter ϕ ₀ corresponds to the phase at the center of the lens [53], [59]. Rather than using dielectric or metallic nano-elements to realize this phase, an atomic metalens instead relies on the use of properly positioned, two-level, solid-state emitters (see Figure 1b).

Although the theory that we present will be rather general, from an experimental perspective color centers in diamond can offer a promising framework for its implementation, as they stand out for their excellent optical properties [60]. Specifically, they behave as atom-like emitters with well-defined selection rules and a dipolar response aligned along one of the four possible tetrahedral directions of the diamond lattice [60], [61], [62], [63], [64], [65]. Current fabrication technologies, moreover, offer good control over their spatial position [66], up to < 10 nm [67], [68]. At the same time, recent works have explored ways to fix the dipole orientations along a well-defined axis [69], [70], [71], [72] or create exactly one emitter at a target position [73], [74]. Although the full combination of these properties into either 2D [75], [76] or 3D [77] large-scale arrays remains a challenge, recent experimental efforts show promising results toward that direction [74]. Concretely, to achieve high-resolution 3D structures, one can envision a process of patterned ion implantation followed by near-field enhanced laser writing, to deterministically write defects at the desired locations. This would be followed by overgrowth of subsequent array layers [78]. Alternatively, processes utilizing block-polymers can also realize the desired resolution [79]. In principle, all these processes could be modified to create desired arrays of point dipoles.

More specifically, we focus on the case of Silicon Vacancy (SiV) centers, which we model as idealized two-level emitters with resonant frequency 2πc/ω ₀ ≈ 737 nm. In this model, we assume that the fabrication process permits to preferentially discriminate over the four possible orientations, so that all the emitters have the same dipole matrix element P 0 = P 0 x ̂ . Moreover, we characterize these emitters with both a coherent, radiative and elastic scattering rate Γ 0 = k 0 3 | P 0 | 2 / ( 3 π ϵ ℏ ) , and an additional broadening Γ′ ≈ 5.75Γ₀ which accounts for losses and other deviations from the ideal case. Here, k ₀ = 2π/λ ₀ = nω ₀/c denotes the resonant wavevector within the bulk diamond of refractive index n ≃ 2.4. Further details on the definition of Γ′ are discussed in Appendix A. Although the ratio Γ 0 / Γ 0 + Γ ′ ≈ 0.15 is relatively low, the optical response of an atomic metalens is protected by the collective behavior, thus allowing for higher efficiencies. To conclude, although we focus on this illustrative level of detrimental broadening, in Sec. IVB, we study the behavior of our system when increasing Γ′ by orders of magnitude.

3 Global control of transmission

We now introduce the theoretical framework to capture the linear optical response of a collection of N quantum emitters in response to a monochromatic classical field, allowing for arbitrary positions. For intensities below the saturation threshold, the nonlinear behavior of a quantum emitter is negligible, and each SiV linearly responds to near-resonant light with a characteristic polarizability α 0 ( Δ , Γ 0 ) = − 3 π Γ 0 / ( Δ + i Γ 0 + Γ ′ / 2 ) k 0 3 , where Δ = ω − ω 0 corresponds to the detuning between the input ω and resonant ω ₀ frequencies [80].

The total field at any point in space consists of the sum between the incident field E _in(r) and the field rescattered by the atomic emitters, reading

(2) E out ( r ) = E in ( r ) + k 0 2 ϵ ∑ j = 1 N G ̄ ̄ ( r − r j ) ⋅ p j ,

where the dyadic Green’s tensor

(3) G ̄ ̄ ( r ) = 1 4 π I ̄ ̄ + ∇ ⊗ ∇ k 0 2 e i k 0 | r | | r | ,

defines the scattering pattern of each atomic dipole p j = p j x ̂ . For simplicity, the Green’s tensor is computed at the resonant frequency ω ₀, making the equations local in time. This approximation is commonly adopted in the context of atomic physics, owing to the small bandwidth of the optical response Γ₀ ≪ ω ₀ [5]. Moreover, this approach becomes exact in the resonant limit of Δ = 0 that will be later considered.

The dipole moments of the emitters are linearly driven by the total field at their position, leading to the self-consistent coupled-dipole equations [81]

(4) p i P 0 = α 0 k 0 3 3 π Ω in ( r i ) Γ 0 + ∑ j ≠ i N − 1 G i j p j P 0 ,

which account for the process of multiple light scattering in a nonperturbative fashion. Here, we defined the parameter G i j = ( 3 π / k 0 ) x ̂ ⋅ G ̄ ̄ r i − r j ⋅ x ̂ , while we introduced the input Rabi frequency Ω in ( r ) = P 0 * ⋅ E in ( r ) / ℏ .

3.1 Transmission of M arrays in series

Our goal is to show how the transverse lattice constants d _x,y and distances d _z of a stack of M ≥ 2, 2D rectangular arrays of atomic emitters can be chosen to impress an arbitrary phase shift, while preserving unit transmission. To do so, it is useful to define the atomic dipoles as p _mj, whose double indices identify the positions as r m j = z m z ̂ + R j , with transverse coordinates R j = x j x ̂ + y j y ̂ .

We first review the cooperative behavior of a single, rectangular 2D array, placed at z = z _m. For simplicity, we assume that the input light is a x ̂ -polarized, plane wave E in ( R , z ) = E 0 e i k 0 z x ̂ , and we focus on the limit where the arrays infinitely extend in the transverse directions x ̂ , y ̂ . Within this regime, any generic solution p m j = ∫ d q x y p m ( q x y ) e i q x y ⋅ R j of Eq. (4) can be written as a superposition of transverse Bloch modes with wavevector q _xy. A plane wave at normal incidence, however, only excites the mode with vanishing transverse wavevector q _xy = 0, meaning that all the dipole moments simplify to p _mj = p _m(q _xy = 0) = p _m. The whole array, then, cooperatively responds to light as a single, collective degree of freedom, with an effective polarizability α _coop = α ₀(Δ − ω _coop, Γ_coop), characterized by the cooperative decay rate Γ_coop(d _x,y) and frequency shift ω _coop(d _x,y) of the excited mode [23], [24]. Physically, these properties come from the single atoms interacting with the fields generated by all the others in the plane; mathematically, when assuming p _mj = p _m in Eq. (4), one obtains the in-plane contribution Γ 0 ∑ j ≠ i G m m i j = − ω coop + i ( Γ coop − Γ 0 ) / 2 , which can be computed with the prescription of Refs. [6], [82], [83].

Once excited, the field coherently scattered by each array can be calculated via Eq. (2). Due to the discrete translational symmetry, the array can add a reciprocal lattice vector k x y ( a , b ) = 2 π ( a x ̂ / d x + b y ̂ / d y ) to the incident field, where a , b ∈ Z are integers. This results in a set of diffraction orders with total wavevector k ( a , b ) = k x y ( a , b ) + k z ( a , b ) z ̂ , where the z-component is k z ( a , b ) = k 0 2 − | k x y ( a , b ) | 2 since energy is conserved |k ^(a,b)| = k ₀. In the relevant subwavelength regime d _x,y < λ ₀, all diffraction orders become evanescent except k z ( 0,0 ) = k 0 . This ensures the selective radiance of the array into the same mode of the input light [84], [85], with a cooperative decay rate Γ coop = 3 Γ 0 λ 0 2 / ( 4 π d x d y ) that scales inversely with the lattice constants, and can thus be significantly greater than the single emitter rate. When stacking M arrays consecutively, the scattered light is then constrained within the normal direction k = k 0 z ̂ , and each array responds with the same polarizability α _coop mentioned before (as pictorially described in Figure 2). At this point, Eq. (4) simplifies into a smaller set of M equations for the dipole amplitudes p _n of each array [16]

(5) p n P 0 = α coop k 0 3 3 π Ω 0 Γ coop e i k 0 z n + ∑ m ≠ n M − 1 G n m p m P 0 ,

where Ω₀ = Ω_in(0, 0), while the terms G n m = G n m rad + G n m ev are related to the field scattered by an array at z _m and probed by the array at z _n. Its radiative part is given by G n m rad = ( i / 2 ) e i k 0 | z m − z n | , while G n m ev is the sum of the evanescent diffraction orders with imaginary wavevectors k z ( a , b ) , whose value is reported in Appendix B.

$Figure 2: 1D, cooperative model for a 3D atomic array, illuminated at normal incidence. We consider a stack of M subwavelength, rectangular 2D arrays of atomic emitters with constant d x,y , separated by a longitudinal distance d z . The emitters are identical two-level systems, with a resonant frequency ω 0 and spontaneous emission rate Γ0, which identify the polarizability α 0(Δ = ω − ω 0, Γ0). The layers are illuminated at normal incidence and can scatter light only in this direction (red, wavy arrows), since the other diffraction orders are evanescent (blue, shaded, wavy arrows). Within each 2D array, the optical response is characterized by a single-mode, collective transition, with cooperative resonant frequency ω coop and decay rate Γ coop = 3 Γ 0 λ 0 2 / ( 4 π d x d y ) ${{\Gamma}}_{\text{coop}}=3{{\Gamma}}_{0}{\lambda }_{0}^{2}/\left(4\pi {d}_{x}{d}_{y}\right)$ , characterizing the cooperative polarizability α coop = α 0(Δ − ω coop, Γcoop).$

Figure 2:

1D, cooperative model for a 3D atomic array, illuminated at normal incidence. We consider a stack of M subwavelength, rectangular 2D arrays of atomic emitters with constant d _x,y, separated by a longitudinal distance d _z. The emitters are identical two-level systems, with a resonant frequency ω ₀ and spontaneous emission rate Γ₀, which identify the polarizability α ₀(Δ = ω − ω ₀, Γ₀). The layers are illuminated at normal incidence and can scatter light only in this direction (red, wavy arrows), since the other diffraction orders are evanescent (blue, shaded, wavy arrows). Within each 2D array, the optical response is characterized by a single-mode, collective transition, with cooperative resonant frequency ω _coop and decay rate Γ coop = 3 Γ 0 λ 0 2 / ( 4 π d x d y ) , characterizing the cooperative polarizability α _coop = α ₀(Δ − ω _coop, Γ_coop).

After solving the set of collective coupled-dipole equations Eq. (5), one can use Eq. (2) to reconstruct the field. Since each array can only selectively radiate into the same mode of the input light, it is straightforward to define the far-field transmission and reflection coefficients [16]

(6) t M L = 1 + i Γ coop 2 Ω 0 ∑ m = 1 M p m P 0 e − i k 0 z m , r M L = i Γ coop 2 Ω 0 ∑ m = 1 M p m P 0 e i k 0 z m .

We notice that these equations can be solved without fixing any value of Ω₀, due to the linearity of the optical response p _m ∝ Ω₀. Similarly, Eq. (5) can be directly solved for the dimensionless ratios p m / P 0 , so that the value of the dipole matrix element P 0 does not have to be specified.

To conclude, for the following calculations, we find it favorable to restrict to a regime where the evanescent fields G n m ev ∼ 0 are negligible. For a subwavelength, rectangular lattice, an approximate rule of thumb that guarantees this condition is that all the diffraction orders are at least exponentially suppressed by a factor ∼1/e², which happens when d _z ≳ d _x,y/π. As discussed in Appendix B, further caution is required when approaching d _y ∼ λ ₀, due to perfect interference effects that make G n m ev nominally diverge in the limit of infinitely extended 2D arrays.

3.2 Phase control

A metalens is typically composed of nanostructures as wide as ≲λ ₀, which transmit the majority of light and impress a tunable phase shift. We now show how the lattice constants of a stack of atomic arrays can be similarly engineered, aiming to use them as the building blocks of an atomic metalens. Hereafter, we define the phase of transmission as ϕ _ML = arg t _ML ∈ (−π, π], and we explicitly focus on the resonant case Δ = 0, although the same method can be extended to near-resonant light.

We begin by considering the simplest scenario, corresponding to a single atomic layer in the lossless regime of Γ′ = 0. The complex value of t _1L depends on the difference between the collective resonance frequency ω _coop(d _x,y) and the frequency of the incoming light, which we fixed to the resonance frequency of a single emitter (i.e., Δ = 0). In principle, this means that the transmission phase ϕ _1L = arg t _1L is itself tunable via the choice of lattice constants d _x,y. Nonetheless, using Eqs. (5) and (6), it is easy to show that high transmission and arbitrary phase cannot be achieved with one layer of atoms, as the conditions of reciprocity r 1 L = ( t 1 L − 1 ) e 2 i k 0 z m and energy conservation |t _1L|² + |r _1L|² = 1 impose |t _1L| = cos(ϕ _1L), which limits the phase range to |ϕ _1L| ≤ π/2 and allows unit transmission only in the trivial case of far-detuned driving, where no phase is imprinted ϕ _1L = 0. On the contrary, the largest phase shifts |ϕ _1L| ∼ π/2 are obtained near resonance, where the transmittance drops sharply to zero (i.e., the input field is strongly reflected). Moreover, the range of achievable phases is particularly fragile to the addition of small losses Γ′/Γ_coop ≪ 1, decreasing as | ϕ 1 L | ≲ π / 2 − 2 Γ ′ / Γ coop .

For an ideal system, we can achieve perfect transmission with arbitrary phase by considering a bi-layer (M = 2) array. As long as G m j ev ∼ 0 , this system is equivalent to a Fabry–Perot cavity, composed of two atomic mirrors with the complex reflectivity r _1L and transmission t _1L mentioned above [33], [86]. In the lossless regime Γ′ = 0, it is well known that such an interferometer ensures unit transmission t _2L = exp(2iϕ _1L) when the distance between the mirrors matches the Airy condition k ₀ d _z = πl − arg(r _1L), with l ∈ N [87]. Due to this reason, a proper choice of d _x,y,z allows to keep unit transmission while arbitrarily designing the total phase ϕ _2L = 2ϕ _1L over the full (−π, π] range. This property is represented in Figure 3a, where we independently vary both the subwavelength lattice constants d _x,y and layer spacing d _z, plotting ϕ _2L as a function of d _z and the single-layer parameter ω _coop(d _x,y)/Γ_coop(d _x,y). As expected, we observe full phase tunability with sufficient transmittance, as quantified by the nonshaded, brightly colored regions where |t _2L|² > 0.5.

Figure 3:

Transmission of a multilayer atomic array, as a function of ω _coop(d _x,y)/Γ_coop(d _x,y) and d _z. (a, b) Colorbar representation of the phase shift ϕ _2L = arg t _2L of two atomic layers, given either Γ′ = 0 (a) or Γ′ = 5.75Γ₀ (b). The transverse lattice constants are varied within the range λ ₀ > d _x,y ≥ d _min = 0.03λ ₀, which means that Γ′/Γ_coop(d _x,y) ≳ 0.03. When different choices of d _x and d _y are associated to the same value of ω _coop(d _x,y)/Γ_coop(d _x,y), the pair with the highest cooperative decay is selected. The region where |t _2L|² < 0.5 is represented by a white shaded area, while the insets show the relevant case of ϕ _2L ≡ arg t _2L ∼ ±π and |t _2L|² ≥ 0.5. (c, d) Same structure of subfigures (a) and (b), but for the three-layer case. The white dashed lines represent the chosen branch d _z(d _x,y) that maximizes the transmittance. Along this path, the insets show that both the phase ϕ _3L = ±π and the transmission |t _3L|² ≥ 0.5 can be simultaneously obtained over a much broader bandwidth (c), becoming more resistant to the losses (d).

However, the phase range contracts as | ϕ 2 L | = 2 | ϕ 1 L | ≲ π − 4 Γ ′ / Γ coop in presence of small losses, preventing the achievement of |ϕ _2L| = π, regardless of how small Γ′ > 0 is. As shown by the inset of Figure 3a, this can be related to the asymptotically small bandwidth associated to both |ϕ _2L| ∼ π and |t _2L|² ≥ 0.5 [86], [88], which makes the system more fragile against Γ′. To better quantify this statement, we must first set a minimum interatomic distance d _min = 10 nm, whose value is inspired by the discussion of Section 2. This translates into d _x,y,z ≥ d _min = 0.03λ ₀, which prevents the cooperative response Γ coop ∝ Γ 0 λ 0 2 / ( d x d y ) 2 to become arbitrarily large and overtake any sources of broadening Γ′ > 0. In Figure 3b, we then use the conventional value Γ′ = 5.75Γ₀ of Section 2, observing that both a sufficient transmission |t _2L|² ≥ 0.5 and full phase control can no longer be simultaneously achieved.

In general, M − 1 transparency conditions d _z(d _x,y) similar to that of a Fabry–Perot cavity can be found for arbitrary values of M [89], and the addition of more atomic layers M > 2 is important to restore the resistance to losses around |ϕ _ML| ∼ π. This can be intuitively understood for even number of layers M, as a proper choice of d _z(d _x,y) can make the system act as M/2 cascaded cavities, so that | ϕ M L | = M | ϕ 1 L | ≲ π M / 2 − 2 M Γ ′ / Γ coop . For odd number of layers M, less intuitive conditions for perfect interference hold, but still we show that M = 3 layers are enough to provide resistance to losses.

To define the proper relations d _z(d _x,y) between the longitudinal and transverse lattice constants, we introduce a closed-form solution of Eq. (6), which reads [90]

(7) t M L = e i ( 1 − M ) k 0 d z t 1 L u M ( k , d z ) − u M − 1 ( k , d z ) e i k 0 d z t 1 L ,

where the function u _M(k, d _z) = sin(Mkd _z)/sin(kd _z) relates the finite-size behavior to the dispersion relation k(d _x,y,z) = k(ω _coop(d _x,y)/Γ_coop(d _x,y), d _z) of an infinite chain [16]. In the lossless regime of Γ′ = 0, the unit transmission t _ML = (−1)^a exp(iMk ₀ d _z) is ensured by fixing d _z(d _x,y) to fulfill k(d _x,y,z) = aπ/(Md _z), where the natural number a = 1, …, M − 1 identifies the M − 1 possible solutions within the first Brillouin zone. With this choice, the field acquires a total phase shift of ϕ _ML(d _x,y) = Mk ₀ d _z(d _x,y) + aπ with respect to propagation in the bulk environment.

In our M = 3 case, we choose the branch of d _z(d _x,y) with a = 2, as represented in Figure 3c and 3d by a dashed, white line. When spanning d _x,y, this is associated to high transmittance and complete phase control, in both the lossless (Figure 3c) and lossy Γ′ = 5.75Γ₀ (Figure 3d) regimes. More specifically, we scan the transverse lattice constants d _x,y along the two straight lines (d _x = d _min) ∪ (d _min ≤ d _y < λ ₀) and (d _min ≤ d _x < λ ₀) ∪ (d _y = d _min), which allows to associate a unique set of spacings d _x,y,z to any value of ϕ _3L(d _x,y) = arg t _3L(d _x,y, d _z(d _x,y)). This correspondence is represented in Figure 4a, showing that only a limited set of distances λ ₀/6 ≤ d _z ≤ λ ₀/3 is required, thus implying a maximum thickness of 2 d z max = 2 λ 0 / 3 , which translates to ≈ 205 nm for the case of SiV centers. To conclude, in Figure 4b, we explicitly prove that this scheme allows, in presence of broadening Γ′ = 5.75Γ₀, to maintain a sufficient transmittance |t _3L|² > 0.6 for any relevant value of ϕ _3L.

$Figure 4: Lattice constants d x,y,z and transmittance |t 3L|2 as a function of phase ϕ 3L, given Γ′ = 5.75Γ0. (a) We scan the transverse lattice constants along the two straight lines (d x = d min) ∪ (d min ≤ d y < λ 0) and (d min ≤ d x < λ 0) ∪ (d y = d min), with d min = 0.03λ 0 (black, dashed line). At the same time, the choice of d z (d x,y ) that maximizes the transmittance allows to associate a unique set of lattice constants (colored lines) to any phase ϕ 3L = arg t 3L(d x,y , d z (d x,y )) (horizontal axis). (b) Transmittance |t 3L|2 as a function of the phase ϕ 3L (gray line). The colored points are associated to the rings composing the illustrative atomic metalens discussed in Section 4. Their colors are associated to the relative power of the input light over their area, i.e., P in j ∝ ∫ R j − 1 R j | E in | 2 d R ${P}_{\text{in}}^{j}\propto {\int }_{{R}_{j-1}}^{{R}_{j}}\vert {\mathbf{E}}_{\text{in}}{\vert }^{2}\text{d}\mathbf{R}$ .$

Figure 4:

Lattice constants d _x,y,z and transmittance |t _3L|² as a function of phase ϕ _3L, given Γ′ = 5.75Γ₀. (a) We scan the transverse lattice constants along the two straight lines (d _x = d _min) ∪ (d _min ≤ d _y < λ ₀) and (d _min ≤ d _x < λ ₀) ∪ (d _y = d _min), with d _min = 0.03λ ₀ (black, dashed line). At the same time, the choice of d _z(d _x,y) that maximizes the transmittance allows to associate a unique set of lattice constants (colored lines) to any phase ϕ _3L = arg t _3L(d _x,y, d _z(d _x,y)) (horizontal axis). (b) Transmittance |t _3L|² as a function of the phase ϕ _3L (gray line). The colored points are associated to the rings composing the illustrative atomic metalens discussed in Section 4. Their colors are associated to the relative power of the input light over their area, i.e., P in j ∝ ∫ R j − 1 R j | E in | 2 d R .

We notice that those phases within the interval of 0 ≲ ϕ _3L ≲ 0.01π cannot be engineered, due to the limited value of max ω _coop ≈ 28Γ_coop for d _min ≤ d _x,y ≤ λ ₀. Nonetheless, for practical applications such as a metalens, this range can be approximated with exactly ϕ _3L = 0 (i.e., no emitters), as its span is negligible compared to typical discretization scales.

4 Atomic metalens

To design an atomic metalens out of three-layer atomic arrays, one needs to spatially tune the lattice constants d _x,y, to make the phase shift ϕ _3L(d _x,y) match that of an ideal lens, i.e., the value ϕ _lens(R) specified in Eq. (1). To define a concrete recipe, we divide the transverse plane into concentric rings j = 1, 2… of radius R _j = jΔR (see Figure 5a), and we associate to each ring the central phase shift ϕ _j ≡ ϕ _lens(R _j−1/2 + R _j/2), by using Eq. (1). Here, we recall that the initial phase ϕ ₀ is a free parameter. At this point, we impose ϕ 3 L d x , y j = ϕ j and extract the lattice constants d x , y j by numerically inverting the solid line of Figure 4a. The transparency condition of Figure 4b can then be used to define the longitudinal constant d z j = d z d x , y j . The final metalens is then the union of these discrete building blocks, as shown in Figure 5. By choosing ΔR ≲ λ ₀, we ensure a discretization scale with the same order of magnitude of that of usual metalenses [34].

$Figure 5: Structure of an atomic metalens, with focal length f = 20λ 0 and radius R lens = 10λ 0. (a) 3D representation of the atomic metalens, where each point depicts the position of one atom. This atomic metalens is composed of 15 concentric rings of thickness ΔR ≈ 2λ 0/3, with a buffer-zone parameter α ≈ 0.2. The lens has a width of Δz ≈ 2λ 0/3, much thinner than the total diameter of 20λ 0. The atoms belonging to the j-th ring have the same lattice constants d x , y , z j ${d}_{x,y,z}^{j}$ , which are uniquely associated to the phase shift ϕ j = ϕ lens(ΔR(2j − 1)/2) of Eq. (1) (with ϕ 0 ≃ −2.06), through the curves ϕ j = ϕ 3 L d x , y j ${\phi }_{j}={\phi }_{3\text{L}}\left({d}_{x,y}^{j}\right)$ and d z j = d z d x , y j ${d}_{z}^{j}={d}_{z}\left({d}_{x,y}^{j}\right)$ shown in Figure 4. The color of the atoms in each ring reflects the value of ϕ j , as described by the colorbar at the bottom. (b) Focusing of a x ̂ $\hat{\mathbf{x}}$ -polarized, resonant, input Gaussian beam with w 0 = 4λ 0, by the action of the atomic metalens. The orange, shaded area shows the textbook beam waist w(z) during the focusing process. The metalens is designed to focus the beam at a distance z f ≃ 17λ 0. This defines the focal plane, where we numerically reconstruct the total relative intensity |E out(R, z f )/E 0|2 via the input–output formalism of Eq. (2), in the lossy regime of Γ′ = 5.75Γ0. The value of |E out(R, z f )/E 0|2 is portrayed with the color scheme shown by the colorbar at the bottom. Further results from the coupled-dipole simulations are shown in Figure 6.$

Figure 5:

Structure of an atomic metalens, with focal length f = 20λ ₀ and radius R _lens = 10λ ₀. (a) 3D representation of the atomic metalens, where each point depicts the position of one atom. This atomic metalens is composed of 15 concentric rings of thickness ΔR ≈ 2λ ₀/3, with a buffer-zone parameter α ≈ 0.2. The lens has a width of Δz ≈ 2λ ₀/3, much thinner than the total diameter of 20λ ₀. The atoms belonging to the j-th ring have the same lattice constants d x , y , z j , which are uniquely associated to the phase shift ϕ _j = ϕ _lens(ΔR(2j − 1)/2) of Eq. (1) (with ϕ ₀ ≃ −2.06), through the curves ϕ j = ϕ 3 L d x , y j and d z j = d z d x , y j shown in Figure 4. The color of the atoms in each ring reflects the value of ϕ _j, as described by the colorbar at the bottom. (b) Focusing of a x ̂ -polarized, resonant, input Gaussian beam with w ₀ = 4λ ₀, by the action of the atomic metalens. The orange, shaded area shows the textbook beam waist w(z) during the focusing process. The metalens is designed to focus the beam at a distance z _f ≃ 17λ ₀. This defines the focal plane, where we numerically reconstruct the total relative intensity |E _out(R, z _f)/E ₀|² via the input–output formalism of Eq. (2), in the lossy regime of Γ′ = 5.75Γ₀. The value of |E _out(R, z _f)/E ₀|² is portrayed with the color scheme shown by the colorbar at the bottom. Further results from the coupled-dipole simulations are shown in Figure 6.

At the interface between the finite rings, the abrupt change of lattice constants can potentially scatter light into unwanted diffraction modes. To soften these detrimental effects, in the x ̂ , y ̂ -plane we introduce a small buffer zone between two consecutive rings, with atoms placed at intermediate positions. These zones extend over the first fraction 0 ≤ α ≤ 1/2 of each ring, and their definition is not strict, with many possible variants. Our approach is described in Appendix C, and we numerically associate it to a small efficiency increase, up to an additional factor ∼ 0.02 in the estimated efficiency.

To conclude, we remark that for each target focal length f, our atomic metalens is defined up to three free parameters, which are an overall phase shift −π < ϕ ₀ ≤ π, the ring thickness d _min ≪ ΔR ≲ λ ₀, and the buffer fraction 0 ≤ α ≤ 1/2.

4.1 Numerical simulations

To check our design, we want to estimate the efficiency of an atomic metalens with focal length f and centered around z = 0. To this aim, we fix the atomic positions, and we illuminate the system at normal incidence with a x ̂ -polarized, resonant, input Gaussian beam focused at z = 0, which has beam waist w ₀ and focal intensity |E ₀|² (see Appendix D). We then perform exact simulations of the linear optical response, reconstructing the total field E _out(R, z) via Eqs. (2) and (4). We want to compare it with the theoretical prediction of the field transmitted by an ideal, thin lens of focal length f. This is given by the Gaussian beam E _f(R, z), characterized by the beam waist w f = w 0 / M , the focal position z f = ( 1 − M − 2 ) f , and the focal intensity | E f ( 0 , z f ) / E 0 | 2 = M 2 . Here, the parameter M = 1 + k 0 w 0 2 / ( 2 f ) 2 ≥ 1 is the so-called magnification of the lens, which quantifies the focusing ability and ensures the conservation of energy ∫|E _f|²dR = ∫|E _in|²dR ∝ P _in.

To characterize the metalens performance, we quantify the fraction η = P _η/P _in of power P _η that is correctly transmitted into the target, ideal Gaussian mode E _f, divided by the total input power P _in [91]. Operatively, this efficiency can be obtained by analytically projecting E _out into the target mode E _f, namely η = |⟨E _f|E _out⟩|². This projection has a simple, closed-form expression, which is detailed in Appendix D. Another quantity of interest is the overlap between the transmitted field and the input field ϵ = |⟨E _in|E _out⟩|². Obviously, one would aim to operate in a regime where η ∼ 1, while ϵ ≪ 1, with the latter inequality signifying that the lens performs some non-negligible transformation. Finally, we notice that, for certain applications, the main requirement is the identification of the focal spot over the background of transmitted light. In view of that, we define the signal-to-background ratio η ̃ = P η / P t , which divides the power transmitted into the target mode P _η by the total transmitted power P _t, rather than by the total input power. Here, one has P _η = ηP _in, while P _t ∝ ∫|E _out|²dR is numerically computed from the total field at the focal plane z = z _f.

To show the potential of our scheme, we can now discuss an illustrative full-scale simulation of a metalens with focal length f = 20λ ₀ and radius R _lens = 10λ ₀, illuminated by an input Gaussian beam of waist w ₀ = 4λ ₀. In this illustrative scenario, the ideal magnification would read M = w 0 / w f ≃ 2.7 , associated to an ideal intensity enhancement of | E f ( 0 , z f ) / E 0 | 2 = M 2 ≃ 7.32 . These simulations involve a substantial number of atoms N ∼ 5 × 10⁵, and the techniques by which we accomplish this result are described in the Methods. All the codes are written in Julia [92] and are available at Ref. [50]. The free parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≃ −2.06, and α ≈ 0.2 are chosen to maximize η in the lossy regime Γ′ = 5.75Γ₀. This was first accomplished via a brute-force optimization and then confirmed through a particle-swarm, global-optimization algorithm [50].

The numerical results are shown in Figure 6, where we plot the relative intensity of the total field |E _out(R, z)/E ₀|², calculated on the horizontal plane y = 0 (top row) and at the expected focal plane z = z _f ≃ 17λ ₀ (bottom row). The column on the left (Figure 6a and 6d) shows the ideal values that one would expect for a textbook, ideal lens, i.e., E _f(R, z). This is compared to the numerical simulations of the atomic metalens, calculated for the lossy case Γ′ = 5.75Γ₀ (right column, Figure 6b and 6e). Very similar plots are obtained when studying the lossless case Γ′ = 0, or when plotting the intensity on the plane x = 0.

$Figure 6: Illustrative case of an atomic metalens with focal length f = 20λ 0, radius R lens = 10λ 0, and parameters ΔR ≈ 2λ 0/3, ϕ 0 ≈ −2.06, and α ≈ 0.2, illuminated by a resonant Gaussian beam with waist w 0 = 4λ 0. The figures show the relative intensity of the total field |E out(R, z)/E 0|2, calculated on the planes y = 0 (top row, subfigures a, b) and z = z f ≃ 17λ 0 (bottom row, subfigures d, e). The subplots (a, d) represent the ideal case of a textbook lens, while the subplots (b, e) show the results of the numerical simulations with Γ′ = 5.75Γ0. The dashed, white lines represent the ideal value of the beam waist w(z), while the dot-dashed, white lines show the waist of the input beam if no lens were present. The efficiency of the lossy Γ′ = 5.75Γ0 case, estimated from the simulations, reads η ≃ 0.90, while the signal-to-background ratio reads η ̃ > 0.98 $\tilde {\eta }{ >}0.98$ . The number of simulated atoms is N ≃ 4.6 × 105. Finally, the subplots (c, f) show two line-cuts of the intensity profile, along either the z axis in the x = y = 0 plane (c) or the x axis in the y = 0, z = z f plane (f). This quantity is depicted for both the ideal (dashed, red line) and lossy (solid, blue line) cases. The gray area in subfigure (c) depicts the space occupied by the atomic metalens.$

Figure 6:

Illustrative case of an atomic metalens with focal length f = 20λ ₀, radius R _lens = 10λ ₀, and parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≈ −2.06, and α ≈ 0.2, illuminated by a resonant Gaussian beam with waist w ₀ = 4λ ₀. The figures show the relative intensity of the total field |E _out(R, z)/E ₀|², calculated on the planes y = 0 (top row, subfigures a, b) and z = z _f ≃ 17λ ₀ (bottom row, subfigures d, e). The subplots (a, d) represent the ideal case of a textbook lens, while the subplots (b, e) show the results of the numerical simulations with Γ′ = 5.75Γ₀. The dashed, white lines represent the ideal value of the beam waist w(z), while the dot-dashed, white lines show the waist of the input beam if no lens were present. The efficiency of the lossy Γ′ = 5.75Γ₀ case, estimated from the simulations, reads η ≃ 0.90, while the signal-to-background ratio reads η ̃ > 0.98 . The number of simulated atoms is N ≃ 4.6 × 10⁵. Finally, the subplots (c, f) show two line-cuts of the intensity profile, along either the z axis in the x = y = 0 plane (c) or the x axis in the y = 0, z = z _f plane (f). This quantity is depicted for both the ideal (dashed, red line) and lossy (solid, blue line) cases. The gray area in subfigure (c) depicts the space occupied by the atomic metalens.

We benchmark the optical response of the atomic metalens from our simulations, finding an efficiency η ≃ 0.95 and an intensity enhancement at the focal point of |E _out(0, z _f)/E ₀|² ≃ 6.03, in the lossless regime of Γ′ = 0. Similarly, in the lossy case of Γ′ = 5.75Γ₀, we obtain the values η ≃ 0.90 and |E _out(0, z _f)/E ₀|² ≃ 5.60. This value can be appreciated in Figure 6c and 6f, where we compare the ideal (red, dashed line) and numerical (blue, solid line) field intensity along, respectively, either the z axis (in the x = y = 0 plane) or the x axis (at the focal plane). These high efficiencies stand out when considering the much lower overlap ϵ ≃ 0.42 between the output field and the input beam, which means that the atomic metalens is nontrivially acting on the input beam. Finally, both the lossy and the lossless cases exhibit a high signal-to-background ratio, reading η ̃ > 0.98 . To understand how the broadening Γ′ = 5.75Γ₀ affects the efficiency, we recall from Figure 4b that the transmittance |t _3L|² highly depends on ϕ _3L, meaning that some rings can transmit more light than others. Considering our illustrative metalens, the complex transmission associated to each ring is represented with a colored point in Figure 4b. The overall reduction of the efficiency due to the losses (i.e., the ratio between the lossy Γ′ > 0 and lossless Γ′ = 0 efficiencies) agrees well with the average transmittance |t _3L(ϕ _j)|² of the rings, each weighted by the relative power of the input light illuminating their area (corresponding to the color of the points in Figure 4b). Notably, this intuitive model explains why the efficiency η can strongly depend on the choice of ϕ ₀.

Although the atomic metalens was designed to operate for resonant light at Δ = 0, a similar reasoning allows to qualitatively predict the spectral bandwidth where the efficiency remains high. To show this, we calculate the cooperative decay rates Γ coop j for all the rings that compose the metalens and weight them by the corresponding fraction of input light, to define the average value ⟨ Γ coop j ⟩ ≈ 96 Γ 0 (of the order of ∼ 2 π × 10 GHz for SiVs [93], [94]). As detailed in Appendix E, we observe that the efficiency remains as high as η ≳ 0.8 as long as | Δ | ≤ ⟨ Γ coop j ⟩ / 2 , while quickly decreasing outside.

To conclude, it is interesting to investigate how the response is modified when increasing the focusing ability of the lens, as quantified by the magnification M . Specifically, in Figure 7, we fix w ₀ = 4λ ₀ and scan different focal lengths f, plotting the efficiency η (blue points) and the signal-to-background ratio η ̃ (green points) as a function of 1 ≤ M ≲ w 0 / λ 0 ≪ k 0 w 0 . Here, the maximum magnification is associated to the limit k ₀ w _f ≫ 1 imposed by the paraxial approximation, while the choice of w ₀ = 4λ ₀ represents the largest beam waist that we can compute, due to the numerical complexity of the simulation. In presence of broadening Γ′ = 5.75Γ₀, we observe that the efficiency remains as high as η ≳ 0.82 (dotted, black line) up to M = 4 , where the overlap with the input field is as low as ϵ ≈ 0.26. Overall, we find the empirical scalings of η ≈ 1.06 − 0.06 M , η ̃ ≈ 1.05 − 0.03 M , and ϵ ≈ − 0.04 + 1.23 / M (colored dashed lines). Assuming that these scalings would hold true for larger values of w ₀, they would predict efficiencies as high as η ≈ 0.5 up to M ≈ 10 (where the overlap with the input field is as low as ϵ ≈ 0.08), and signal-to-background ratios larger than η ̃ ≳ 0.5 up to M ≈ 20 (where ϵ ≈ 0.02).

$Figure 7: Efficiency of an atomic metalens as a function of the magnification, given Γ′ = 5.75Γ0. We fix both the waist of the input beam to w 0 = 4λ 0, and the radius of the lens R lens = 10λ 0, while showing the efficiency η (blue points), signal-to-background ratio η ̃ $\tilde {\eta }$ (green points), and input-field overlap ϵ (orange points) as a function of the magnification 1 ≤ M ≲ w 0 / λ 0 = 4 $1\le \mathcal{M}< sim {w}_{0}/{\lambda }_{0}=4$ . For each point, we perform a particle-swarm optimization of the free parameters ϕ 0, α, and ΔR to maximize the efficiency η [50]. By fitting the data, we infer the empirical scalings η ≈ 1.06 − 0.06 M $\eta \approx 1.06-0.06\mathcal{M}$ , η ̃ ≈ 1.05 − 0.03 M $\tilde {\eta }\approx 1.05-0.03\mathcal{M}$ and ϵ ≈ − 0.04 + 1.23 / M ${\epsilon}\approx -0.04+1.23/\mathcal{M}$ (colored, dashed lines). The black, dotted line shows the reference value of 0.8.$

Figure 7:

Efficiency of an atomic metalens as a function of the magnification, given Γ′ = 5.75Γ₀. We fix both the waist of the input beam to w ₀ = 4λ ₀, and the radius of the lens R _lens = 10λ ₀, while showing the efficiency η (blue points), signal-to-background ratio η ̃ (green points), and input-field overlap ϵ (orange points) as a function of the magnification 1 ≤ M ≲ w 0 / λ 0 = 4 . For each point, we perform a particle-swarm optimization of the free parameters ϕ ₀, α, and ΔR to maximize the efficiency η [50]. By fitting the data, we infer the empirical scalings η ≈ 1.06 − 0.06 M , η ̃ ≈ 1.05 − 0.03 M and ϵ ≈ − 0.04 + 1.23 / M (colored, dashed lines). The black, dotted line shows the reference value of 0.8.

4.2 Losses and imperfections

Up to now, the presence of experimental losses and imperfections has been modeled by the addition of a detrimental broadening Γ′ ≈ 5.75Γ₀, whose value was chosen to qualitatively capture some key properties of state-of-the-art experiments with color centers in diamond. While our studies up to now represent an optimistic scenario, here we investigate the performance of the metalens as the broadening rate Γ′ increases, or when the atoms are subject to increasing spatial disorder.

First, we study the resistance to increasing levels of broadening Γ′, which we compare with the maximum cooperative decay rate Γ coop max = Γ coop ( d x , y = d min ) ≈ 225 Γ 0 allowed in the system. To this aim, it is instructive to focus on the single building blocks of the metalens. In Figure 8a, we show the relation between the phase ϕ _3L (on the horizontal axis) and transmittance |t _3L|² (color scheme), when considering increasing values of Γ′ (vertical axis, in log scale). This corresponds to the extension of Figure 4b (which coincides with the black dashed line in Figure 8a) to arbitrary values of Γ′. Notably, when Γ ′ ≳ 0.15 Γ coop max ≈ 30 Γ 0 , some phases cannot be realized anymore (black areas in the plot). We recall that the addition of further atomic layers is expected to drastically increase the resistance to losses, although presenting the drawback of adding more atomic emitters, and increasing the overall thickness of the metalens. Reducing the minimum lattice constant d _min would similarly work, by increasing the maximum cooperative rate Γ coop max .

$Figure 8: Resistance to nonradiative losses. (a) Transmission of a three-layer array, given increasing levels of Γ′. Similarly to Figure 4b, we use our definition of d x,y,y to associate a unique transmittance |t 3L|2 (color scheme) to any target phase ϕ 3L (horizontal axis). We then vary Γ′ (vertical axis) to track the change in the transmittance. We notice that an almost identical plot is obtained when numerically optimizing the choices of d x,y,z ≥ d min = 0.03λ 0 to maximize transmittance, proving the validity of our scheme. The black, dashed line highlights the particular case Γ′ = 5.75Γ0. The black areas (bounded by dotted, white lines) identify regions of the parameter space that cannot be obtained with any choice of d x,y,z . (b) Efficiency as a function of Γ′, given an atomic metalens with focal length f = 20λ 0, radius R lens = 10λ 0, and construction parameters ΔR ≈ 2λ 0/3, ϕ 0 ≈ −2.06, and α ≈ 0.2, illuminated by a Gaussian beam with w 0 = 4λ 0. The lines show the efficiency η (blue), signal-to-background ratio η ̃ $\tilde {\eta }$ (green), input-field overlap ϵ (orange), and base-line efficiency η in = |⟨E f |E in⟩|2 (black, dashed line). The colored, dotted lines represent the values at Γ′ = 0, while the colored points show the case of Γ′ = 5.75Γ0. The black, dotted line depicts a threshold value of 0.9, while the shaded, gray region portrays the regime where some phases cannot be engineered anymore, corresponding to the appearance of black areas in subfigure (a). Finally, the blue asterisks show the efficiencies in case the structural parameters ΔR, ϕ 0, and α are changed to be optimal for the corresponding value of Γ′.$

Figure 8:

Resistance to nonradiative losses. (a) Transmission of a three-layer array, given increasing levels of Γ′. Similarly to Figure 4b, we use our definition of d _x,y,y to associate a unique transmittance |t _3L|² (color scheme) to any target phase ϕ _3L (horizontal axis). We then vary Γ′ (vertical axis) to track the change in the transmittance. We notice that an almost identical plot is obtained when numerically optimizing the choices of d _x,y,z ≥ d _min = 0.03λ ₀ to maximize transmittance, proving the validity of our scheme. The black, dashed line highlights the particular case Γ′ = 5.75Γ₀. The black areas (bounded by dotted, white lines) identify regions of the parameter space that cannot be obtained with any choice of d _x,y,z. (b) Efficiency as a function of Γ′, given an atomic metalens with focal length f = 20λ ₀, radius R _lens = 10λ ₀, and construction parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≈ −2.06, and α ≈ 0.2, illuminated by a Gaussian beam with w ₀ = 4λ ₀. The lines show the efficiency η (blue), signal-to-background ratio η ̃ (green), input-field overlap ϵ (orange), and base-line efficiency η _in = |⟨E _f|E _in⟩|² (black, dashed line). The colored, dotted lines represent the values at Γ′ = 0, while the colored points show the case of Γ′ = 5.75Γ₀. The black, dotted line depicts a threshold value of 0.9, while the shaded, gray region portrays the regime where some phases cannot be engineered anymore, corresponding to the appearance of black areas in subfigure (a). Finally, the blue asterisks show the efficiencies in case the structural parameters ΔR, ϕ ₀, and α are changed to be optimal for the corresponding value of Γ′.

To get further insights, it is instructive to explicitly focus on the illustrative atomic metalens of Figure 6, with focal length f = 20λ ₀, radius R _lens = 10λ ₀, and parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≃ −2.06, and α ≈ 0.2. In Figure 8b, we discuss the overall response of this metalens, for broadening levels up to Γ ′ ≃ 3 × 1 0 2 Γ coop max ≈ 1 0 5 Γ 0 . The blue line depicts the efficiency η, the orange line the input-field overlap ϵ, and the green line the signal-to-background ratio η ̃ . Roughly, the system becomes ineffective above the threshold Γ ′ ≳ ⟨ Γ coop j ⟩ ≈ 0.5 Γ coop max ≈ 1 0 2 Γ 0 . Notably, the efficiency remains acceptable η ≳ 0.7 as long as Γ′ ≳ 60Γ₀ although, in principle, this corresponds to a regime where some phases around |ϕ| ∼ π cannot be engineered anymore (gray, shaded region). At the same time, the signal-to-background ratio η ̃ remains relatively high up to much higher losses, so that η ̃ ≳ 0.9 up to Γ ′ ≈ 0.8 Γ coop max ≈ 1 0 2 Γ 0 and η ̃ ≳ 0.5 up to Γ ′ ≈ 5 Γ coop max ≈ 1 0 3 Γ 0 . We note that these efficiencies are calculated for a fixed choice of ΔR, ϕ ₀, and α, which are optimal only for Γ′ = 5.75Γ₀. This reasoning well describes a situation where the amount of losses is unknown. On the other hand, higher efficiencies (blue asterisk in Figure 8b) are obtained by choosing optimal parameters tailored on the broadening Γ′, as computed via particle-swarm optimization [50].

Finally, we discuss the effect of disorder in the atomic positions, defined by randomly displacing each atomic emitter inside a 3D sphere of radius δd, with a uniform distribution. In Figure 9, we represent with colored points the same quantities of Figure 8b, as a function of increasing disorder δd. As intuitively expected, when the displacement is comparable to d _min, then the efficiency is strongly undermined, with η ∼ 0. In that regime, the transmitted light is so randomly altered, that it does not overlap anymore with the input field either, and one gets ϵ ∼ 0. Nonetheless, we notice that the signal-to-background ratio exhibits more robust properties, with η ̃ ≳ 0.6 up to δd ∼ 0.7d _min. We relate these results to the overall drop of transmitted light that occurs in the disordered regime.

$Figure 9: Resistance to additional position disorder. The data are calculated for the atomic metalens with focal length f = 20λ 0, radius R lens ≃ 9λ 0, and construction parameters ΔR ≈ 2λ 0/3, ϕ 0 ≈ −2.06, and α ≈ 0.2, illuminated by a Gaussian beam with w 0 = 4λ 0. The horizontal axis represents the random displacement radius δd in units of the minimum lattice constant d min. The points represent the average efficiency (η, blue), signal-to-background ratio ( η ̃ $\tilde {\eta }$ , green), and overlap with the input beam (ϵ, orange). Each point is calculated by averaging over 10 random configurations, and the error bars represent one standard deviation. The simulation is performed for the lossy case Γ′ = 5.75Γ0. The lines represent the theoretical prediction when replacing the random displacement with the additional inelastic rate ∼ 2.5 Γ dis ′ ( δ d , d min ) = 2.5 ( π / 2 ) ( δ d / d min ) 2 Γ coop max $\sim 2.5{{\Gamma}}_{\text{dis}}^{\prime }\left(\delta d,{d}_{\mathrm{min}}\right)=2.5\left(\pi /2\right){\left(\delta d/{d}_{\mathrm{min}}\right)}^{2}{{\Gamma}}_{\text{coop}}^{\mathrm{max}}$ , where the numerically inferred prefactor stems from the additional complexity of the metalens, compared to stacks of infinite arrays.$

Figure 9:

Resistance to additional position disorder. The data are calculated for the atomic metalens with focal length f = 20λ ₀, radius R _lens ≃ 9λ ₀, and construction parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≈ −2.06, and α ≈ 0.2, illuminated by a Gaussian beam with w ₀ = 4λ ₀. The horizontal axis represents the random displacement radius δd in units of the minimum lattice constant d _min. The points represent the average efficiency (η, blue), signal-to-background ratio ( η ̃ , green), and overlap with the input beam (ϵ, orange). Each point is calculated by averaging over 10 random configurations, and the error bars represent one standard deviation. The simulation is performed for the lossy case Γ′ = 5.75Γ₀. The lines represent the theoretical prediction when replacing the random displacement with the additional inelastic rate ∼ 2.5 Γ dis ′ ( δ d , d min ) = 2.5 ( π / 2 ) ( δ d / d min ) 2 Γ coop max , where the numerically inferred prefactor stems from the additional complexity of the metalens, compared to stacks of infinite arrays.

As detailed in Appendix A.2, small displacements in a 2D array (or in a stack of arrays) can be well described by a supplementary broadening Γ dis ′ ( δ d , d x , y ) ≈ π δ d 2 / ( 2 d x d y ) Γ coop ( d x , y ) , whose scaling ensures the dependence of the optical response only on the relative displacement δ d / d x d y . For the more complex case of an atomic metalens, we numerically find that the position disorder can be still characterized by a supplementary rate ∼ 2.5 Γ dis ′ ( δ d , d min ) , where the empirical prefactor can be attributed to the more fragile interference patterns involved in the metalens response, as well as to the attempt of capturing the overall behavior of different rings with only one unique rate, calculated for d _x,y = d _min. To show this, we consider a metalens with perfect spatial positioning but with an additional broadening rate ∼ 2.5 Γ dis ′ ( δ d , d min ) , and we then use the results of Figure 8b to obtain the curves shown in Figure 9. As long as the displacement is small (solid part of the curves), these approximated predictions are in good agreement with the numerical points.

5 Discussion

Complete wavefront shaping requires the simultaneous achievement of high transmittance and full phase control. Usually, metamaterials achieve these requirements by engineering the local properties of the individual scatterers, such as, for example, the shape of nanoresonators. Solid-state, atom-like emitters, however, do not provide the same manufacturing flexibility, and theoretical proposals of atom-based metasurfaces rely on external drives with subwavelength intensity profiles to locally change the emitter properties [28], [29], [30], [37]. Still, the possibility of engineering a complex optical response by solely implanting atomic-scale scatterers in a solid-state environment represents an interesting perspective on device integration and miniaturizability [95], especially when considering the thick substrate that is usually required by standard metasurfaces (typically ∼ 1 mm [34]).

In this work, we showed that stacks of two or more consecutive arrays of solid-state emitters can be engineered to fulfill the necessary requirements of transmittance and phase control, by only choosing proper lattice constants that ensure their correct collective response. Via large-scale numerical simulations [50], we argued that these elements can be combined as the building blocks of a metalens, whose efficiency is robust to losses and other imperfections, due to the collective enhancement of the optical response. This is achieved within a maximum thickness of ∼ 2 λ 0 / 3 , which might be potentially reduced even further, by properly addressing the more complicated regime of evanescent interactions. Notably, the perfect tunability of these building blocks and the possibility of their combination can in principle guarantee arbitrary wavefront shaping, which suggests the extension of this mechanism to more articulated applications, such as phase-only holograms [96].

The core design of our atomic metalens is based on an analytic map between any discretized phase pattern and the corresponding set of lattice constants. Although this scheme is intrinsically scalable, the design is complete only up to three macroscopic free variables, given by the overall phase shift ϕ ₀ of the metalens, the discretization size of the rings ΔR, and the fraction α of “buffer zones”. The scalability of this optimization step is not trivial, as it involves large-scale coupled-dipole simulations. To facilitate it, one possible strategy would consist of investigating whether each ring made of discrete atoms could be modeled by smooth, flat mirrors, with proper transmission and reflection coefficients. This would enable simulation via optical commercial software, with a computational complexity decoupled from the number of dipoles [97], [98], [99], [100], [101]. Alternatively, it would be interesting to explore if a target, collective optical response could be obtained with far fewer emitters, by inverse-designing their positions through proper optimization algorithms [102]. Some preliminary numerical simulations suggest that adjoint methods might be a promising path in this direction [103].

With our scheme, the total efficiency is protected by the collective response, even if the losses of the individual scatterers are non-negligible Γ′ ≫ Γ₀. Similar considerations apply beyond the case of atom-like emitters, to any set of optical scatterers with a well-defined resonant, dipolar response, and a ratio between scattering and total cross section equating Γ 0 / Γ 0 + Γ ′ [36]. This would be the case of plasmonic nanoparticles, for example, which are indeed known to become more resistant to their intrinsic losses when collectively (i.e., nonlocally) responding to light in a 2D, subwavelength array [104]. Our work, based on the idea of combining different arrays together, can then provide additional insights and tools to the context of nonlocal metasurfaces [105].

Finally, it is interesting to mention some specific features of color centers in diamond, whose two-level nature provides nontrivial properties both at the classical and at the quantum level. For example, an atomic metalens based on SiVs would be extremely narrowband and polarization sensitive, finding possible applications in terms of spectral filtering [106], [107], [108], tunability [109], or polarization control [110], [111]. Furthermore, color centers are highly saturable objects, due to their intrinsic nonlinearity, and this behavior would automatically limit the metalens response up to a threshold intensity of light.

At the quantum level, it is known that color centers can be embedded inside a metasurface to enhance some of their functionalities, for example as single-photon sources [112]. It would be interesting to explore if enhanced, collective properties of an ensemble of color centers could be more easily designed by engineering the emitters to act as a nonlocal metasurface. Some evidence exist, for example, that stacks of two atomic arrays can exhibit enhanced nonlinear correlations [86]. More generally, a metasurface based on color centers could provide a possible playground for the emerging contexts of quantum metasurfaces [113] and quantum holography [114], [115].

Methods: We numerically simulate the optical response of the system by solving the coupled-dipole equations of Eqs. (2) and (4), whose computational time scales as ∼N ³, where N is the number of atomic dipoles. The input Gaussian beam must have a waist w ₀ much smaller than the radius R _lens of the atomic metalens, to avoid scattering from the edges or a non-negligible fraction of light passing outside the lens. Due to the paraxial approximation, however, this imposes the constraint λ ₀ ≪ w ₀ ≪ R _lens. Furthermore, to counteract the effects of the broadening Γ′, one must work with small lattice constants down to d _min ≈ 0.03λ₀, thus explaining the necessity of simulating up to N ∼ 5 × 10⁵ atomic dipoles. To accomplish this task, we exploit the fact that the system is symmetric for x ̂ → − x ̂ and y ̂ → − y ̂ , which implies that the each dipole d _j is equal to those of the atoms at the mirrored positions. The actual degrees of freedom are given by the number of atoms satisfying x _j ≥ 0 and y _j ≥ 0, which are roughly N ̃ ∼ N / 4 . The coupled dipole equations can be then simplified by accounting only for these atoms, and then considering as if each of them scattered light from the mirrored positions as well. A supplementary problem is the amount of Random Access Memory (RAM) needed to perform the simulation. We design the code in such a way that the maximum allocation of memory is given by the construction of the N ̃ × N ̃ Green’s function matrix. By defining it as a matrix of Complex{Float32} (64 bit) rather than the custom Complex{Float64} (128 bit), we cut the memory consumption to ∼ 200 – 300 GB of RAM. We checked that we were still working with enough numerical precision, by comparing the simulations of smaller systems, performed with both choices of the variable definition. Finally, the overall computational time was sped up by using the native, multicore implementation of linear algebra in Julia as well as its vectorized treatment of tensor operations [92], while other relevant computations were split over multiple threads. More information is available in the Github repository provided at Ref. [50].

Corresponding author: Francesco Andreoli, ICFO – Institut de Ciències Fotòniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels, Spain, E-mail: francesco.andreoli@alumni.icfo.eu

Funding source: HORIZON EUROPE Marie Sklodowska-Curie Actions

Award Identifier / Grant number: 101068503

Funding source: Generalitat de Catalunya

Award Identifier / Grant number: 2021 SGR 01442

Award Identifier / Grant number: CERCA

Funding source: Horizon 2020 Framework Programme

Award Identifier / Grant number: 101017733

Funding source: European Research Council

Award Identifier / Grant number: 101002107

Funding source: Government of Spain

Award Identifier / Grant number: CEX2019- 000910-S

Award Identifier / Grant number: NextGenerationEU/PRTR PCI2022-132945

Funding source: HORIZON EUROPE European Innovation Council

Award Identifier / Grant number: 101115420

Award Identifier / Grant number: 899275

Funding source: H2020 Marie Skłodowska-Curie Actions

Award Identifier / Grant number: 713729

Research funding: ICFOstepstone – PhD Programme funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 713729. Marie Sklodowska-Curie Actions Postdoctoral Fellowship ATOMAG (grant agreement No. 101068503). European Research Council grant agreement No. 101002107 (NEWSPIN), FET-Open grant agreement No 899275 (DAALI), EIC Pathfinder Grant No. 101115420 (PANDA). Severo Ochoa Grant CEX2019-000910-S [MCIN/AEI/10.13039/501100011033]. QuantERA II project QuSiED, co-funded by the European Union Horizon 2020 research and innovation programme (No. 101017733) and the Government of Spain (European Union NextGenerationEU/PRTR PCI2022-132945 funded by MCIN/AEI/10.13039/501100011033). Generalitat de Catalunya (CERCA program and AGAUR Project No. 2021 SGR 01442), Fundacio Cellex, and Fundacio Mir-Puig.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: Authors state no conflicts of interest.
Data availability: No experimental datasets were generated or analyzed during the current study. All codes used to perform the numerical simulations are available upon request.

Appendix A: Coherent scattering by silicon vacancy centers

In this appendix, we explain more in detail our model of SiV centers in diamond as two-level, dipole emitters. We stress that similar considerations apply to other group IV color centers [60]. The Zero-Phonon-Line (ZPL) of a SiV is centered around 2πc/ω ₀ ≈ 737 nm and is composed of four resonances, associated to the spin–orbit splitting into two ground g ± and two excited states e ± [116]. At cryogenic temperatures ∼ 4 K, these resonances become spectrally resolved [117], and one can target the brightest spectral line (between the lowest ground g − and the lowest excited state e − ) to obtain an effective two-level emitter, with a well-defined dipole moment aligned along the axis between the silicon atom and the carbon vacancy [116]. Although part of the initial population can be in g + rather than g − , in principle this problem can be solved by optical pumping [117], or by further lowering the temperature [118]. At the same time, the inelastic excitation of e + from e − , via phonon coupling, is strongly suppressed already at ∼ 4 K [117]. Still, the target excited state e − can inelasticaly decay into the upper ground state g + or nonradiatively decay out of the ZPL, eventually returning to g − by phononic relaxation [118], [119].

In view of all these considerations, we model the system by considering that the lifetime τ of the excited state e − defines a transform-limited linewidth 1 / τ = Γ 0 + Γ inel ′ + Γ nr ′ composed of several terms [117]. First, we identify the elastic, radiative component Γ 0 = k 0 3 | P 0 | 2 / ( 3 π ϵ ℏ ) , which describes the radiative decay into g − , with P 0 associated to the corresponding dipole matrix element. Second, we include the inelastic Γ inel ′ and nonradiative Γ nr ′ processes by considering that, at ∼ 4 K, they should account for roughly half of the photonic decay [120], leading to Γ inel ′ + Γ nr ′ ≈ Γ 0 . On top of that, the linewidth of the target resonance can undergo homogeneous broadening, which can be caused by nonradiative decay [117], or phonon-induced dephasing [118] and depolarization [117]. In our model, we group these phenomena into an additional rate Γ hom ′ ≲ Γ 0 , whose value qualitatively accounts for the fact that nearly transform-limited linewidths have been observed at cryogenic temperatures (we also notice that encouraging paths have been suggested to extend this property up to much higher temperatures [121]). At the same time, we consider the possibility that local properties (such as strain, or spectral diffusion) randomly shift the resonance frequencies of the individual emitters, thus resulting in an additional inhomogeneous broadening. As detailed in Appendix A1, we model this process with a supplementary rate Γ inhom ′ ≈ 2.75 Γ 0 , where its value is inspired by the experimental results of Ref. [120]. Finally, in Appendix A.2, we show that small disorder in the positions of the solid-state emitters can be similarly modeled by a supplementary inelastic rate, where we take the value of Γ dis ′ ≈ Γ 0 . Given the set of parameters considered, this will roughly correspond to a random displacement within a radius of ∼ 0.1 times the average lattice constants. Overall, this defines the total, additional broadening Γ ′ = Γ inel ′ + Γ nr ′ + Γ hom ′ + Γ inhom ′ + Γ dis ′ ≈ 5.75 Γ 0 .

For each emitter, the quantity Γ 0 / Γ 0 + Γ ′ ≈ 0.15 then estimates, on average, the fraction of resonant scattering events that are coherent and elastic [122].

A.1 Inhomogeneous broadening

To model the presence of inhomogeneous broadening, we assume that each atom of the array has a shifted resonance frequency ω ̃ i , randomly distributed according to the probability distribution P i n h o m ( ω ̃ ) . Here, for simplicity, we focus on a Lorentzian distribution of full-width-half-maximum 2σ _Lorentz, i.e., P i n h o m ( ω ̃ ) = σ L o r e n t z / π σ L o r e n t z 2 + ω ̃ 2 . The system is still described by the coupled-dipole equations Eq. (4), with the difference that each atom now exhibits the shifted polarizability α ′ ( Δ ) = α 0 ( Δ − ω ̃ j ) . In our model, we assume that we can average the atomic response over disorder first, before solving the multiple-scattering problem. We obtain an average atomic polarizability, which reads

(A.1) α 0 = ⟨ α ′ ( ω ̃ ) ⟩ = − 3 π ϵ k 0 3 ∫ Γ 0 P i n h o m ( ω ̃ ) Δ − ω ̃ + i Γ 0 / 2 d ω ̃ = − 3 π ϵ k 0 3 Γ 0 Δ + i Γ 0 + Γ i n h o m ′ / 2 ,

where we defined Γ i n h o m ′ = 2 σ L o r e n t z .

In Figure A.1, we numerically check the soundness of this assumption by evaluating the spectrum of transmission |t(Δ)|² of a 2D square array of size L = 6.4λ ₀ and lattice constant d _x,y = 0.2λ ₀, illuminated by a Gaussian beam of waist w ₀ = L/4. Here, t(Δ) is calculated by projecting the output field of Eq. (2) onto the same mode as the input beam [80], as also detailed in Appendix D. We compare the analytic model of Γ′ = Γ_inhom = 5Γ₀ (red line), with the numerical results obtained by considering a Lorentzian distribution of resonant frequencies with σ _Lorentz = 2.5Γ₀ (blue points), observing a remarkable agreement. As a reference, the green points show the case of resonances ω ̃ i sampled from a Gaussian distribution of standard deviation σ _Gauss = 5Γ₀, whose qualitative similarity leads to the approximate estimation Γ inhom ′ ≈ σ Gauss .

$Figure A.1: Effects of inhomogeneous broadening on a 2D atomic array. Transmission spectrum of a finite 2D square lattice with transverse size L = 6.4λ 0 and lattice constants d x,y = 0.2λ 0, illuminated by a Gaussian beam of waist w 0 = L/4. The blue points (green squares) are calculated by solving the inhomogeneous version of the coupled-dipole equations Eq. (4) with randomly shifted polarizabilities α 0 → α ′ ( ω ̃ i ) ${\alpha }_{0}\to {\alpha }^{\prime }\left({\tilde {\omega }}_{i}\right)$ , and considering atomic resonance frequencies ω ̃ i ${\tilde {\omega }}_{i}$ randomly sampled from a Lorentzian (Gaussian) distribution of half-width-half-maximum σ Lorentz = 2.5Γ0 (standard deviation σ Gauss = 5Γ0). The red line shows the analytic model of a nonradiative decay rate Γ i n h o m ′ = 5 Γ 0 ${{\Gamma}}_{\mathrm{i}\mathrm{n}\mathrm{h}\mathrm{o}\mathrm{m}}^{\prime }=5{{\Gamma}}_{0}$ . The data are averaged over ∼ 100 $\sim 100$ randomly sampled configurations. Similar plots can be derived when calculating the phase of transmission, or the reflection properties.$

Figure A.1:

Effects of inhomogeneous broadening on a 2D atomic array. Transmission spectrum of a finite 2D square lattice with transverse size L = 6.4λ ₀ and lattice constants d _x,y = 0.2λ ₀, illuminated by a Gaussian beam of waist w ₀ = L/4. The blue points (green squares) are calculated by solving the inhomogeneous version of the coupled-dipole equations Eq. (4) with randomly shifted polarizabilities α 0 → α ′ ( ω ̃ i ) , and considering atomic resonance frequencies ω ̃ i randomly sampled from a Lorentzian (Gaussian) distribution of half-width-half-maximum σ _Lorentz = 2.5Γ₀ (standard deviation σ _Gauss = 5Γ₀). The red line shows the analytic model of a nonradiative decay rate Γ i n h o m ′ = 5 Γ 0 . The data are averaged over ∼ 100 randomly sampled configurations. Similar plots can be derived when calculating the phase of transmission, or the reflection properties.

Rather than an exact model of a specific set of experimental data, we aim to capture a reasonable, qualitative description of the effects of inhomogeneous broadening. To this aim, we consider the results of Ref. [120], where they observe a set of 14 SiVs with the same polarization, which exhibit frequencies spanning an interval of Δω ≈ 9.5Γ₀. Assuming that these resonances are uniformly distributed within that bandwidth, we consider a Gaussian distribution with the same standard deviation, which leads to the rough estimation of Γ inhom ′ ≈ Δ ω / 12 ≈ 2.75 Γ 0 . We notice that Lorentzian distributions do not have a well-defined standard deviation, which prompted us to consider the Gaussian case.

A.2 Position disorder

Here, we discuss how small random displacements in the positions of the atomic emitters can affect the optical response. To this aim, we focus on a 2D, square array of constant d _xy, and we uniformly sample the displacement within small spheres of radius δd.

Specifically, we aim to define an effective broadening Γ dis ′ , which should describe the average transmission and reflection of the array. To do so, we consider the reflection r _1L of an idealized infinite array, as expressed in Eq. (6), and we assume that all the losses are contained into Γ dis ′ . At the resonant condition Δ = ω _coop, this allows to define

(A.2) Γ dis ′ Γ coop = 1 | r 1 L ( ω coop ) | − 1 .

In Figure A.2, we thus consider a finite, square array of size L and subwavelength lattice constants d _xy < λ ₀, illuminated by an input Gaussian beam of waist λ ₀ ≪ w ₀ ≪ L, and we numerically compute the reflection r by projecting onto the same Gaussian mode as the input [80]. For each value of d _xy, we randomly displace the positions uniformly within a sphere of radius δd, and define the average reflection ⟨r⟩.

$Figure A.2: Average inelastic scattering due to the disorder in the atomic positions. For each value of the lattice constant d x = d y = d xy of a square, 2D array, we randomly displace the atomic positions within a sphere of radius δd. We then compute the average resonant reflection ⟨r⟩ at Δ = ω coop, upon illumination by an input Gaussian beam of waist w 0 = L/4 ≥ λ 0, where L is the size of the array. Each point is obtained by averaging over 50 random sets of displaced positions. The value of L varied to keep the number of atomic emitters to N = 500. For each configuration, we define the inelastic rate from Eq. (A.2), as Γ dis ′ / Γ coop = | ⟨ r ⟩ | − 1 − 1 ${{\Gamma}}_{\text{dis}}^{\prime }/{{\Gamma}}_{\text{coop}}=\vert \langle r\rangle {\vert }^{-1}-1$ , and we scan increasing radii of disorder δd. The black, dashed line represents the empirical scaling of Eq. (A.3). The insets show the average spectral behavior of reflection (circles) and transmission (triangles) for the cases d xy = 0.25λ 0 and δd ≈ 0.12d xy (up-left inset) or d xy = 0.1λ 0 and δd ≈ 0.34d xy (bottom-right inset). The black, dashed curves represent the predictions of Eq. (A.3), which are in large agreement.$

Figure A.2:

Average inelastic scattering due to the disorder in the atomic positions. For each value of the lattice constant d _x = d _y = d _xy of a square, 2D array, we randomly displace the atomic positions within a sphere of radius δd. We then compute the average resonant reflection ⟨r⟩ at Δ = ω _coop, upon illumination by an input Gaussian beam of waist w ₀ = L/4 ≥ λ ₀, where L is the size of the array. Each point is obtained by averaging over 50 random sets of displaced positions. The value of L varied to keep the number of atomic emitters to N = 500. For each configuration, we define the inelastic rate from Eq. (A.2), as Γ dis ′ / Γ coop = | ⟨ r ⟩ | − 1 − 1 , and we scan increasing radii of disorder δd. The black, dashed line represents the empirical scaling of Eq. (A.3). The insets show the average spectral behavior of reflection (circles) and transmission (triangles) for the cases d _xy = 0.25λ ₀ and δd ≈ 0.12d _xy (up-left inset) or d _xy = 0.1λ ₀ and δd ≈ 0.34d _xy (bottom-right inset). The black, dashed curves represent the predictions of Eq. (A.3), which are in large agreement.

Using Eq. (A.2) as an operative definition, we are able to estimate the value of Γ dis ′ from ⟨r⟩. We reasonably expect that the optical response should be a function of the ratio δd/d _xy. Due to this reason, we plot Γ dis ′ as a function of δd/d _xy in log-log scale, which we numerically fit to obtain the equation

(A.3) Γ dis ′ ≈ π 2 δ d d x y 2 Γ coop = 3 8 δ d λ 0 d x y 2 2 Γ 0 ,

which is represented by a black, dashed line. The numerics confirm our intuition, proving that Eq. (A.3) describes well the reflection properties, at least as long as δd ≲ 0.7d _xy and δd ≪ λ ₀. From a similar analysis, we found that this prediction can capture the spectral behavior of both the reflection ⟨r(Δ)⟩ and the transmission ⟨t(Δ)⟩ coefficients as a function of the detuning (e.g., see the insets of Figure A.2, for the case with either d _xy = 0.025λ ₀ and δd ≈ 0.12d _xy or d _xy = 0.1λ ₀ and δd ≈ 0.34d _xy). The result in Eq. (A.3) should in principle be extendable to the case of rectangular arrays with d _x ≠ d _y, by using the substitution rule d x y 2 → d x d y . We also notice that this result extends the simplified scaling Γ dis ′ ∝ ( δ d / λ 0 ) 2 mentioned in Ref. [17], to the case of arbitrary lattice constants. Finally, we report that the scaling in Eq. (A.3) is confirmed by similar calculations performed on M = 2, 3 square arrays in series.

Appendix B: Evanescent interaction

In this appendix, we further investigate the role of evanescent interactions between 2D atomic arrays, with the goal of justifying the assumption that they are negligible in our regime of interest. We recall that we deal with rectangular, 2D arrays of constants d _x,y ≤ λ ₀, placed at a distance of d _z, and that the dipole matrix elements of the emitters are P 0 = P 0 x ̂ . The evanescent interaction G n m ev results from the evanescent diffraction orders of the field scattered by the atomic layer at z _m, when probed by the atoms at z _n. For subwavelength arrays, its value reads

(A.4) G n m ev = ∑ k x y ( a , b ) ≠ 0 ξ a b 2 k 0 k 0 2 − k x y ( a , b ) ⋅ x ̂ 2 e − | m − n | d z / ξ a b ,

where the diffraction orders are labeled by the integer numbers (a, b), which identify the corresponding wavevector k x y ( a , b ) = 2 π ( a x ̂ / d x + b y ̂ / d y ) and characteristic distance of exponential suppression ξ a b = 1 / | k x y ( a , b ) | 2 − k 0 2 .

The evanescent interaction is stronger for nearest neighbor layers, so we focus on G 12 ev . Moreover, the leading contributions are given by the first two diffraction orders (a, b) = (1, 0) and (a, b) = (0, 1), which are exponentially suppressed by a factor of ∼ 1 / e 2 roughly when d _z = 2max(ξ ₁₀, ξ ₀₁) ≈ max d _x,y/π. The last step is valid for very subwavelength arrays with max d _x,y ≫ λ ₀ and can serve as a simple rule of thumb to roughly identify the regime where G 12 ev ∼ 0 .

Going beyond this rough estimate, in Figure A.3, we numerically calculate the ratio of evanescent to radiative interaction strength | G 12 ev / G 12 rad | , as a function of the lattice constants d _x,y,z. Specifically, on the horizontal axis, we vary the transverse constants along one of the two paths (d _y = d _min) ∪ (d _min ≤ d _x < λ ₀) (Figure A.3a) or (d _x = d _min) ∪ (d _min ≤ d _y < λ ₀) (Figure A.3b). The distance is spanned on the vertical axis within the range d _min ≤ d _z ≤ λ ₀/2, while the white, dotted lines show the specific choice d _z(d _x,y) that we used to define an atomic metalens. The white, solid line shows the rule of thumb d _z ≈ max d _x,y/π. As long as d _y = d _min, d _x ≲ λ ₀/4 or d _x = d _min, d _y ≲ λ ₀/4, the evanescent interaction is completely negligible, being | G 12 ev | / | G 12 rad | < 0.01 (black region). The specific sets of lattice constants used to define the illustrative metalens in the main text (white points) genuinely fall in that regime. By comparing Figure A.3 with Figure 4a, we can infer that almost all phases ϕ _3L can be engineered, except the small range −0.03π ≲ ϕ _3L ≲ 0.06π around ϕ _3L ∼ 0. Two possible ways exist to address this issue. First, one can think of leaving the related ring empty (which would correspond to approximating the phase with ϕ _3L = 0). Otherwise, one can consider larger distances d _z, given that t _3L is invariant (ignoring evanescent interactions) for d _z → d _z + aλ ₀/2 (with a = 1, 2, …), while the effects of evanescent fields rapidly diminish with increasing d _z.

$Figure A.3: Strength of the evanescent interaction between two nearest neighbor layers of atoms. The color legend identifies the relative magnitude | G 12 ev | / | G 12 rad | $\vert {\mathcal{G}}_{12}^{\text{ev}}\vert /\vert {\mathcal{G}}_{12}^{\text{rad}}\vert $ as a function of the lattice constant d x,y and distance d z . We recall that the radiative contribution has a constant magnitude of | G 12 rad | = 1 / 2 $\vert {\mathcal{G}}_{12}^{\text{rad}}\vert =1/2$ , since in Eq. (5) we define the interactions in units of the cooperative rate Γcoop. The red color describes the region where the evanescent field dominates, i.e., | G 12 ev | / | G 12 rad | > 1 $\vert {\mathcal{G}}_{12}^{\text{ev}}\vert /\vert {\mathcal{G}}_{12}^{\text{rad}}\vert { >}1$ . On the contrary, the black area is associated to negligible evanescent interaction | G 12 ev | / | G 12 rad | < 0.01 $\vert {\mathcal{G}}_{12}^{\text{ev}}\vert /\vert {\mathcal{G}}_{12}^{\text{rad}}\vert {< }0.01$ . In the two panels, we explore the two branches of the path chosen for our scheme, reading (d y = d min) ∪ (d min ≤ d x < λ 0) (a) and (d x = d min) ∪ (d min ≤ d y < λ 0) (b), where we recall that d min ≈ 0.03λ 0. The evanescent interaction is calculated from the full equation Eq. (A.4). The white, dotted line represents the possible range of values d z (d x,y ) that guarantee high transmission and full phase control in a three-layer scheme. The white points show the actual values that we used to design the lens of Figure 6, which all fall in a regime where G 12 ev ∼ 0 ${\mathcal{G}}_{12}^{\text{ev}}\sim 0$ . Finally, the white, solid lines show the approximated rule of thumb dz =max dx,y /π.$

Figure A.3:

Strength of the evanescent interaction between two nearest neighbor layers of atoms. The color legend identifies the relative magnitude | G 12 ev | / | G 12 rad | as a function of the lattice constant d _x,y and distance d _z. We recall that the radiative contribution has a constant magnitude of | G 12 rad | = 1 / 2 , since in Eq. (5) we define the interactions in units of the cooperative rate Γ_coop. The red color describes the region where the evanescent field dominates, i.e., | G 12 ev | / | G 12 rad | > 1 . On the contrary, the black area is associated to negligible evanescent interaction | G 12 ev | / | G 12 rad | < 0.01 . In the two panels, we explore the two branches of the path chosen for our scheme, reading (d _y = d _min) ∪ (d _min ≤ d _x < λ ₀) (a) and (d _x = d _min) ∪ (d _min ≤ d _y < λ ₀) (b), where we recall that d _min ≈ 0.03λ ₀. The evanescent interaction is calculated from the full equation Eq. (A.4). The white, dotted line represents the possible range of values d _z(d _x,y) that guarantee high transmission and full phase control in a three-layer scheme. The white points show the actual values that we used to design the lens of Figure 6, which all fall in a regime where G 12 ev ∼ 0 . Finally, the white, solid lines show the approximated rule of thumb d_z =max d_x,y /π.

These conclusions apply for all sets of lattice constants, excluding the limit of d _y → λ ₀. In that specific case, indeed, the diffraction order with (a, b) = (0, 1) would give rise, in Eq. (A.4), to a nominally infinite evanescent contribution arising from the constructive interference between an infinite number of atoms in each 2D layer, associated with an infinite range ξ ₀₁ → ∞ of interaction.

Appendix C: Buffer zones

Here, we describe in detail our definition of the buffer zone between consecutive rings of an atomic metalens. This scheme explicitly takes advantage of the fact that, in our approach, often one of the two lattice constants d _x,y does not change between two consecutive rings. The full algorithm is described below.

Given each ring j, its first 0 ≤ α < 1 fraction is reserved as a buffer zone (green and orange regions of Figure A.4), aimed to connect the array inside the j-th ring with the previous, in a smoother way. Hereafter, we describe how a generic j-th buffer (separating the (j − 1)-th and the j-th ring) is constructed.
First, the system checks if either d x j = d x j − 1 = d min or d y j = d y j − 1 = d min is satisfied. If none of them is fulfilled, then the algorithm ignores that buffer (as in the orange regions of Figure A.4).
Let us assume that one has d y j = d y j − 1 = d min , as in the green regions of Figure A.4. The opposite case is a straightforward extension, which can be described by simply reversing the references to the vertical and horizontal coordinates.
In this regime, the lattices are organized in columns spaced by either d x j − 1 and d x j . The algorithm defines x max = max x j − 1 + ( 3 / 4 ) d x j − 1 , where x _j−1 identify the horizontal positions of the columns of the (j − 1)-th ring. If there are columns of the j-th ring having x _j > x _max, then those columns are ignored in the following steps (as in the black box of Figure A.4).
At this point, the algorithm counts the number of columns in either the j-th or the (j − 1)-th ring, satisfying the condition 0 ≤ x _j,j−1 ≤ x _max. Then, it identifies which of the two rings has less columns. For the sake of simplicity, we will assume it to be the j-th ring, but the algorithm deals with the opposite case in a similar manner. For each column i of this ring, the code searches the horizontally nearest column k among the ones of the (j − 1)-th ring, i.e., the one minimizing the quantity | x j i − x j − 1 k | .
Given this pair of columns, the algorithm connects them by drawing a straight line, and then placing atoms with a vertical spacing d y j = d y j − 1 = d min . For a line to be drawn, the condition y j i > y j − 1 k must be fulfilled. When the number of columns in the two original rings are different, some columns must remain unconnected, as highlighted by the purple box in Figure A.4.
For what concerns the z ̂ position, all the atoms of the j-th buffer are associated to the lattice constant d z j , meaning that the columns are “connected” only in the x ̂ , y ̂ -plane. We tested the idea of fully connecting them in 3D, without noticing significant improvements in the efficiencies.

$Figure A.4: Example of “buffer zones” between two consecutive rings, in the x ̂ , y ̂ $\hat{\mathbf{x}},\hat{\mathbf{y}}$ -plane. The blue points show the atomic positions, while each ring is identified by a red line, as well as an ordinal number, still in red. The first α = 0.2 fraction of each ring is dedicated to the buffer zones, which are represented by either green or orange regions. In particular, the green areas describe the case where one of the two conditions d x j = d x j − 1 = d min ${d}_{x}^{j}={d}_{x}^{j-1}={d}_{\mathrm{min}}$ or d y j = d y j − 1 = d min ${d}_{y}^{j}={d}_{y}^{j-1}={d}_{\mathrm{min}}$ are satisfied, which allows to smoothly connect the neighboring rings. On the contrary, the case where none of these two conditions is fulfilled is shown by the orange zones, which are simply treated as normal parts of the corresponding ring. The black and purple boxes identify two peculiar instances, as described in the main text.$

Figure A.4:

Example of “buffer zones” between two consecutive rings, in the x ̂ , y ̂ -plane. The blue points show the atomic positions, while each ring is identified by a red line, as well as an ordinal number, still in red. The first α = 0.2 fraction of each ring is dedicated to the buffer zones, which are represented by either green or orange regions. In particular, the green areas describe the case where one of the two conditions d x j = d x j − 1 = d min or d y j = d y j − 1 = d min are satisfied, which allows to smoothly connect the neighboring rings. On the contrary, the case where none of these two conditions is fulfilled is shown by the orange zones, which are simply treated as normal parts of the corresponding ring. The black and purple boxes identify two peculiar instances, as described in the main text.

Appendix D: Definition of the efficiency

In our simulations of an atomic metalens, we consider a finite ensemble of N, x ̂ -polarizable atomic emitters, with resonant frequency ω ₀ and embedded in a nonabsorbing, bulk material of index n, so that the resonant wavevector reads k ₀ = 2π/λ ₀ = nω ₀/c. The system is illuminated by a resonant, x ̂ -polarized Gaussian beam of waist w ₀, which reads E _in(R, z) = E _gauss(R, z, w ₀), with

(A.5) E gauss ( R , z , w 0 ) = E 0 w 0 w ( z , w 0 ) exp − | R | 2 w ( z , w 0 ) 2 + i k 0 z + i φ ( | R | , z , w 0 ) x ̂ ,

where w ( z , w 0 ) = w 0 1 + 2 z / k 0 w 0 2 2 is the waist of the beam, while we have φ ( R , z , w 0 ) = − arctan ( 2 z / k 0 w 0 2 ) + k 0 R 2 / [ 2 ρ ( z , w 0 ) ] , with radius of curvature ρ ( z , w 0 ) = z 1 + [ k 0 w 0 2 / ( 2 z ) ] 2 . The total field E _out(R, z) is given by Eqs. (2) and (4) and must be compared to the theoretical output field that one would expect for an ideal lens of focal length f [87]

(A.6) E f ( R , z ) = E gauss ( R , z − z f , w f ) e i k 0 z f M ,

where one has

(A.7) w 0 w f = M = 1 + k 0 w 0 2 2 f 2 ,

and

(A.8) z f = 1 − M − 2 f .

Here, M is the so-called magnification of the lens and ensures energy conservation in the form of P in ∝ ∫ d R | E in | 2 = ∫ d R | E f | 2 = π | E 0 | 2 w 0 2 / 2 . The ideal increase in the beam intensity at the focal point (over the peak input intensity |E ₀|²) is instead given by | E f ( 0 , z f ) | 2 / | E 0 | 2 = M 2 . We can calculate the efficiency η of the atomic metalens by evaluating the overlap between this ideal solution and the total field. In the paraxial limit, this reads η = ⟨ E f | E out ⟩ 2 , where [7], [80]

(A.9) ⟨ E f | E out ⟩ = ∫ R 2 E f * ( R , z ) ⋅ E out ( R , z ) d R ∫ R 2 | E f ( R , z ) | 2 d R = t 0 + 3 i ( k 0 w 0 ) 2 Γ 0 Ω 0 ∑ j = 1 N E f * ( R j , z j ) E 0 * p j P 0 ,

where we have

(A.10) t 0 = ⟨ E f | E in ⟩ = k 0 w 0 w f k 0 w 0 2 + w f 2 / 2 + i z f .

Here, we recall that P 0 = P 0 x ̂ is the dipole matrix element of the emitters, while we have the Rabi frequency Ω 0 = P 0 * E 0 / ℏ . With this definition, the value of η describes the fraction of input power that is transmitted into the desired spatial mode of light. We have made use of the relation ∫ E f * ( R , z ) ⋅ G ̄ ̄ ( R − R j , z − z j ) d R = i E f * ( R j , z j ) / ( 2 k 0 ) , which is true in the far field and as long as the paraxial condition of w _f ≳ λ ₀ (or, more exactly, k ₀ w _f ≫ 1) is satisfied [7], [38], [80]. Similarly, we define the overlap ϵ = ⟨ E in | E out ⟩ 2 between the output and the input mode, where

(A.11) ⟨ E in | E out ⟩ = 1 + 3 i ( k 0 w 0 ) 2 Γ 0 Ω 0 ∑ j = 1 N E in * ( R j , z j ) E 0 * p j P 0 .

Appendix E: Spectral behavior of the metalens

We described a method to engineer an atomic metalens, designed to optimally focus resonant light Δ = ω − ω ₀ = 0. Nonetheless, it is interesting to explore the bandwidth where the efficiency remains high. To address this question, we consider the illustrative example of the main text, corresponding to a metalens with focal length f = 20λ ₀, radius R _lens = 10λ ₀, and constitutive parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≃ −2.06, and α ≈ 0.2, which acts on an input beam of waist w ₀ = 4λ ₀.

Intuitively, we expect the largest bandwidth of nonvanishing optical response to be of the same order of the maximum cooperative decay rate allowed in our system, i.e., Γ coop max = Γ coop ( d x , y = d min ) ≃ 225 Γ 0 . This intuition matches well with what we numerically observe in Figure A.5, where we plot the spectrum of efficiency η (blue), signal-to-background ratio η ̃ (green), and overlap with the input mode ϵ (orange). This is calculated when illuminating our illustrative atomic metalens with a Gaussian beam of waist w ₀ = 4λ ₀, in the lossy regime of Γ′ = 5.75Γ₀. As expected, when | Δ / Γ coop max | ≫ 1 , the metalens shows the features of a transparent system, i.e., E _out ∼E _in, meaning that ϵ ∼ 1, while the efficiency tends to the overlap between the ideal and the input mode, i.e., η ∼ η _in = |⟨E _f|E _in⟩|² ≈ 0.4 (approximately marked with a dark gray region in the plot).

$Figure A.5: Spectral response of the atomic metalens, with focal length f = 20λ 0, radius R lens = 10λ 0, and parameters ΔR ≈ 2λ 0/3, ϕ 0 ≃ −2.06, and α ≈ 0.2. The curves represent the efficiency η (blue), signal-to-background ratio η ̃ $\tilde {\eta }$ (green), and overlap ϵ (orange) with the input beam. The dashed, black, horizontal line shows the value of the overlap between the input and the ideal field η in = |⟨E f |E in⟩|2. The simulation is performed for the lossy case Γ′ = 5.75Γ0. The detuning Δ = ω − ω 0 is expressed either in units of Γ0 (label below) or in units of Γ coop max = Γ coop ( d min ) ≃ 225 Γ 0 ${{\Gamma}}_{\text{coop}}^{\mathrm{max}}={{\Gamma}}_{\text{coop}}\left({d}_{\mathrm{min}}\right)\simeq 225{{\Gamma}}_{0}$ (label above). The dark gray region empirically corresponds to the regime where the atomic emitters become transparent, which roughly reads Δ ≲ − 2 Γ coop max ${\Delta}< sim -2{{\Gamma}}_{\text{coop}}^{\mathrm{max}}$ and Δ ≳ Γ coop max ${\Delta}\gtrsim {{\Gamma}}_{\text{coop}}^{\mathrm{max}}$ . On the contrary, the white region corresponds to the bandwidth | Δ | ≤ ⟨ Γ coop j ⟩ / 2 $\vert {\Delta}\vert \le \langle {{\Gamma}}_{\text{coop}}^{j}\rangle /2$ , where the efficiency remains high η ≳ 0.8. Here, ⟨ Γ coop j ⟩ $\langle {{\Gamma}}_{\text{coop}}^{j}\rangle $ is the average decay rate within the rings, weighted by the fraction of light power illuminating each ring.$

Figure A.5:

Spectral response of the atomic metalens, with focal length f = 20λ ₀, radius R _lens = 10λ ₀, and parameters ΔR ≈ 2λ ₀/3, ϕ ₀ ≃ −2.06, and α ≈ 0.2. The curves represent the efficiency η (blue), signal-to-background ratio η ̃ (green), and overlap ϵ (orange) with the input beam. The dashed, black, horizontal line shows the value of the overlap between the input and the ideal field η _in = |⟨E _f|E _in⟩|². The simulation is performed for the lossy case Γ′ = 5.75Γ₀. The detuning Δ = ω − ω ₀ is expressed either in units of Γ₀ (label below) or in units of Γ coop max = Γ coop ( d min ) ≃ 225 Γ 0 (label above). The dark gray region empirically corresponds to the regime where the atomic emitters become transparent, which roughly reads Δ ≲ − 2 Γ coop max and Δ ≳ Γ coop max . On the contrary, the white region corresponds to the bandwidth | Δ | ≤ ⟨ Γ coop j ⟩ / 2 , where the efficiency remains high η ≳ 0.8. Here, ⟨ Γ coop j ⟩ is the average decay rate within the rings, weighted by the fraction of light power illuminating each ring.

On the contrary, the behavior inside the light-gray area is irregular, but we can identify a bandwidth (white area) of ∼ ± ⟨ Γ c o o p j ⟩ / 2 where the efficiency remains as high as η ≳ 0.8. Here, we defined the average decay rate ⟨ Γ c o o p j ⟩ ≈ 96 Γ 0 by calculating the decay rates Γ c o o p j within each ring that compose the metalens, and then computing the mean value, after weighting each element with the fraction of input light that illuminates the area of the corresponding ring. The value of these weights is illustrated by the colors of the points in Figure 4b. Finally, we stress that the values of Γ coop max and ⟨ Γ c o o p j ⟩ are related to our particular choice of d _min. In general, we can identify a trade-off between the tightness of the bandwidth and the resistance to losses, meaning that some applications that require smaller bandwidths, but can tolerate lower efficiencies, can opt for higher values of d _min.

References

[1] M. Gross and S. Haroche, “Superradiance: an essay on the theory of collective spontaneous emission,” Phys. Rep., vol. 93, no. 5, pp. 301–396, 1982. https://doi.org/10.1016/0370-1573(82)90102-8.Search in Google Scholar

[2] A. S. Sheremet, et al.., “Waveguide quantum electrodynamics: collective radiance and photon-photon correlations,” Rev. Mod. Phys., vol. 95, no. 1, p. 015002, 2023. https://doi.org/10.1103/revmodphys.95.015002.Search in Google Scholar

[3] S. D. Jenkins and J. Ruostekoski, “Controlled manipulation of light by cooperative response of atoms in an optical lattice,” Phys. Rev. A, vol. 86, no. 3, p. 031602, 2012. https://doi.org/10.1103/physreva.86.031602.Search in Google Scholar

[4] G. Facchinetti, S. D. Jenkins, and J. Ruostekoski, “Storing light with subradiant correlations in arrays of atoms,” Phys. Rev. Lett., vol. 117, no. 24, p. 243601, 2016. https://doi.org/10.1103/physrevlett.117.243601.Search in Google Scholar PubMed

[5] A. Asenjo-Garcia, et al.., “Exponential improvement in photon storage fidelities using subradiance and “selective radiance” in atomic arrays,” Phys. Rev. X, vol. 7, no. 3, p. 31024, 2017. https://doi.org/10.1103/physrevx.7.031024.Search in Google Scholar

[6] J. Perczel, et al.., “Photonic band structure of two-dimensional atomic lattices,” Phys. Rev. A, vol. 96, no. 6, p. 063801, 2017. https://doi.org/10.1103/physreva.96.063801.Search in Google Scholar

[7] M. T. Manzoni, et al.., “Optimization of photon storage fidelity in ordered atomic arrays,” New J. Phys., vol. 20, no. 8, p. 83048, 2018. https://doi.org/10.1088/1367-2630/aadb74.Search in Google Scholar PubMed PubMed Central

[8] P. O. Guimond, et al.., “Subradiant bell states in distant atomic arrays,” Phys. Rev. Lett., vol. 122, no. 9, p. 093601, 2019. https://doi.org/10.1103/physrevlett.122.093601.Search in Google Scholar

[9] R. J. Bettles, et al.., “Quantum and nonlinear effects in light transmitted through planar atomic arrays,” Commun. Phys., vol. 3, no. 141, pp. 1–9, 2020. https://doi.org/10.1038/s42005-020-00404-3.Search in Google Scholar

[10] K. Brechtelsbauer and D. Malz, “Quantum simulation with fully coherent dipole-dipole interactions mediated by three-dimensional subwavelength atomic arrays,” Phys. Rev. A, vol. 104, no. 1, p. 013701, 2021. https://doi.org/10.1103/physreva.104.013701.Search in Google Scholar

[11] Z. Y. Wei, et al.., “Generation of photonic matrix product states with Rydberg atomic arrays,” Phys. Rev. Res., vol. 3, no. 2, p. 023021, 2021. https://doi.org/10.1103/physrevresearch.3.023021.Search in Google Scholar

[12] C. C. Rusconi, T. Shi, and J. I. Cirac, “Exploiting the photonic nonlinearity of free-space subwavelength arrays of atoms,” Phys. Rev. A, vol. 104, no. 3, p. 033718, 2021. https://doi.org/10.1103/physreva.104.033718.Search in Google Scholar

[13] T. L. Patti, et al.., “Controlling interactions between quantum emitters using atom arrays,” Phys. Rev. Lett., vol. 126, no. 22, p. 223602, 2021. https://doi.org/10.1103/physrevlett.126.223602.Search in Google Scholar

[14] E. Sierra, S. J. Masson, and A. Asenjo-Garcia, “Dicke superradiance in ordered lattices: dimensionality matters,” Phys. Rev. Res., vol. 4, no. 2, p. 023207, 2022. https://doi.org/10.1103/physrevresearch.4.023207.Search in Google Scholar

[15] S. J. Masson and A. Asenjo-Garcia, “Universality of Dicke superradiance in arrays of quantum emitters,” Nat. Commun., vol. 13, no. 2285, pp. 1–7, 2022. https://doi.org/10.1038/s41467-022-29805-4.Search in Google Scholar PubMed PubMed Central

[16] F. Andreoli, et al.., “The maximum refractive index of an atomic crystal – from quantum optics to quantum chemistry,” arXiv:2303.10998, 2023.Search in Google Scholar

[17] Y. Solomons, R. Ben-Maimon, and E. Shahmoon, “Universal approach for quantum interfaces with atomic arrays,” PRX Quantum, vol. 5, no. 2, p. 020329, 2024.10.1103/PRXQuantum.5.020329Search in Google Scholar

[18] M. Moreno-Cardoner, D. Goncalves, and D. Chang, “Quantum nonlinear optics based on two-dimensional rydberg atom arrays,” Phys. Rev. Lett., vol. 127, no. 26, p. 263602, 2021. https://doi.org/10.1103/physrevlett.127.263602.Search in Google Scholar PubMed

[19] O. Rubies-Bigorda, et al.., “Photon control and coherent interactions via lattice dark states in atomic arrays,” Phys. Rev. Res., vol. 4, no. 1, p. 013110, 2022. https://doi.org/10.1103/physrevresearch.4.013110.Search in Google Scholar

[20] K. E. Ballantine and J. Ruostekoski, “Subradiance-protected excitation spreading in the generation of collimated photon emission from an atomic array,” Phys. Rev. Res., vol. 2, no. 2, p. 023086, 2020. https://doi.org/10.1103/physrevresearch.2.023086.Search in Google Scholar

[21] K. E. Ballantine and J. Ruostekoski, “Quantum single-photon control, storage, and entanglement generation with planar atomic arrays,” PRX Quantum, vol. 2, no. 4, p. 040362, 2021. https://doi.org/10.1103/prxquantum.2.040362.Search in Google Scholar

[22] J. Ruostekoski, “Cooperative quantum-optical planar arrays of atoms,” Phys. Rev. A, vol. 108, no. 3, p. 030101, 2023. https://doi.org/10.1103/physreva.108.030101.Search in Google Scholar

[23] R. J. Bettles, S. A. Gardiner, and C. S. Adams, “Enhanced optical cross section via collective coupling of atomic dipoles in a 2D array,” Phys. Rev. Lett., vol. 116, no. 10, p. 103602, 2016. https://doi.org/10.1103/physrevlett.116.103602.Search in Google Scholar

[24] E. Shahmoon, et al.., “Cooperative resonances in light scattering from two-dimensional atomic arrays,” Phys. Rev. Lett., vol. 118, no. 11, p. 113601, 2017. https://doi.org/10.1103/physrevlett.118.113601.Search in Google Scholar

[25] J. Rui, et al.., “A subradiant optical mirror formed by a single structured atomic layer,” Nature, vol. 583, no. 7816, pp. 369–374, 2020. https://doi.org/10.1038/s41586-020-2463-x.Search in Google Scholar PubMed

[26] N. Nefedkin, M. Cotrufo, and A. Alù, “Nonreciprocal total cross section of quantum metasurfaces,” Nanophotonics, vol. 12, no. 3, pp. 589–606, 2023. https://doi.org/10.1515/nanoph-2022-0596.Search in Google Scholar PubMed PubMed Central

[27] R. Alaee, et al.., “Quantum metamaterials with magnetic response at optical frequencies,” Phys. Rev. Lett., vol. 125, no. 6, p. 063601, 2020. https://doi.org/10.1103/physrevlett.125.063601.Search in Google Scholar PubMed

[28] K. E. Ballantine and J. Ruostekoski, “Optical magnetism and huygens’ surfaces in arrays of atoms induced by cooperative responses,” Phys. Rev. Lett., vol. 125, no. 14, p. 143604, 2020. https://doi.org/10.1103/physrevlett.125.143604.Search in Google Scholar PubMed

[29] K. E. Ballantine, D. Wilkowski, and J. Ruostekoski, “Optical magnetism and wavefront control by arrays of strontium atoms,” Phys. Rev. Res., vol. 4, no. 3, p. 033242, 2022. https://doi.org/10.1103/physrevresearch.4.033242.Search in Google Scholar

[30] K. E. Ballantine and J. Ruostekoski, “Cooperative optical wavefront engineering with atomic arrays,” Nanophotonics, vol. 10, no. 7, pp. 1901–1909, 2021. https://doi.org/10.1515/nanoph-2021-0059.Search in Google Scholar

[31] B. X. Wang, et al.., “Design of metasurface polarizers based on two-dimensional cold atomic arrays,” Opt. Express, vol. 25, no. 16, pp. 18760–18773, 2017. https://doi.org/10.1364/oe.25.018760.Search in Google Scholar PubMed

[32] N. S. Baßler, et al.., “Linear optical elements based on cooperative subwavelength emitter arrays,” Opt. Express, vol. 31, no. 4, pp. 6003–6026, 2023. https://doi.org/10.1364/oe.476830.Search in Google Scholar

[33] N. S. Baßler, et al.., “Metasurface-based hybrid optical cavities for chiral sensing,” Phys. Rev. Lett., vol. 132, no. 4, p. 043602, 2024. https://doi.org/10.1103/physrevlett.132.043602.Search in Google Scholar

[34] J. Engelberg and U. Levy, “The advantages of metalenses over diffractive lenses,” Nat. Commun., vol. 11, no. 1991, 2020. https://doi.org/10.1038/s41467-020-15972-9.Search in Google Scholar PubMed PubMed Central

[35] W. T. Chen, A. Y. Zhu, and F. Capasso, “Flat optics with dispersion-engineered metasurfaces,” Nat. Rev. Mater., vol. 5, no. 8, pp. 604–620, 2020. https://doi.org/10.1038/s41578-020-0203-3.Search in Google Scholar

[36] Z. Li, et al.., “Atomic optical antennas in solids,” Nat. Photonics, vol. 18, no. 10, pp. 1113–1120, 2024. https://doi.org/10.1038/s41566-024-01456-5.Search in Google Scholar

[37] M. Zhou, et al.., “Optical metasurface based on the resonant scattering in electronic transitions,” ACS Photonics, vol. 4, no. 5, pp. 1279–1285, 2017. https://doi.org/10.1021/acsphotonics.7b00219.Search in Google Scholar

[38] L. Chomaz, et al.., “Absorption imaging of a quasi-two-dimensional gas: a multiple scattering analysis,” New J. Phys., vol. 14, no. 5, p. 055001, 2012. https://doi.org/10.1088/1367-2630/14/5/055001.Search in Google Scholar

[39] J. Javanainen, et al.., “Shifts of a resonance line in a dense atomic sample,” Phys. Rev. Lett., vol. 112, no. 11, 2014. https://doi.org/10.1103/physrevlett.112.113603.Search in Google Scholar

[40] J. Javanainen and J. Ruostekoski, “Light propagation beyond the mean-field theory of standard optics,” Opt. Express, vol. 24, no. 2, p. 993, 2016. https://doi.org/10.1364/oe.24.000993.Search in Google Scholar PubMed

[41] B. Zhu, et al.., “Light scattering from dense cold atomic media,” Phys. Rev. A, vol. 94, no. 2, p. 023612, 2016. https://doi.org/10.1103/physreva.94.023612.Search in Google Scholar

[42] N. J. Schilder, et al.., “Polaritonic modes in a dense cloud of cold atoms,” Phys. Rev. A, vol. 93, no. 6, p. 063835, 2016. https://doi.org/10.1103/physreva.93.063835.Search in Google Scholar

[43] N. J. Schilder, et al.., “Homogenization of an ensemble of interacting resonant scatterers,” Phys. Rev. A, vol. 96, no. 1, p. 013825, 2017. https://doi.org/10.1103/physreva.96.013825.Search in Google Scholar

[44] N. Schilder, et al.., “Near-resonant light scattering by a subwavelength ensemble of identical atoms,” Phys. Rev. Lett., vol. 124, no. 7, p. 073403, 2020. https://doi.org/10.1103/physrevlett.124.073403.Search in Google Scholar

[45] S. Jennewein, et al.., “Propagation of light through small clouds of cold interacting atoms,” Phys. Rev. A, vol. 94, no. 5, 2016, https://doi.org/10.1103/physreva.94.053828.Search in Google Scholar

[46] S. Jennewein, et al.., “Coherent scattering of near-resonant light by a dense, microscopic cloud of cold two-level atoms: experiment versus theory,” Phys. Rev. A, vol. 97, no. 5, 2018. https://doi.org/10.1103/physreva.97.053816.Search in Google Scholar

[47] L. Corman, et al.., “Transmission of near-resonant light through a dense slab of cold atoms,” Phys. Rev. A, vol. 96, no. 5, p. 53629, 2017. https://doi.org/10.1103/physreva.96.053629.Search in Google Scholar

[48] S. D. Jenkins, et al.., “Collective resonance fluorescence in small and dense atom clouds: comparison between theory and experiment,” Phys. Rev. A, vol. 94, no. 2, p. 023842, 2016. https://doi.org/10.1103/physreva.94.023842.Search in Google Scholar

[49] H. Dobbertin, R. Löw, and S. Scheel, “Collective dipole-dipole interactions in planar nanocavities,” Phys. Rev. A, vol. 102, no. 3, p. 031701, 2020. https://doi.org/10.1103/physreva.102.031701.Search in Google Scholar

[50] The codes are available at the Github repositories: https://github.com/frandreoli/atoms_optical_response (optical simulations) and https://github.com/frandreoli/optimization_atoms_metalens (metalens efficiency optimization). The optical simulations allow to reconstruct the field at any point in space and calculate reflection, transmission and optical efficiencies, for a set of atomic emitters at arbitrary positions, including those of an atomic metalens. The optimization toolbox can implement many different algorithms, with the particle swarm empirically outperforming the others.Search in Google Scholar

[51] D. Fattal, et al.., “Flat dielectric grating reflectors with focusing abilities,” Nat. Photonics, vol. 4, no. 7, pp. 466–470, 2010. https://doi.org/10.1038/nphoton.2010.116.Search in Google Scholar

[52] A. B. Klemm, et al.., “Experimental high numerical aperture focusing with high contrast gratings,” Opt. Lett., vol. 38, no. 17, p. 3410, 2013. https://doi.org/10.1364/ol.38.003410.Search in Google Scholar PubMed

[53] M. Khorasaninejad, et al.., “Metalenses at visible wavelengths: diffraction-limited focusing and subwavelength resolution imaging,” Science, vol. 352, no. 6290, pp. 1190–1194, 2016. https://doi.org/10.1126/science.aaf6644.Search in Google Scholar PubMed

[54] Z. Zhou, et al.., “Efficient silicon metasurfaces for visible light,” ACS Photonics, vol. 4, no. 3, pp. 544–551, 2017. https://doi.org/10.1021/acsphotonics.6b00740.Search in Google Scholar

[55] H. Liang, et al.., “Ultrahigh numerical aperture metalens at visible wavelengths,” Nano Lett., vol. 18, no. 7, pp. 4460–4466, 2018. https://doi.org/10.1021/acs.nanolett.8b01570.Search in Google Scholar PubMed

[56] A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Planar photonics with metasurfaces,” Science, vol. 339, no. 6125, pp. 12320091–12320096, 2013. https://doi.org/10.1126/science.1232009.Search in Google Scholar PubMed

[57] N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater., vol. 13, no. 2, pp. 139–150, 2014. https://doi.org/10.1038/nmat3839.Search in Google Scholar PubMed

[58] W. T. Chen and F. Capasso, “Will flat optics appear in everyday life anytime soon?,” Appl. Phys. Lett., vol. 118, no. 10, p. 100503, 2021. https://doi.org/10.1063/5.0039885.Search in Google Scholar

[59] S. Shrestha, et al.., “Broadband achromatic dielectric metalenses,” Light Sci. Appl., vol. 7, no. 85, 2018. https://doi.org/10.1038/s41377-018-0078-x.Search in Google Scholar PubMed PubMed Central

[60] C. Bradac, et al.., “Quantum nanophotonics with group IV defects in diamond,” Nat. Commun., vol. 10, no. 5625, pp. 1–13, 2019. https://doi.org/10.1038/s41467-019-13332-w.Search in Google Scholar PubMed PubMed Central

[61] V. M. Acosta, et al.., “Diamonds with a high density of nitrogen-vacancy centers for magnetometry applications,” Phys. Rev. B, vol. 80, no. 11, 2009. https://doi.org/10.1103/physrevb.80.115202.Search in Google Scholar

[62] C. Hepp, et al.., “Electronic structure of the silicon vacancy color center in diamond,” Phys. Rev. Lett., vol. 112, no. 3, p. 036405, 2014. https://doi.org/10.1103/physrevlett.112.036405.Search in Google Scholar

[63] T. Müller, et al.., “Optical signatures of silicon-vacancy spins in diamond,” Nat. Commun., vol. 5, no. 3328, pp. 1–7, 2014. https://doi.org/10.1038/ncomms4328.Search in Google Scholar PubMed

[64] T. Iwasaki, et al.., “Tin-vacancy quantum emitters in diamond,” Phys. Rev. Lett., vol. 119, no. 25, p. 253601, 2017. https://doi.org/10.1103/physrevlett.119.253601.Search in Google Scholar PubMed

[65] M. E. Trusheim, et al.., “Lead-related quantum emitters in diamond,” Phys. Rev. B, vol. 99, no. 7, p. 075430, 2019. https://doi.org/10.1103/physrevb.99.075430.Search in Google Scholar

[66] J. M. Smith, et al.., “Colour centre generation in diamond for quantum technologies,” Nanophotonics, vol. 8, no. 11, pp. 1889–1906, 2019. https://doi.org/10.1515/nanoph-2019-0196.Search in Google Scholar

[67] K. Ohno, et al.., “Three-dimensional localization of spins in diamond using 12C implantation,” Appl. Phys. Lett., vol. 105, no. 5, p. 52406, 2014. https://doi.org/10.1063/1.4890613.Search in Google Scholar

[68] T. Y. Hwang, et al.., “Sub-10 nm precision engineering of solid-state defects via nanoscale aperture array mask,” Nano Lett., vol. 22, no. 4, pp. 1672–1679, 2022. https://doi.org/10.1021/acs.nanolett.1c04699.Search in Google Scholar PubMed

[69] J. Michl, et al.., “Perfect alignment and preferential orientation of nitrogen-vacancy centers during chemical vapor deposition diamond growth on (111) surfaces,” Appl. Phys. Lett., vol. 104, no. 10, p. 102407, 2014. https://doi.org/10.1063/1.4868128.Search in Google Scholar

[70] M. Lesik, et al.., “Perfect preferential orientation of nitrogen-vacancy defects in a synthetic diamond sample,” Appl. Phys. Lett., vol. 104, no. 11, p. 113107, 2014. https://doi.org/10.1063/1.4869103.Search in Google Scholar

[71] T. Fukui, et al.., “Perfect selective alignment of nitrogen-vacancy centers in diamond,” Appl. Phys. Express, vol. 7, no. 5, p. 055201, 2014. https://doi.org/10.7567/apex.7.055201.Search in Google Scholar

[72] H. Ozawa, et al.., “Formation of perfectly aligned nitrogen-vacancy-center ensembles in chemical-vapor-deposition-grown diamond (111),” Appl. Phys. Express, vol. 10, no. 4, p. 045501, 2017. https://doi.org/10.7567/apex.10.045501.Search in Google Scholar

[73] J. L. Pacheco, et al.., “Ion implantation for deterministic single atom devices,” Rev. Sci. Instrum., vol. 88, no. 12, 2017. https://doi.org/10.1063/1.5001520.Search in Google Scholar PubMed

[74] Y.-C. Chen, et al.., “Laser writing of individual nitrogen-vacancy defects in diamond with near-unity yield,” Optica, vol. 6, no. 5, pp. 662–667, 2019. https://doi.org/10.1364/optica.6.000662.Search in Google Scholar

[75] Y. Zhou, et al.., “Direct writing of single germanium vacancy center arrays in diamond,” New J. Phys., vol. 20, no. 12, p. 125004, 2018. https://doi.org/10.1088/1367-2630/aaf2ac.Search in Google Scholar

[76] D. Scarabelli, et al.., “Nanoscale engineering of closely-spaced electronic spins in diamond,” Nano Lett., vol. 16, no. 8, pp. 4982–4990, 2016. https://doi.org/10.1021/acs.nanolett.6b01692.Search in Google Scholar PubMed

[77] C. J. Stephen, et al.., “Deep three-dimensional solid-state qubit arrays with long-lived spin coherence,” Phys. Rev. Appl., vol. 12, no. 6, p. 064005, 2019. https://doi.org/10.1103/physrevapplied.12.064005.Search in Google Scholar

[78] X. Guo, et al.., “Tunable and transferable diamond membranes for integrated quantum technologies,” Nano Lett., vol. 21, no. 24, pp. 10392–10399, 2021. https://doi.org/10.1021/acs.nanolett.1c03703.Search in Google Scholar PubMed PubMed Central

[79] J. Ren, et al.., “Three-dimensional superlattice engineering with block copolymer epitaxy,” Sci. Adv., vol. 6, no. 24, 2020. https://doi.org/10.1126/sciadv.aaz0002.Search in Google Scholar PubMed PubMed Central

[80] F. Andreoli, et al.., “Maximum refractive index of an atomic medium,” Phys. Rev. X, vol. 11, no. 1, p. 011026, 2021. https://doi.org/10.1103/physrevx.11.011026.Search in Google Scholar

[81] L. Novotny and B. Hecht, Principles of Nano-Optics, Cambridge, Cambridge University Press, 2009.Search in Google Scholar

[82] M. Antezza and Y. Castin, “Spectrum of light in a quantum fluctuating periodic structure,” Phys. Rev. Lett., vol. 103, no. 12, p. 123903, 2009. https://doi.org/10.1103/physrevlett.103.123903.Search in Google Scholar PubMed

[83] M. Antezza and Y. Castin, “Fano-Hopfield model and photonic band gaps for an arbitrary atomic lattice,” Phys. Rev. A, vol. 80, no. 1, p. 013816, 2009. https://doi.org/10.1103/physreva.80.013816.Search in Google Scholar

[84] C.-R. Mann, et al.., “Selective radiance in super-wavelength atomic arrays,” arXiv:2402.06439, 2024.Search in Google Scholar

[85] R. Ben-Maimon, et al.., “Quantum interfaces with multilayered superwavelength atomic arrays,” arXiv:2402.06839, 2024.Search in Google Scholar

[86] S. P. Pedersen, L. Zhang, and T. Pohl, “Quantum nonlinear metasurfaces from dual arrays of ultracold atoms,” Phys. Rev. Res., vol. 5, no. 1, p. L012047, 2023. https://doi.org/10.1103/physrevresearch.5.l012047.Search in Google Scholar

[87] B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics, New York, John Wiley & Sons, Inc., 1991.10.1002/0471213748Search in Google Scholar

[88] H. Zheng and H. U. Baranger, “Persistent quantum beats and long-distance entanglement from waveguide-mediated interactions,” Phys. Rev. Lett., vol. 110, no. 11, p. 113601, 2013. https://doi.org/10.1103/physrevlett.110.113601.Search in Google Scholar

[89] H. van de Stadt and J. M. Muller, “Multimirror Fabry–Perot interferometers,” JOSA A, vol. 2, no. 8, p. 1363, 1985. https://doi.org/10.1364/josaa.2.001363.Search in Google Scholar

[90] I. H. Deutsch, et al.., “Photonic band gaps in optical lattices,” Phys. Rev. A, vol. 52, no. 2, pp. 1394–1410, 1995. https://doi.org/10.1103/physreva.52.1394.Search in Google Scholar PubMed

[91] R. Menon and B. Sensale-Rodriguez, “Inconsistencies of metalens performance and comparison with conventional diffractive optics,” Nat. Photonics, vol. 17, no. 11, pp. 923–924, 2023. https://doi.org/10.1038/s41566-023-01306-w.Search in Google Scholar

[92] J. Bezanson, et al.., “Julia: a fresh approach to numerical computing,” SIAM Rev., vol. 59, no. 1, pp. 65–98, 2017. https://doi.org/10.1137/141000671.Search in Google Scholar

[93] R. E. Evans, et al.., “Narrow-linewidth homogeneous optical emitters in diamond nanostructures via silicon ion implantation,” Phys. Rev. Appl., vol. 5, no. 4, p. 044010, 2016. https://doi.org/10.1103/physrevapplied.5.044010.Search in Google Scholar

[94] T. Schröder, et al.., “Scalable focused ion beam creation of nearly lifetime-limited single quantum emitters in diamond nanostructures,” Nat. Commun., vol. 8, no. 15376, pp. 1–7, 2017. https://doi.org/10.1038/ncomms15376.Search in Google Scholar PubMed PubMed Central

[95] M. Pan, et al.., “Dielectric metalens for miniaturized imaging systems: progress and challenges,” Light Sci. Appl., vol. 11, no. 195, pp. 1–32, 2022. https://doi.org/10.1038/s41377-022-00885-7.Search in Google Scholar PubMed PubMed Central

[96] L. Huang, S. Zhang, and T. Zentgraf, “Metasurface holography: from fundamentals to applications,” Nanophotonics, vol. 7, no. 6, pp. 1169–1190, 2018. https://doi.org/10.1515/nanoph-2017-0118.Search in Google Scholar

[97] C. M. Lalau-Keraly, et al.., “Adjoint shape optimization applied to electromagnetic design,” Opt. Express, vol. 21, no. 18, pp. 21693–21701, 2013. https://doi.org/10.1364/oe.21.021693.Search in Google Scholar PubMed

[98] R. Paniagua-Domínguez, et al.., “A metalens with a near-unity numerical aperture,” Nano Lett., vol. 18, no. 3, pp. 2124–2132, 2018. https://doi.org/10.1021/acs.nanolett.8b00368.Search in Google Scholar PubMed

[99] R. E. Christiansen, et al.., “Inverse design of nanoparticles for enhanced Raman scattering,” Opt. Express, vol. 28, no. 4, pp. 4444–4462, 2020. https://doi.org/10.1364/oe.28.004444.Search in Google Scholar PubMed

[100] H. Chung and O. D. Miller, “High-NA achromatic metalenses by inverse design,” Opt. Express, vol. 28, no. 5, pp. 6945–6965, 2020. https://doi.org/10.1364/oe.385440.Search in Google Scholar

[101] O. A. Abdelraouf, et al.., “Recent advances in tunable metasurfaces: materials, design, and applications,” ACS Nano, vol. 16, no. 9, pp. 13339–13369, 2022. https://doi.org/10.1021/acsnano.2c04628.Search in Google Scholar PubMed

[102] I. Volkov, et al.., “Non-radiative configurations of a few quantum emitters ensembles: evolutionary optimization approach,” Appl. Phys. Lett., vol. 124, no. 8, 2024. https://doi.org/10.1063/5.0189405.Search in Google Scholar

[103] O. D. Miller, “Photonic design: from fundamental solar cell physics to computational inverse design,” arXiv:1308.0212, 2013.Search in Google Scholar

[104] M. S. Bin-Alam, et al.., “Ultra-high-Q resonances in plasmonic metasurfaces,” Nat. Commun., vol. 12, no. 974, pp. 1–8, 2021. https://doi.org/10.1038/s41467-021-21196-2.Search in Google Scholar PubMed PubMed Central

[105] K. Shastri and F. Monticone, “Nonlocal flat optics,” Nat. Photonics, vol. 17, no. 1, pp. 36–47, 2022. https://doi.org/10.1038/s41566-022-01098-5.Search in Google Scholar

[106] C. Chen, et al.., “Spectral tomographic imaging with aplanatic metalens,” Light Sci. Appl., vol. 8, no. 1, pp. 1–8, 2019. https://doi.org/10.1038/s41377-019-0208-0.Search in Google Scholar PubMed PubMed Central

[107] A. Arbabi, et al.., “Subwavelength-thick lenses with high numerical apertures and large efficiency based on high-contrast transmitarrays,” Nat. Commun., vol. 6, no. 1, pp. 1–6, 2015. https://doi.org/10.1038/ncomms8069.Search in Google Scholar PubMed

[108] A. McClung, et al.., “Snapshot spectral imaging with parallel metasystems,” Sci. Adv., vol. 6, no. 38, pp. 7646–7664, 2020. https://doi.org/10.1126/sciadv.abc7646.Search in Google Scholar PubMed PubMed Central

[109] J. van de Groep, et al.., “Exciton resonance tuning of an atomically thin lens,” Nat. Photonics, vol. 14, no. 7, pp. 426–430, 2020. https://doi.org/10.1038/s41566-020-0624-y.Search in Google Scholar

[110] K. Ou, et al.., “Mid-infrared polarization-controlled broadband achromatic metadevice,” Sci. Adv., vol. 6, no. 37, pp. 711–722, 2020. https://doi.org/10.1126/sciadv.abc0711.Search in Google Scholar PubMed PubMed Central

[111] S. Gao, et al.., “Twofold polarization-selective all-dielectric trifoci metalens for linearly polarized visible light,” Adv. Opt. Mater., vol. 7, no. 21, p. 1900883, 2019. https://doi.org/10.1002/adom.201900883.Search in Google Scholar

[112] J. Ma, et al.., “Engineering quantum light sources with flat optics,” Adv. Mater., p. 2313589, 2024, https://doi.org/10.1002/adma.202313589.Search in Google Scholar PubMed

[113] A. S. Solntsev, G. S. Agarwal, and Y. Y. Kivshar, “Metasurfaces for quantum photonics,” Nat. Photonics, vol. 15, no. 5, pp. 327–336, 2021. https://doi.org/10.1038/s41566-021-00793-z.Search in Google Scholar

[114] J.-Z. Yang, et al.., “Quantum metasurface holography,” Photonics Res., vol. 10, no. 11, pp. 2607–2613, 2022. https://doi.org/10.1364/prj.470537.Search in Google Scholar

[115] Q. Y. Wu, et al.., “Quantum process tomography on holographic metasurfaces,” Npj Quantum Inf., vol. 8, no. 46, pp. 1–6, 2022. https://doi.org/10.1038/s41534-022-00561-z.Search in Google Scholar

[116] L. J. Rogers, et al.., “Electronic structure of the negatively charged silicon-vacancy center in diamond,” Phys. Rev. B, vol. 89, no. 23, p. 235101, 2014. https://doi.org/10.1103/physrevb.89.235101.Search in Google Scholar

[117] K. D. Jahnke, et al.., “Electron-phonon processes of the silicon-vacancy centre in diamond,” New J. Phys., vol. 17, no. 4, p. 043011, 2015. https://doi.org/10.1088/1367-2630/17/4/043011.Search in Google Scholar

[118] D. D. Sukachev, et al.., “Silicon-vacancy spin qubit in diamond: a quantum memory exceeding 10 ms with single-shot state readout,” Phys. Rev. Lett., vol. 119, no. 22, p. 223602, 2017. https://doi.org/10.1103/physrevlett.119.223602.Search in Google Scholar PubMed

[119] A. Sipahigil, et al.., “Indistinguishable photons from separated silicon-vacancy centers in diamond,” Phys. Rev. Lett., vol. 113, no. 11, p. 113602, 2014. https://doi.org/10.1103/physrevlett.113.113602.Search in Google Scholar

[120] L. J. Rogers, et al.., “Multiple intrinsically identical single-photon emitters in the solid state,” Nat. Commun., vol. 5, 2014, https://doi.org/10.1038/ncomms5739.Search in Google Scholar PubMed

[121] P. Wang, et al.., “Transform-limited photon emission from a lead-vacancy center in diamond above 10 K,” Phys. Rev. Lett., vol. 132, no. 7, p. 073601, 2024. https://doi.org/10.1103/physrevlett.132.073601.Search in Google Scholar PubMed

[122] L. Novotny and N. Van Hulst, “Antennas for light,” Nat. Photonics, vol. 5, no. 2, pp. 83–90, 2011. https://doi.org/10.1038/nphoton.2010.237.Search in Google Scholar

Received: 2024-11-04

Accepted: 2024-12-20

Published Online: 2025-01-31

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/nanoph-2024-0603

Keywords for this article

metalens; metasurface; quantum optics; atomic emitters; color centers; atomic arrays

Creative Commons

BY 4.0