Entanglement is an important feature of quantum mechanics that underpins various applications of quantum technologies^[1]. In particular, photon pairs that are entangled in the high-dimensional spatial space represent an essential resource in a broad range of quantum applications, including imaging^[2,3], communications^[4], and computations^[5]. The performance of such applications is related to the degree of entanglement. The most common way to generate spatially entangled photon pairs is based on the spontaneous parametric downconversion (SPDC) in quadratic nonlinear crystals, where a pump photon spontaneously splits into two lower-energy photons in two different directions^[6]. The properties of generated states, including the emission pattern and spatial entanglement, have been extensively investigated^[7–9]. For nonlinear crystals with a typical thickness on the scale of millimeters to centimeters, the stringent phase-matching condition limits the emission directions of the photon pairs to a certain predefined angle range, making it difficult to flexibly tune the spatial pattern and entanglement of the photon pairs while maintaining generation efficiency. Although control of spatial correlations of the photon pairs was recently reported by engineering the pump beam profile^[10], the specific tuning range of the spatial entanglement was unknown. It was shown that strong multimode entanglement can be achieved in thin nonlinear films^[11]; however, the generation efficiency was much weaker compared to conventional schemes. Thereby, it remained a challenge on how to efficiently generate tunable spatially entangled photon pairs.

Recently, it was shown experimentally that metasurfaces in the form of nanostructured layers supporting optical resonances^[12–17] can boost photon pair generation^[18–21]. Furthermore, it was demonstrated that a lithium niobate metasurface featuring nonlocal guided mode resonances can generate spatially entangled photon pairs in a broad angle range^[20]. In this work, we provide a theoretical description of the two-photon state generated from metasurfaces and a quantitative study on its degree of spatial entanglement, based on the Schmidt decomposition approach. By taking advantage of the strong angular dispersion of the nonlocal metasurface, we predict that the emission pattern of the photon pairs can be varied by simply changing the wavelength or beam size of the pump laser. This represents a simple and effective way to generate photon states with tunable spatial entanglement from ultrathin metasurfaces, which may find applications in advanced quantum imaging and communications.

We consider the metasurface design following the experimental platform described in Ref. [20]. As shown in Figs. 1(a) and 1(b), it is based on an $x$-cut lithium niobate (${\mathrm{LiNbO}}_3$) thin film (thickness $h_{\mathrm{LN}}=304\mathrm{nm}$) covered by a silicon dioxide (${\mathrm{SiO}}_2$) grating with a period $a=890\mathrm{nm}$, width $w_{\mathrm{gt}}=550\mathrm{nm}$, and thickness $h_{\mathrm{gt}}=200\mathrm{nm}$. The optical axis of the ${\mathrm{LiNbO}}_3$ and the grating are along the in-plane $z$ direction. In the SPDC process, a pump photon splits into two photons, called signal and idler photons. This metasurface supports nonlocal guided mode resonance at the signal and idler wavelengths around 1570 nm, which facilitates an increased density of states and boosts the SPDC process.

