Entangled Photon-pair Generation in Nonlinear Thin-films

We develop a fully vectorial and non-paraxial formalism to describe spontaneous parametric down-conversion in nonlinear thin films. The formalism is capable of treating slabs with a sub-wavelength thickness, describe the associated Fabry-P\'erot effects, and even treat absorptive nonlinear materials. With this formalism, we perform an in-depth study of the dynamics of entangled photon-pair generation in nonlinear thin films, to provide a needed theoretical understanding for such systems that have recently attracted much experimental attention as sources of photon pairs. As an important example, we study the far-field radiation properties of photon pairs generated from a high-refractive-index nonlinear thin-film with Zinc-Blende structure, that is deposited on a linear low-refractive-index substrate. In particular, we study the thickness-dependent effect of Fabry-P\'erot interferences on the far-field radiation pattern of the photon pairs. We also pay special attention to study of entanglement generation, and find the conditions under which maximally polarization-entangled photon pairs can be generated and detected in such nonlinear thin-films.


Introduction
The optical process of spontaneous parametric downconversion (SPDC) in  (2) -materials [1], in which a pump (p) photon of wavelength  p can split into a pair of entangled signal (s) and idler (i) photons of lower energy, is an instrumental tool in optical quantum technologies.Predominant applications are quantum communication [2], quantum cryptography [3], [4], quantum imaging [5], [6] and sensing [7], and testing fundamental quantum effects [8].The broad applicability of SPDC is due to its efficiency and versatility in generating pairs of entangled photons with high control over various optical degrees of freedom.
Nonlinear systems involving SPDC constitute to date the dominant approach to generate entangled photon-pair states and promise tangible real-world applications in the near future [9].Entangled photon pairs have been generated in single nonlinear bulk crystals [10], crossed bulk crystals [11], or bulk crystals in a Sagnac loop [12].More recently, several experiments of photon-pair generation in much smaller nonlinear systems such as single GaAs nanowires [13], dielectric nanoantennas [14], and metasurfaces with subwavelength thicknesses [15]- [18] have been demonstrated.In particular, there has been a recent growing interest in pair generation in sub-wavelength-thin nonlinear crystals like lithium niobate and "zinc-blende" structures like gallium arsenide (GaAs) or gallium phosphide (GaP) [19]- [21], as well as in van der Waals [22] or transition metal dichalcogenide (TMD) crystals [23], [24].Such nonlinear thin sources are interesting because they are not restricted by the longitudinal phase-matching condition; hence they are capable of generating photon pairs in a very wide spectral and angular range.This also motivates their use in quantum imaging applications, for pushing the resolution of such schemes to the limit [25], [26].
However, to the best of our knowledge, there has not yet been a comprehensive theoretical analysis of photon-pair generation in a thin and unstructured nonlinear slab, where the dynamics of entanglement generation is investigated and the conditions for generating a maximally entangled photon-pair state are identified.Such an analysis, which would be facilitated by a general theoretical framework, is an essential missing link for further development of extremely thin and broadband sources of entangled photon pairs.This theoretical framework should incorporate the Fabry-Pérot-type interferences due to multiple reflections of signal, idler and also pump photons within the nonlinear slab.Moreover, such a theory should be able to treat absorptive nonlinear layers, particularly in light of the more prevalent use of TMD systems which provide an enormous enhancement of the nonlinear coefficients near their excitonic lines that at the same time can also cause enhanced material absorption in the system [27].Importantly, the formalism should also be able to describe the directional, spectral, and polarization properties of the generated photon pairs, including emission angles up to 90 • .
In this work, we develop such a theoretical framework for describing SPDC in nonlinear thin-films, and use it to study the far-field properties of the generated pairs.Our formalism is based on the Green's function (GF) quantization method [28].In this method, SPDC is described using a quantization scheme of the electric field operator that utilizes the classical GF of the system and local bosonic excitations in the medium-assisted field [29], [30].The GF, as a classical quantity, contains all the linear properties of the system.In our work, we construct this classical GF using a developed theory for the GF of multilayered optical systems [31].Due to the generality of the GF method, we develop a comprehensive, fully vectorial framework capable of treating any thickness of the slab, ultra-thin or thick.It can further be used to incorporate lossy layers, near-and far-field radiation in the non-paraxial regime, as well as the generation of guided photonic modes inside the source.In the current work, we focus on far-field properties of the photons that are generated outside the slab.Furthermore, our model allows keeping track of any polarization and directionality effects in the pair-generation process, which we will exploit in this paper to reconstruct the polarization states of entangled photon pairs.
In this work we are treating SPDC in the low-gain regime of photon-pair generation.However, recently, the application of the GF quantization method has also been extended to treat the high-gain regime of the interaction [32], where the same classical GF can be used for description of high-gain effects.This further emphasizes the versatility of the GF approach in description of parametric down-conversion under different conditions.We also note that there are other methods for description of SPDC [1], [33]- [37], some of which were also developed for description of low-gain SPDC in multilayer systems [35]- [37].Nonetheless, the GF method provides a general yet unified formulation for treatment of such systems, which we employ to perform a detailed investigation of SPDC and identify entanglement properties in thin-films.Specially, the GF method can naturally include the effect of internal losses, and also it provides a way towards description of high-gain effects [32], which allows for study of thin-film systems in more advanced scenarios beyond this work.
In this investigation, we will utilize numerical calculations to explore the far-field characteristics of photon-pair generation via SPDC within a GaAs thin film.This is part of a larger category of nonlinear materials characterized by a zinc-blende structure, a type of nonlinear material with a cubic lattice and point group 43 m.This crystal structure is shared by various other prominent nonlinear materials, including III-V semiconductors such as gallium phosphide (GaP) and indium phosphide (InP), as well as II-VI semiconductors like zinc telluride (ZnTe), zinc selenide (ZnSe), and zinc sulfide (ZnS).Due to the cross-polarized nature of their nonlinear tensor, these materials are naturally predisposed for generating polarization-entangled Bell states [21], [38].This holds significant promise for advancing various quantum information processing applications, ranging from quantum cryptography to quantum teleportation.By investigating the unique properties of photon-pair generation in GaAs and related nonlinear materials, we aim to contribute to the broader understanding and potential utilization of entangled photon pairs in cutting-edge quantum technologies.
The manuscript is organized as follows.In Section 2 we introduce a general framework to calculate the far-field joint detection probability of photon pairs generated in a nonlinear slab in a layered geometry using the GF formalism, and compute the far-field joint probability distribution.In Section 3, we explicitly write the general expression of the GF's of a dielectric slab that contains multiple reflections of the signal and idler fields.In Section 4, the classical pump beam is treated.In Section 5, we numerically compute the far-field radiation for a nonlinear source made of GaAs deposited on a silicon dioxide (SiO 2 ) substrate, and study its far-field properties.In Section 6, we use a quantum tomography method to extract the polarization state of the generated photon pairs and show that one can generate maximally polarization-entangled photon pairs in such a simple system.Lastly, we summarize and conclude in Section 7.