Importantly, due to the subwavelength thickness of the metasurface, the longitudinal phase matching is relaxed, allowing SPDC in a broad anglular range, as shown in Fig. 1(c). The frequencies and wave vectors of the photons satisfy the following energy and transverse phase-matching conditions, as illustrated in Figs. 1(d) and 1(e), $${\omega }_s+{\omega }_i={\omega }_p,$$(1)$${\mathbf{k}}_{\perp ,s}+{\mathbf{k}}_{\perp ,i}+m·2\pi /a\widehat{\mathbf{y}}={\mathbf{k}}_{\perp ,p},$$(2)where ${\omega }_p$, ${\omega }_s$, and ${\omega }_i$ are the angular frequencies of the pump, signal, and idler photons, ${\mathbf{k}}_{\perp ,p}$, ${\mathbf{k}}_{\perp ,s}$, and ${\mathbf{k}}_{\perp ,i}$ are their transverse wave vectors in the $y–z$ plane, and $m$ is an integer. For a specific pump, matching these two conditions means that if the frequency or angle of one photon is known, the information for another photon is determined automatically. This essentially leads to the frequency and spatial correlations of the photon pairs. In the following, we consider the case of $m=0$ corresponding to the experimental conditions in Ref. [20]. For grating periods larger than the signal/idler wavelengths and allowing for high-order diffraction with nonzero $m$, one can anticipate new features in spatial entanglement, which goes beyond the scope of the current paper. The resulting two-photon state has a wave function, $$\begin{gathered} |\mathrm{\Psi }\rangle =\int \mathrm{d}{\omega }_s\mathrm{d}{\omega }_i\mathrm{d}{\mathbf{k}}_{\perp ,s}\mathrm{d}{\mathbf{k}}_{\perp ,i}S({\omega }_s+{\omega }_i,{\mathbf{k}}_{\perp ,s}+{\mathbf{k}}_{\perp ,i}) \\ ·\mathrm{\Xi }({\omega }_s,{\mathbf{k}}_{\perp ,s};{\omega }_i,{\mathbf{k}}_{\perp ,i})|{\omega }_s,{\mathbf{k}}_{\perp ,s}\rangle |{\omega }_i,{\mathbf{k}}_{\perp ,i}\rangle , \end{gathered}$$(3)where $S$ is the normalized frequency-angular spectrum of the pump light. $\mathrm{\Xi }({\omega }_s,{\mathbf{k}}_{\perp ,s};{\omega }_i,{\mathbf{k}}_{\perp ,i})$ is the SPDC efficiency function.

The wave function is controllable by the spectrum of the pump, similar to the case in conventional nonlinear crystals^[10]. Note that this control is weighted by the SPDC efficiency function. The difference in metasurfaces is that the SPDC is free from longitudinal phase matching, and its efficiency is fully determined by the frequency and angular dispersion of the optical resonances. This feature enables control of the spatial properties of the photon pairs by engineering optical resonances supported by the metasurface. For example, in the proposed metasurface with angular-dependent nonlocal resonances, the photon-pair emission pattern is very sensitive to the frequency and angular spectra of the pump. In contrast, the SPDC efficiency in bulky nonlinear crystals is limited by the longitudinal phase matching, being proportional to $\mathrm{sinc}(|\delta {\mathbf{k}}_{\parallel }|L/2)$, where $L$ is the thickness of the crystal and $\delta {\mathbf{k}}_{\parallel }={\mathbf{k}}_{\parallel ,p}-({\mathbf{k}}_{\parallel ,s}+{\mathbf{k}}_{\parallel ,i})$ is the longitudinal phase mismatch^[22,23].

The theoretical modeling of quantum photon-pair generation through SPDC in quadratically nonlinear metasurfaces can be performed through the general Green’s function formalism^[24]. A mathematically equivalent approach, which can be more convenient for numerical modeling, is based on the quantum-classical correspondence [Fig. 2(a)], where the SPDC efficiency is derived from its classical reverse process called sum-frequency generation (SFG)^{[18,20,25–28]}, $$\mathrm{\Xi }({\omega }_s,{\mathbf{k}}_{\perp ,s};{\omega }_i,{\mathbf{k}}_{\perp ,i})=\sqrt{\frac{{\omega }_s{\omega }_i}{{\omega }_p^2{(2\pi )}^3}}{\xi }_{\mathrm{SFG}}({\omega }_s,-{\mathbf{k}}_{\perp ,s};{\omega }_i,-{\mathbf{k}}_{\perp ,i}).$$(4)

Here, ${\xi }_{\mathrm{SFG}}({\omega }_s,-{\mathbf{k}}_{\perp ,s};{\omega }_i,-{\mathbf{k}}_{\perp ,i})=A_p/(A_s{A}_i)$ is the SFG efficiency in the reverse direction, with $A_p$ being the output SFG and $A_s$ ($A_i$) being the input signal (idler) complex amplitudes of unit intensity plane waves. In this work, we are interested in the angular or spatial properties of the photons and consider a continuous-wave pump. Here we focus on the case of degenerate SPDC (${\omega }_s={\omega }_i=0.5{\omega }_p$), since the frequency spectrum of photons was found to be narrow^[20]. This is because the degenerate SPDC efficiency is much stronger than in the nondegenerate case, as we show in Section 3. The corresponding spatial wave function becomes $$|\mathrm{\Psi }\rangle =\int \mathrm{d}{\mathbf{k}}_{\perp ,s}\mathrm{d}{\mathbf{k}}_{\perp ,i}\frac{S({\mathbf{k}}_{\perp ,s}+{\mathbf{k}}_{\perp ,i})}{4\pi }{\xi }_{\mathrm{SFG}}(-{\mathbf{k}}_{\perp ,s};-{\mathbf{k}}_{\perp ,i})|{\mathbf{k}}_{\perp ,s}\rangle |{\mathbf{k}}_{\perp ,i}\rangle .$$(5)

After calculating this two-photon wave function, we can reveal the emission brightness and pattern of the SPDC and the spatial entanglement of the photon pairs. Apparently, they are controllable through the angular spectrum of the pump and the angular-dependent SFG efficiency determined by the optical resonance of the metasurface. The complexity of Eq. (5) comes from the fact that for a specific pump, the photon pairs from the metasurface can emit to all directions satisfying the transverse phase-matching conditions in Eq. (2). As a result, one needs to simulate the SFG process in all directions to simulate the full SPDC process, which is a computationally demanding task.

To efficiently model the SPDC process in the lithium niobate metasurface, we have proposed a coupled mode theory (CMT) that can accurately calculate the frequency and angular dispersion of the guided mode resonances in the metasurface (see detailed theory in the Supplementary Material of Ref. [20]). Figure 2(b) shows the CMT predicted resonance wavelengths of the metasurface at different transverse wave vectors along the $y$ direction. Based on the CMT, we can also calculate the SFG efficiencies for different $k_{y,s}$ and pump wavelengths ${\lambda }_p$. Note that the guided mode resonance and SFG efficiency has a weak dependence on $k_{z,s}$, and thus we only consider the case of $k_{z,s}=0$ here. In the SFG simulations, we calculate the generated pump with a fixed frequency (${\omega }_p$) and emitting along the normal direction (${\mathbf{k}}_{\perp ,p}=0$) such that ${\omega }_i={\omega }_p-{\omega }_s$ and $k_{y,i}=-k_{y,s}$ in each simulation. As shown in Figs. 2(c)–2(f), the CMT results show a good agreement with the full-wave simulations conducted in COMSOL Multiphysics, confirming its high accuracy in modeling SFG and the corresponding SPDC process based on Eq. (5). Importantly, the COMSOL simulation takes around 1 h to calculate 500 points of the signal wavelength, while the CMT takes only 14 ms. Note that the SFG efficiencies in Figs. 2(c) and 2(f) are about an order of magnitude higher than the ones in Figs. 2(d) and 2(e). This is because in Figs. 2(c) and 2(f), the degenerate signal and idler photons are both at resonance wavelength, while only either a signal or idler photon is resonant for the nondegenerate case in Figs. 2(d) and 2(e). Due to the quantum-classical correspondence, we expect that the degenerate SPDC efficiency is the strongest in the specific direction determined by the pump wavelength. Therefore, only the frequency-degenerate SPDC is considered in the following analysis of the spatial entanglement of the photon pairs.