Far-field joint spatial probability in Fourier domain
Consider a general nonlinear thin-film photon-pair source, as depicted in Figure 1.It comprises three layers with relative permittivities  1 ,  2 , and  3 , where the nonlinear material (of thickness a) constitutes medium 2. Typically, medium 3 is a linear substrate and medium 1 is a linear cladding covering the nonlinear material (often air).Here, the system is illuminated by a pump beam from medium 3 propagating along the positive z-direction, and the down-converted photon pairs can be detected in media 1 or 3.It is important to note that the generated photon-pair amplitudes undergo multiple reflections inside the nonlinear slab.Although this paper does not delve into the dynamics of photon pairs that could be generated in the guided slab modes of medium 2, our formalism can handle such scenarios too.This is left to future works.In this work, we focus on studying the dynamics of photons generated into the free propagating modes of medium 1.
Using the GF method, we analyze SPDC in a nonlinear structure excited by a single-frequency pump beam dence detection rate for a pair of signal and idler photons of different frequencies as a function of the spatial coordinates follows [28], [38]: (r) The detection rate,  ≡ , is the number of photon pairs given per units of solid angles dΩ i (for the idler photon) and dΩ s (for the signal photon), and per unit of signal-photon angular frequency  s .The idler angular frequency is fixed by conservation of energy to  i =  p −  s .The detection rate is also a function of the idler and signal photons' detection polarization, which are along the unitary e s and e i directions.The , , ,  s and  i indices run over the x, y, and z cartesian coordinates.G ij (r, r ′ , ) are the tensor-components of the electric GF,  (2)    (r) are the components of the second-order nonlinear tensor and the integral runs over the volume of nonlinearity V characterized by the position vector r.r i and r s are the positions for the detection of the idler and signal photons in the far-field, respectively, with n i and n s being the refractive indices of the detection medium at the idler and signal frequencies.
In our system, one can take advantage of the homogeneity in the x and y direction to work in the spatial frequencies domain.This approach offers an intuitive understanding of the far-field detection in terms of the angular modes of the signal and idler photons, simplifying numerical calculations.To do this, the GF, which is a solution of , with (r, ) being the relative permittivity of the system, can be represented as an inverse Fourier transform of the form [39], where q = k x x + k y ŷ and r ⊥ = x x + yŷ are the two-dimensional real wave-vector and spatial vector in the x-y plane, respectively.We define Similarly, for a classical pump field treated in the undepleted pump approximation we have, ( Then, by introducing Eqs. ( 2) and (3) into (1), we can separate the spatial integral into dr → dzdr ⊥ .After rearranging terms, we can solve for dr ⊥ the integral dr ⊥ e i(q p −q s −q i )⋅r ⊥ = (2) 2 (q p − q s − q i ), (4) which results into the joint detection rate 4 -E. A. Santos et al.: Photon-pair generation in nonlinear thin-films (2)   (z) E p, (q s + q i , z) Here, Eq. ( 6) is understood as a joint angular probability (JAP) amplitude.Equation ( 5) represents the detection rate for signal and idler photon-pairs generated through SPDC in a transversally homogeneous nonlinear source, accounting for spatial coordinates, frequencies, and detection polarization.Notably,  in the form of (5) resembles a Fourier transform integral, enabling Fourier analysis for interpreting the probability in the more intuitive angular domain.
However, since our focus is mainly on the far-field properties of the generated photon pairs, we will apply a far-field approximation to the probability rate similar to the classical Fraunhofer approximation for field diffraction [40].We aim to detect both photon pairs in the far field in medium 1 with refractive index n 1 .Here, k 1 =  c n 1 is the wave number of plane waves in medium 1, and denotes the z-component of its wave vector k 1 = q + k z,1 ẑ.
Note that k x and k y are the transversal components of the wave-vector for a photon generated in medium 2, which are also conserved in medium 1 due to the system's transversal invariant nature.The far-field approximation requires joint detection of photons at sufficiently long distances from the source, i.e. r s ≫ r and r i ≫ r.Additionally, we assume a finitesized pump beam, implying that detection occurs at a much greater distance compared to the extent of the pump illumination across the nonlinear slab (see Appendix A for details).
Under this approximation the joint detection rate in Eq. ( 5) is simplified to (see Appendix A) Here, q s = k x,s x + k y,s ŷ, where k x,s = k 1s sin s cos s and k y,s = k 1s sin s sin s , refer to the transverse k-vector of the signal photon detected in the far-field at propagation angles { s ,  s }, and similarly for the idler photon.Thus, Eq. ( 7) describes the joint detection rate in the far-field as a function of the detection angles { s ,  i ,  s ,  i } in spherical coordinates.Note that k z,1s and k z,1i in Eq. ( 7) correspond to the detection medium 1.If the detection medium changes, these factors should be adjusted accordingly (see Appendix A).Also note that,  in Eq. ( 7) is z-independent, where (q s , q i , z s , z i ,  s ,  p −  s ) = (q s , q i ,  s ,  p −  s )e ik z,1s z s e ik z,1i z i (see Appendix A).
In the far-field, the joint detection rate is proportional to the squared absolute value of  multiplied by the signal and idler z-components of the wave vector, k z,1s and k z,1i .
In the paraxial regime (q ≪ k 1 with q = |q|), k z,1s and k z,1i tend to k 1s and k 1i , resulting in the far-field joint detection rate being directly proportional to the squared absolute value of .This calculation extends the concept of classical Fraunhofer approximation [40] to two-photon jointdetection.As our model is capable of handling emission angles up to 90 • , it highlights the importance of the factors k z,1s and k z,1i in the nonparaxial regime.The importance of the longitudinal wave-vectors in nonparaxial description of quantum-light propagation was also shown in studying the fundamental resolution limit of quantum imaging schemes [25].Equations ( 5)-( 7) constitute the main results of this section, to be numerically solved in Section 5.

Green's function of a dielectric slab
The optical properties of the system at the signal and idler wavelengths are fully contained in their respective Green's functions (GFs).To implement Eq. ( 7), we employ a GF suitable for the multilayered scheme depicted in Figure (1).
In this scenario, the GFs for signal and idler, ↔ G(r, r s ) and ↔ G(r, r i ) connect a point r within the  (2) -nonlinearity region (medium 2, see Eq. ( 1)) to the detection positions r s and r i in medium 1.
If a plane wave of the form E(r) = Ee ik⋅r , with c n satisfies Maxwell's equation in a homogeneous medium of refractive index n, two possible wave vectors k for a given q are obtained, representing upward and downward propagation.These wave vectors take the form k + = q + k z ẑ (upward) and k − = q − k z ẑ (downward).When n is real and q > n∕c, k z becomes imaginary, resulting in an evanescent wave.Additionally, at a given q, a solution can be either s-polarized (transverse electric) or p-polarized (transverse magnetic), defined by unit polarization vectors ŝ and p as These vectors are normalized according to ŝ ⋅ ŝ = 1 and p± ⋅ p± = 1.Note that ŝ is real, rests on the x-y plane, and is medium independent, while p can be complex and is medium dependent.Also, ŝ does not depend on the direction of k, while p changes for upward-or downward-generated fields.
Using the sand p-polarized fields, the GF in Fourier domain for such a system, where the field is generated in medium 2 and propagates outside the slab to medium 1, can be constructed as [31], [39 21 (q, z ′ , ) ŝŝ where T 21 (q, z ′ , ) is the generalized transmission coefficient of the sand p-polarized fields generated from the source plane at z ′ ∈ {−a, 0} and transmitted to medium 1 (see Appendix B for analytical expression).The unit vector ŝ is parallel to the interface and p is perpendicular to the wave vector k and ŝ.The subscript numbers in k z indicate the corresponding medium.
The GF in Eq. ( 10) describes from right to left (apart from the phase factor e ik z,1 z ) the contributions of a downward p-polarized wave generated in medium 2 that transform into an upward p-polarized wave in medium 1, an upward p-polarized wave in medium 2 that transform into an upward p-polarized wave in medium 1 and both the contribution of the upward/downward s-polarized waves generated in medium 2 that transfers to medium 1.All of them after multiple reflections and transmission, whose information is contained in the coefficients T 21 .

Pump field with multiple reflections
The pump field undergoes diffraction and multiple reflections within the slab.Illustrated in Figure 1, the pump originating from medium 3 propagates to z = −a, enters the slab, and reflects multiple times.Consequently, the amplitude of the pump field within medium 2 varies with z and is influenced by the slab's thickness a.
In many theoretical calculations, it is common to assume a pump beam polarization parallel to the interface (e.g.x-polarized) for SPDC modeling.However, true transverse electromagnetic waves are idealizations and physical beams of light also contain a longitudinal polarization component.This longitudinal component may influence the photon-pair generation process if the material properties permit such excitation.For our calculations, we consider a pump field polarized in the x-direction right before the interface (z = −a) with a longitudinal component in the z-direction, given by (11) where k z,3 p represents the z-component of the pump wave vector in medium 3, and E p,z is calculated from E p,x using the transversality constraint k p ⋅ E p = 0 (see Appendix C.2). Taking y ẑ, which allows us to express the pump field in the form where we take ê = x for k x = k z = 0. U p (q, z = −a) is the angular spectrum of the pump field at the interface.We assume it follows a Gaussian distribution, where  is the width of the Gaussian angular spectrum.
The condition q 2 ≤ k 2 3 p ensures that Eq. ( 3) comprises only propagating waves in medium 3.Moreover, the spectrum width influences the contribution of the z-component of the pump beam: larger  results in greater E p,z .This is especially relevant in tightly focused pump beam scenarios (refer to Appendix C.1 for details).
At the interface, the pump beam can also be described using sand p-polarized fields.After transmission into the slab, the z-dependent pump beam can be calculated as (see Appendix C.2 for details) ) , (14) where Q (s) (q, z) is the transmission coefficient from medium 3 to medium 2, while E (s) p (q, z = −a) and E ( p) p (q, z = −a) are sand p-components of the pump before the slab.Unit vectors ŝ, p2+ , and p2− are defined in Eqs. ( 8) and ( 9) for medium 2.
For degenerate photon-pairs, this motivates investigating  (far) in terms of propagation angles { s ,  s } for the signal photon, where the idler is detected symmetrically at { i =  s ,  i =  +  s }.This configuration is shown in Figure 2(a) and is called -symmetric throughout the manuscript.The -symmetric setup has already been proven crucial for entangled state generation in nanoresonators [38] and will also show to be essential for thin films in this work.The coincidence rate in the -symmetric configuration is shown in Figure 2(b) for an ultra-thin source of thickness a = 0.01 p , emphasizing again that both photons are detected in medium 1.Here, the  (2)  yzx and  (2)   zyx dominate the pair generation process due to weak yand z-polarized and dominant x-polarized components of the pump field inside the slab.Note that y-components of the pump field are created inside the slab due to refraction, despite being zero outside (see Eq. ( 14)).Hence, all components of the GaAs nonlinear tensor participate in the generation process.
Yet, the dominance of  (2)  yzx and  (2)   zyx components implies that signal and idler fields within the slab are mainly generated by yand z-polarized point-like nonlinear sources, resulting in no radiation in the y-z plane, as observed in Figure 2(b).Additionally, emission at large angles is significantly reduced, despite the GaAs nonlinear tensor favoring stronger emission near  s = 90 • for such an ultra-thin crystal.This reduction is mainly due to the strong decrease in transmittance from inside to outside of the slab at larger emission angles, ultimately reaching zero for outside emission angles of  s = 90 • due to total internal reflection experienced by photons inside the slab.
In Figure 2(c), we display cuts along the xz-plane for various subwavelength thicknesses of the slab, a = 0.01 p , a = 0.1 p , a = 0.5 p , and a = 0.8 p . Figure 2(d) illustrates xz-plane cuts for larger slab thicknesses, a = 3.61 p , a = 3.67 p , a = 3.72 p , and a = 3.82 p , which are special points for understanding the Fabry-Pérot dynamics, marked in Figure 2(e)-(g).All graphs are normalized to one for better observation of angular dependencies.Figure 2(c) reveals minimal angle variations in joint emission radiation for different subwavelength thicknesses, primarily influenced by the nonlinear tensor structure and angle-dependent transmission at interfaces.In contrast, Figure 2(d) demonstrates more significant angular changes for larger slab thicknesses.This is mainly due to the Fabry-Pérot effect, which enhances or reduces radiation in some angles due to constructive or destructive interferences, respectively.Notably, the angular pattern is highly sensitive to slab thicknesses in this regime.To better understand the dynamics caused by Fabry-Pérot effect, we show in Figure 2(e) the angle of maximum intensity, (f) the intensity at  s = 45 • , and (g) and the angular width of the emission pattern Δ, respectively, all as a function of the slab thickness.In Figure 2(e), the maximum emission angle exhibits a growing oscillatory pattern with increasing the thickness.The oscillatory behavior occurs around a mean angle of ∼45 • .At the detection angle of  s = 45 • , we plot in Figure 2(f), the detection intensity (normalized to maximum) as a function of the slab thickness.We clearly observe an etalon-like effect caused by the FP interferences, where peaks of intensity repeat periodically with increasing slab thicknesses, with a period of Δa ∼ 0.29 p .This period of oscillation is the same on all three plots of Figure 2(e)-(g).The angular uncertainty of emitted photons exhibits a similar thickness-dependent variation.We define the angle uncertainty Δ as the difference between the angles at 50 % of the maximum intensity (full width at half maximum).In Figure 2(g), we again observe increasing oscillations for thicknesses above a =  p , while it remains relatively constant for subwavelength thicknesses, as also seen in Figure 2(c).
It can be verified, that the peaks and valleys in Figure 2(f), correspond to constructive and destructive interference conditions, satisfying the relations Δ = m s and Δ = (m + 1 2 ) s , respectively, where m is a positive integer.Δ = 2an 2s cos( 2s ) is the optical path difference (OPD) between successive plane waves that are transmitted to medium 1 (being air) after a round trip in the slab [43].
2s is the angle in medium 2 such that it yields  1s = 45 • in medium 1 after refraction.The period of the resonance peaks is found by Δa =  s 2 n 2s cos( 2s ) , which results in Δa ≈ 0.29 p , which matches the numerical results.These points of constructive and destructive interferences for  1s = 45 • also coincide with peaks and valleys of Δ, see in Figure 2(g).We point out that the angles of constructive interference oscillate around the  1s = 45 • , as the thickness is increased, which results in an oscillation of the angles of maximum intensity around  1s = 45 • in Figure 2(e).This can also can be seen in Figure 2(d), where at thicknesses of a = 3.61 p or a = 3.72 p , the directionality plots are maximized at angles ∼58.5 • or ∼27.5 • , respectively.
A crucial observation is that due to the inherent losses in GaAs at  p = 500 nm, there are constrictions on longitudinal phase-matching effect.Calculating the pump's decay length inside the GaAs slab, denoted as When compared with the coherence length of our sys- , we obtain L c = 0.6 p for joint detection in medium 1 at  1s =  1i = 45 • .Notably, the decay length is smaller than the coherence length ( p < L c ), which is important, as it limits the range of impact of longitudinal phase matching with increasing slab thicknesses.Consequently, the FP effect at the signal and idler wavelengths has a dominant effect in creating the oscillatory behavior observed in Figure 2(e)-(g), and not the longitudinal phase matching.Hence, the longitudinal phase matching, only affects the height of the first few peaks in Figure 2(f), while subsequent peak intensities stabilize.

Non-degenerate photon-pair detection
Now we consider a frequency non-degenerate scenario for photon-pair detection.Unlike the degenerate case, the symmetric configuration is not expected to dominate in the non-degenerate case.For a normally incident plane-wave pump, with q p = 0, the emitted pairs satisfy the transversal phase-matching condition of q s = −q i .This right away sets  i =  s +  for signal and idler photons in a pair.Then, taking  as the emission angle with respect to the forward z-direction we get k 2 x + k 2 y = k 2 sin 2  for both signal and idler.Using the degeneracy factor r defined as r ≡  s ∕ i , we get as the relation between the angles of emission for signal and idler photons in a pair.In medium one, we have n i = n s = 1.
In Figure 3(a), we calculate the joint detection probability with a fixed signal detector at  s = 45 • ,  s = 180 • for nondegeneracy factors r = 0.8, 1, 1.2, and 1.5 in the ultra-thin case (a = 0.01 p ).Here, smaller idler wavelengths (r > 1) tend to emit the idler closer to the optical axis, while larger idler wavelengths (r < 1) push the idler away from the optical axis.Maximum emission angles calculated using the relation between angles are  i (r = 0.8) = 62.11Here we can see that a better correlation in the non-degenerate case can be achieved by detecting at small angles, however this comes at expenses of the efficiency of the process.Additionally, emission at ∼45 • , dominant in the degenerate case, is no longer prevalent in the nondegenerate scenario, which is now shifted by transverse phase-matching.In the depicted scenario in Figure 3

Absolute pair generation rates
A critical question revolves around the overall efficiency of the pair-generation process.The total number of pairs that can be collected by a lens of numerical aperture NA, is obtained by integrating the differential pair-rate  ≡ d 4 N dtdΩ s dΩ i d s over the entire signal and idler solid angles that fall within the angle  = arcsin(NA) where dΩ = sindd, then summing over all possible detected polarization combinations, and finally integrating over the frequency range of the detected signal photons.To do this, we consider a 0.59 p ≈ 295 nm for the GaAs slab thickness, which corresponds to the maximum intensity peak in Figure 2(f).We assume a pump power of 1 mW spread over the Gaussian beam of width W = 3 μm with spectrum defined in Eq. (13), where the relation W = 2∕ is satisfied in the paraxial regime.The incident pump power P pump is related to the amplitude of the Gaussian spectrum through A = √  0 cP pump W 2 , in the paraxial pump limit.For ⟨100⟩ GaAs we take  (2)  0 ≈ 300 pm/V as an approximate value for its nonlinearity coefficient based on different frequencydependent measurements [44].
After the solid-angle integrations, we obtain the quantity d 2 N pair ∕(dtd s ), where we calculate this for the fixed signal frequency corresponding to the degenerate case, at which  s =  i = 1 μm.This quantity is dimensionless and represents the number of photon pairs per second that one can collect per units of signal photon frequency.For NA values of 0.6, 0.8 and 0.9 we obtain d 2 N pair ∕(dtd s ) ≈ 5.2 × 10 −13 , 1.7 × 10 −12 , and 2.6 × 10 −12 , respectively.
Assuming that the generation efficiency is more or less constant for about 10 nm of bandwidth for the signal photons (i.e. from 1 to 1.01 μm wavelength), which is a good approximation based on the result in Figure 4(b), we can multiply the derived d 2 N pair ∕(dtd s ) by a corresponding Δ ≈ 1.87 × 10 13 .Then, we find the total rate of photon pairs with a bandwidth of 10 nm for the signal photons, collected with NA values of 0.6, 0.8, and 0.9, to be about 10, 32, and 49 Hz.These numbers are comparable to SPDC measurements in sub-micron-thick thin films [20], [23], once the measured rates are scaled appropriately for experimental collection efficiencies and detection bandwidths.Finally, we note that changing the width of the pump beam from 3 μm to 6 μm (carrying the same 1 mW of power), does not change the collected photon-pair rate appreciably, where with an NA = 0.9 for the collection lens and 10 nm bandwidth for the signal photon, we find a photon-pair rate of 51 Hz.

Quantum polarization tomography
So far we have examined the mechanism of photon-pair generation by a thin-film nonlinear source, concentrating on the far-field directionality properties as functions of the slab thickness and photon-pair wavelengths.For this purpose, our investigation has focused only on polarizationindependent joint-detection.However, we are also interested in the polarization state of the photon pairs, which we will study in this section.Although our model based on the GF quantization approach does not directly predict the quantum state of the photon pairs, but rather predicts detection rates, it still allows us to make these detection rates dependent on the specific wavelength, position, and polarization of the detected photon pairs.Utilizing this, one can determine the density matrices of the biphoton polarization states following a tomographic method [38], [45].There, the density matrix ρ, expressed in the {|HH⟩,  [38], [45].The details of the tomographic approach can be found in the Appendix D.
Using this method, we examine the density matrix for several scenarios, starting with an ultra-thin slab (a = 0.01 p ) in a -symmetric configuration detection ( s = 45 • ,  s = 0) for r = 1.The pump is the same as in Section 5.
In Figure 4(a), we display the real and imaginary parts of the density matrix ρ representing the polarization state.
Remarkably, under these conditions, the density matrix corresponds to a maximally polarization-entangled state For every photon-pair state, we assess the level of polarization entanglement by determining the Schmidt number, K.A value of K = 1 indicates a completely untangled polarization state, while a value of K = 2 means a fully entangled state (see Appendix D).We investigate entanglement for two thicknesses (a = 0.01 p and a = 0.1 p ), including nondegeneracy, in Figure 4(b) at the same -symmetric configuration.We plot both the Schmidt number (K) and the detection rate (normalized to one) as a function of the wavelength of the signal (bottom axis) and its corresponding idler photon (top axis).We observe that high-quality entanglement is preserved along an extremely broad range of wavelengths in both thicknesses, despite the joint-detection probability dropping at values r ≠ 1 due to moving away from the transversal phase-matching condition (as in Figure 3 , with a = 0.01 p and a = 1 p .Interestingly, while the angles of maximum entanglement (marked with red " * ") and maximum joint-probability (marked with green "x") coincide for r = 1, they move in opposite directions for r = 1.5.Specifically, the angle of maximum detection rate approaches the optical axis rapidly, following transverse phase matching, whereas the angle of maximum entanglement shifts slightly away from the -symmetric angle.Notably, although entanglement slightly degrades in the non-degenerate -symmetric case (the values of K ∼ 1.98 and 1.99 in Figure 5(b) and (d)), there actually exists a close by angle at which the entanglement is maximum again at K ∼ 2.Moreover, the angles of maximum detection rate and entanglement are significantly apart for the non-degenerate case, indicating a trade-off between efficiency and entanglement with nondegenerate photon pairs, at least under the stated pumping conditions.Additionally, maximally polarization-entangled states exhibit strong robustness to changes in thickness and wavelengths, which parallels findings in Ref. [38], where similar responses were observed in Bell state generation in GaAs nonlinear nanoresonators.

Effects from a finite collection angle
Throughout this work, we focused on an angularly-resolved analysis to best understand the emission dynamics of SPDC from thin films, which experimentally corresponds to detecting photons that are collected over infinitely narrow emission angles.In experimental scenarios, collection over extremely narrow angular ranges is not feasible, especially because it drastically reduces the observed number of photon pairs.Commonly, the two-photon emission from such nonlinear thin films is collected with an objective/lens that collects a finite range of emission angles [21]- [23].In the following, we investigate how collection over finite emission angles can affect the degree of entanglement.
We focus on a case with slab thickness of a = 0.01 p , detecting only frequency-degenerate signal and idler photons.Our numerical analysis is performed in a way that resembles a simplified experimental scenario, in which a lens with numerical aperture NA is collecting all signal and idler photons going in any pairs of emission angles within the cone of angle  = arcsin(NA) (see Appendix D), and the tomography is performed by integrating over the detection probabilities from these collected pairs.Firstly, we find that the polarization density matrices for the two-photon states are no longer fully pure as was before, but have a purity P < 1, as shown in Table 1 for different values of collection NA and a fixed pump width of W = 3 μm.This is due to the pump beam having a finite width, which makes the angular correlation between pairs not perfect (see Figure 3(a) for r = 1), and hence in the process of integrating over the different collected angles, the information about these correlations is lost, which results in a mixed state.We verify this by finding the purity for a fixed NA and varying values of pump width, as shown in Table 2.Here we see that larger pump widths, which produce higher angular correlations, result in a higher purity.Also, the fact that purity rises with larger NA in Table 1, is attributed to the fact that, since in our system the pairs are mainly radiated around  ∼ 45 • (as shown in Figure 2), going towards higher NAs results in capturing the emission at angles near 45 • , which have a dominant effect in the sum over all collected angles, thus improving the systems performance.With a mixed state, a better measure for entanglement is the concurrence, C, a quantity ranging between C = 0 for separable and C = 1 for fully entangled states [46].Here, we also observe a C < 1, yet quite close to 1, which means that the entanglement degree is also slightly reduced.This can also be explained, as we showed in Figure 5, because entanglement is only perfect for pairs of -symmetric angles, yet the fact that the pump has a finite width results in generation of degenerate pairs that deviate from this condition (Figure 3(a) for r = 1).This means that we also collect pairs that are not perfectly polarization entangled, and this reduces the concurrence of the total collected state.This is verified from the data in Table 2, which shows that the concurrence increases with a wider pump beam.Concurrence also improves with collection of larger angles, as seen in Table 1, which we attribute to the same reason that purity increases with higher NA.Overall, we can see that with reasonable values of NA and pump widths, high values for purity and concurrence can be achieved, which can be improved consistently with increasing the pump width, at least for this given system.It is important to mention that the resultant entangled state is always close to |⟩ = |HV⟩ − |VH⟩ state.

Summary and conclusion
In summary, we have developed a theoretical framework based on the classical Green's function of the system, to comprehensively describe spontaneous parametric downconversion (SPDC) in nonlinear thin films.This vectorial formalism is capable of treating subwavelength slab thicknesses, accounts for Fabry-Pérot effects, and can address absorption in the slab.Additionally, we have introduced a simplified formulation for numerical investigation of farfield properties of photon pairs detected within one of the media surrounding the slab.Using a zinc-blende ⟨100⟩ GaAs crystal as a specific example, deposited on a SiO 2 substrate with air as the detecting medium, we have computed probability rates for both degenerate and non-degenerate emissions.Through analytical calculations and numerical simulations, we examined and explained the impact of Fabry-Pérot interferences, characterized by oscillations in emitted angles and probability distributions.
Furthermore, we explored the potential of this system to produce polarization-entangled states.Employing a tomographic approach, we demonstrated that a GaAs slab can generate maximally polarization-entangled states of the form |⟩ = |HV⟩ − |VH⟩, within a -symmetric configuration by analyzing its density matrix ρ.Our findings revealed that the entanglement achieved remains highly resilient against variations in wavelengths and thickness, even in non-degenerate scenarios.However, in the case of non-degenerate photons, achieving maximal entanglement entails a trade-off with the efficiency of the process, where the latter is affected by the transversal phasematching condition.This suggests a nuanced interplay between entanglement quality and detection rates in such systems, which could potentially be manipulated by using spatially structured pump fields to reach maximum entanglement together with high detection rates.
that describes the amplitude propagation of a field with angular spectrum U(k x , k y )∕k z at the initial position r ′ = (x ′ , y ′ , z ′ = 0) to the final position r = (x, y, z).With wave vector component in the y .To do a far-field approximation of this propagation integral, we first use the convolution theorem to write Eq. ( 15) as a convolution of the form where * denotes the convolution operation, U(k x , k y ) =  [u(x, y)] and we have used the Weyl expansion of an outgoing spherical wave as a superposition of plane waves, Equation ( 18) is interpreted as the beam propagation of an initial field u(x ′ , y ′ ) to a distance z with an outgoing spherical wave as the response function h = e ikr ∕r.In the spherical coordinates representation with x = rsincos, y = rsinsin, and z = rcos, R takes the form By taking the condition that r ≫ r ′ the last to terms are approximated to be zero, and using a first order binomial expansion of R, we arrive at Now, when introducing R into (18), we notice that in the response function the phase factor varies rapidly in comparison to its denominator.Therefore we keep the 1st order expansion for the phase factor while using a 0th order as a good enough approximation for the denominator.Thus, we get the far-field-approximated expression for beam propagation, ) This equation tells us that, under the far-field approximation, the amplitude propagation to the far-field of an initial beam field defined at plane z = 0, with angular spectrum U(k x , k y )∕k z is proportional to U(k x , k y ) with k x = k sincos and k y = k sinsin.This approximation, which is similar to the Fraunhofer approximation for the diffraction of classical beams [40], also considers nonparaxial propagation, meaning that it also describes propagation at angles up to 90 • .However, our goal here is to extend this concept to the two-photon joint-detection probability.To this aim, we start with Eq. ( 5), and notice that the integral inside the absolute value can be written as where k z,1 is the z component of the wave vector in the medium of detection (medium 1 in this case).The factors e ik z,1s z s and e ik z,1i z i come from the definition of the GF in Eq. (10), which allows us to write (q s , q i , z s , z i ,  s ,  p −  s ) = (q s , q i ,  s ,  p −  s )e ik z,1s z s e ik z,1i z i . ( The integral in Eq. ( 23) is a two-photon analogue of the propagation integral in Eq. ( 15) and can be interpreted as the propagation of the probability amplitude of the generated signal and idler photon pair with joint angular probability (JAP) amplitude (q s , q i ,  s ,  p −  s ).Therefore, we can perform a two-photon-like far-field approximation where we identify  in Eq. ( 23) as U∕k z in Eq. ( 15), and by using the identity in (22) twice (for signal and idler propagation), we obtain Here, q s = k x,s x + k y,s ŷ, where k x,s = k 1s sin s cos s and k y,s = k 1s sin s sin s , refer to the transverse k-vector of the signal photon detected in the far-field at propagation angles { s ,  s }, and similarly for the idler photon.After inserting (25) in ( 5) we arrive at the expression in (7) which describes a far-field-approximated joint-probability rate of detecting coincidences in the far-field as a function of the detection angles { s ,  i ,  s ,  i }.

Appendix B: Green's function of a dielectric slab
The GFs in Fourier domain for the multilayered system depicted in Figure 1, where a field is generated in medium 2 and propagates outside the slab can be constructed as 21 (q, z ′ , ) ŝŝ and where we have also added the GF connecting medium 2 and 3 for completeness.
Using the definitions for the sand p-polarized fields in Eqs. ( 8) and ( 9) we can compute the dyadics ŝŝ, p1+ p2± and p3− p2± as ŝŝ = The quantities T (s) and T ( p) correspond to the generalized transmission coefficients for the sand p-polarized fields respectively that takes into account all internal reflections in the slab.These coefficients are based on the well-known Fresnel reflection and transmission coefficients from medium i to medium j, r Thus, T (s, p) can be calculated as [31], [39] and where the corresponding Fresnel coefficient for the sand p-polarized vectors must be used.As explained in Section 3, the s-polarized field is direction independent, and there is no distinction in directionality for T (s) .However, the p-polarized field is direction dependant and the coefficient T ( p) needs to be treated differently for upward and downward propagating waves.For the ppolarized field one can write, T ( p) 21 (q, z ′ ) = T ( p+)  21 (q, z ′ ) + T ( p−) 21 (q, z ′ ), (35) which splits the transmission into an upward field generated in medium 2 that transmits to medium 1, namely T ( This distinction becomes important when we want to construct the GF of our system and it was not considered in Ref. [39]. For completeness, the same analysis is applied to the GF in Eq. ( 27) for waves generated in medium 2 and transmitting to medium 3. Here, a downward p-wave generated in medium 2 that also transmits to medium 3 after multiple reflections.

Appendix C: Pump field with multiple reflections C.1 Gaussian angular spectrum
To test the contribution of the longitudinal polarization component of the pump beam to the photon-pair generation process we take the Gaussian angular spectrum defined in (13) and introduce it to the pump field defined in (12).Using  p = 500 nm as the pump wavelength, we can observe in Figure 6 the intensity of the spectrum, xand z-component of the field in Fourier domain for two spectral pump widths .For a width  = 6.6 × 10 5 m −1 Figure 6(a) shows that the spectrum of spatial frequencies are much smaller than the wave-vector (q ≪ k p = 2∕ p ), this is the so-called nonparaxial regime (weakly focused beam in position space).Here, U p (q) corresponds to the angular spectrum of a Gaussian beam of width W = 3 μm defined as u p (x, y) ∝ e −(x 2 + y 2 )∕W 2 .In this regime, the contribution of the longitudinal component of the pump field is very small compared to the transversal component, as it can be seen in ] .In this regime, the individual plane waves in the superposition (3) can be approximated to be only x-polarized, due to the small contribution of the z-component of the pump electric field.For comparison, we show in Figure 6(b) the same quantities for a width  = 2 × 10 7 m −1 .This is a nonparaxial scenario where the angular spectrum saturates to values close to ∼ k p .The Here, the angular spectrum corresponds to a tightly focused beam and the intensity of the z-component of the pump is comparable to the x-component.

C.2 Pump field inside the slab
Right before the interface at z = −a we take an x-polarized electric pump with field components polarized in the direction of propagation, E p (q, z = −a) = (E p,x , 0, E p,z ), where its longitudinal component E p,z is calculated from E p,x using the transverse condition k p ⋅ E p = 0.This results in We express the pump field at the interface in terms of the sand p-polarized fields as E p (q, z = −a) = E p,x (q, z = −a)x + E p,z (q, z = −a)ẑ = E (s)  p (q, z = −a)ŝ + E ( p) p (q, z = −a)p 3+ , (41) where we used the definitions of ŝ and p in ( 8) and ( 9) for medium 3.After solving we find that p (q, z = −a) = k z,3 p k y qk x E p,z (q, z = −a), where Eq. ( 40) still holds.Notice that we have a singularity at k x = k y = 0, which describes a plane wave normal to the interface.For this special point, we take the field to be only xpolarized, that is E P (q = 0, z = −a) = E P,x (q = 0, z = −a)x.
To describe the transmission of the field to the slab, one can calculate a corresponding coefficient associated to the sand p-polarized pump field coming from the outside (medium 3) to the inside of the slab (medium 2), which goes along the lines to the generalized Fresnel transmission coefficients (Eq.( 33)).The calculation of the coefficients is done in the Appendix C.3 and takes the form

Figure 1 :
Figure 1: Multilayered geometry where medium 2 is a nonlinear medium, in which SPDC happens due to a pump field impinging from medium 3. The joint signal and idler field can radiate to outside the slab to medium 1 and medium 3, after experiencing multiple reflections inside medium 2.

Figure 2 :
Figure 2: Far-field joint detection probability rate in the (a) -symmetric configuration, which is considered for all the results shown in this figure.(b) Probability rate for an ultra-thin film of GaAs of thickness a = 0.01 p in the -symmetric configuration.(c, d) Slice cut through far-field radiation pattern in the xz plane for different thicknesses (below and above a =  p ) of the slab with normalized intensities.(e) Angle of maximum joint emission as a function of the thickness of the slab.The four values plotted in (d) are marked by the stars, corresponding to the peak (blue star at ∼58.5 • ), two in the middle (green and purple stars at ∼45 • ) and the valley (red star at ∼27.5 • ).(f) Intensity of joint detection at  s = 45 • which exhibit FP oscillations.(g) Angle uncertainty Δ of joint emission as a function of the thickness of the slab.The values in (d) are also in (f) and (g).The period of oscillation Δa ∼ 0.29 p .
12 • , matching Figure 3(a).For r = 1 the maximum coincidence rate lie along  s =  i , justifying the previously used symmetric configuration.Additionally, in Figure 3(a), we observe that as the degeneracy factor r decreases, the uncertainty in idler angles increases, meaning weaker correlations in emitted photons.Due to transverse phase matching, most pairs are generated under the condition  s =  i + 180 • , which suggests visualizing the joint-detection rate as a function of the angles { s ,  i }, in a particular plane defined by .In

Figure 3 :
Figure 3: Coincidence detection rate in the non-degenerate case for (a) different degeneracy factors r for a fixed signal detector position at  s = 45 • ,  s = 180 • and varying idler detector positions.The maximum angle of detection in the plot are  i (r = 0.8) ∼ 61.2 • ,  i (r = 1) ∼ 45 • ,  i (r = 1.2) ∼ 36.3 • and  i (r = 1.5) ∼ 28.2 • .(b) Joint detection probability as a function of  s and  i in the xz-plane (with fixed  s = 0 and  i = 180 • ) for r = 1.5.The thickness is a = 0.01 p in both graphs.

Figure 3 (
Figure3(b) we plot  along the xz-plane, where  i = 0 and  s = 180 • , as a function of  s and  i for r = 1.5 and a = 0.01 p .Here we can see that a better correlation in the non-degenerate case can be achieved by detecting at small angles, however this comes at expenses of the efficiency (b), these angles are  s = 55.35 • and  i = 33.3• .

Figure 4 :
Figure 4: Quantum polarization tomography.(a) Real and imaginary parts of the elements of the density matrix, ρ, for the polarization state of a degenerate (r = 1) photon pair, generated from an ultra-thin GaAs slab of thickness a = 0.01 p .The detectors are in a -symmetric configuration with  s = 45 • and  s = 0. (b) Schmidt number and joint-coincidence rate as a function of the signal wavelength (bottom axis) and idler wavelength (top axis) for the ultra-thin GaAs slab in the same -symmetric configuration, for two thicknesses a = 0.01 p and a = 0.1 p .
(b)).This demonstrates GaAs's ability to produce polarizationentangled Bell states of the form |⟩ = |HV⟩ − |VH⟩ across a wide range of wavelengths, maintaining near-perfect entanglement at -symmetric detection with minor degradation at non-degenerate wavelengths.It is interesting to explore the preservation of broad entanglement obtained in the -symmetric setup across various emission angles, especially in non-degenerate scenarios where the -symmetric configuration is not dominant.Figure 5 displays the Schmidt number (K) variation with { i ,  i }, with fixed signal detection at  s = 45 • and  s = 0, for both degenerate (Figure 5(a) and (c)) and nondegenerate cases with r = 1.5 (Figure 5(b) and (d))

Figure 5 :
Figure 5: Schmidt number as a function of the idler angles { i ,  i } with fixed  s = 45 • and  s = 0 for (a) a = 0.01 p and r = 1, (b) a = 1 p and r = 1.5, (c) a = 1 p and r = 1, (d) a = 1 p and r = 1.5.The angle of maximum entanglement and the angle of maximum detection rate are marked by a red " * " and a green "x", respectively.

Research funding :
The authors acknowledge funding by the German Ministry of Education and Research (13N14877, 13N16108), the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Collaborative Research Center NOA (CRC 1375, project number 398816777), the German Federal Ministry of Education and Research BMBF through project QOMPLEX (13N15985), the Thuringian State Government through project Quantum Hub (project number 2021 FGI 0043), the Nexus program of the Carl-Zeiss-Stiftung (project MetaNN), by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project-Nr.505897284 and Projekt-Nr.512648189 and the Open Access Publication Fund of the Thüringer Universitäts-und Landesbibliothek Jena.I(x, y, z) = 1 (2) 2 ∬ dk x dk y U(k x , k y ) e ik z z k z e i(k x x+k y y) ,(15)

Figure 6
Figure 6(a) where the maximum Max [ E 2 P,z (q) ]

Table 1 :
Purity and concurrence of the two-photon state for several collection numerical apertures (NAs), with a fixed pump width W = 3 μm.Signal and idler are frequency degenerate and slab thickness is a = 0.01 p .

Table 2 :
Purity and concurrence of the two-photon state for several pump widths W, with a fixed collection numerical aperture NA = 0.4.