We now calculate SFG at the degenerate wavelength for different wave vectors and obtain the two-photon wave function of the photons based on Eq. (5). Specifically, we consider a Gaussian pump beam with an angular spectrum $S({\mathbf{k}}_{\perp })=\exp (-{\sigma }_p^2{|{\mathbf{k}}_{\perp }|}^2/4)$, where ${\sigma }_p$ is the beam waist radius. Note that a normally incident Gaussian pump beam is effectively composed of infinitely obliquely incident plane waves with a nonzero ${\mathbf{k}}_{\perp }$. By integrating over the wave vector, we obtain the spectral brightness of the SPDC, as shown in Fig. 3(a) for different pump wavelengths and beam radii. In the calculation, an upper bound of 0.05 rad/µm has been applied for $k_y$ and $k_z$, which corresponds to a collection angle of 0.7 deg in a previous experimental setup^[20]. One can see that for each beam size, there is a peak for the SPDC at $2{\lambda }_p=1569.1\mathrm{nm}$, which is slightly blueshifted from the resonance at normal incidence. We notice that the highest brightness occurs at ${\sigma }_p=45\mathrm{\mu m}$. This is because the generation rate of photon pairs depends on the incident angle of the plane-wave pump and the maximum happens near $k_{y,p}=0.015\mathrm{rad}/\mathrm{\mu m}$, as shown by Fig. S3d in the Supplementary Material of Ref. [20]. The pump beam with a radius of ${\sigma }_p=45\mathrm{\mu m}$ has the largest weight at the wave vector components near this value.

We perform a Schmidt decomposition of the two-photon wave function to quantify the spatial entanglement. The Schmidt decomposition of the wave function vector is performed, and the related Schmidt coefficients ($s_m$) and modes are obtained by using the open source MATLAB code QETLAB. The related Schmidt numbers are calculated by ${\sum }_m{s}_m^4$ and shown in Fig. 3(b). Obviously, the Schmidt number is strongly dependent on the wavelength and beam size of the pump and can be tuned in a large range from 1.1 (weakly entangled) to 20 (strongly entangled). A larger beam radius leads to a larger Schmidt number. This is reasonable, since a plane-wave pump with an infinite beam radius has a single transverse wave vector, resulting in photon pairs with fully deterministic relation between the signal and idler wave vectors according to the transverse phase matching, and each of a continuum number of such pairs represents a Schmidt mode. For the beam radius of 200 µm, the maximum Schmidt number is found to be 20 when the pump wavelength is 784.15 nm, as marked by point $A$ in Fig. 3(b). At the same pump wavelength, the Schmidt number is only 1.1 for a pump radius of 5 µm, which is point $B$ in Fig. 3(b). The related SPDC emission patterns for these two points are shown in Figs. 3(c) and 3(d), respectively.

To gain more insight into the spatial entanglement in Fig. 3(b), we focus on the Schmidt decomposition results at points $A$ (${\lambda }_p=784.15\mathrm{nm}$, ${\sigma }_p=200\mathrm{\mu m}$) and $B$ (${\lambda }_p=784.15\mathrm{nm}$, ${\sigma }_p=5\mathrm{\mu m}$). Figure 4 depicts the Schmidt coefficients of the lowest 50 Schmidt modes and the mode distributions of the lowest four modes for point $A$. As can be seen, all the 50 Schmidt coefficients are nonzero, and the first four modes have similar magnitudes, indicating a high degree of spatial entanglement. The results for point $B$ in Fig. 3(b) are shown in Fig. 5. The Schmidt coefficients show a fast drop for higher-order modes and are close to zero for the modes with a mode number over 3, confirming a weak entanglement at point $B$.

In conclusion, we have theoretically investigated the generation of spatially entangled photon pairs via SPDC from a lithium niobate metasurface. By quantifying the spatial entanglement with the Schmidt number, we have shown that the degree of spatial entanglement can be diversely tuned by the pump laser wavelength and beam size. The capability to realize an arbitrary degree of spatial entanglement can find applications in quantum imaging, whose resolution is related to the spatial entanglement of the photon source. In the future, pump beams with tailored spatial and temporal profiles and metasurfaces with different optical resonances can be further explored to tune both the emission pattern and spatial entanglement.