# Entangled Photon-pair Generation in Nonlinear Thin-films

Elkin A. Santos<EMAIL_ADDRESS>Institute of Applied Physics, Abbe
Center of Photonics, Friedrich Schiller University Jena, Albert-Einstein-Str.
15, 07745 Jena, Germany  Maximilian A. Weissflog  Institute of Applied
Physics, Abbe Center of Photonics, Friedrich Schiller University Jena, Albert-
Einstein-Str. 15, 07745 Jena, Germany  Max Planck School of Photonics, Hans-
Knöll-Straße 1, 07745 Jena, Germany  Thomas Pertsch  Institute of Applied
Physics, Abbe Center of Photonics, Friedrich Schiller University Jena, Albert-
Einstein-Str. 15, 07745 Jena, Germany  Fraunhofer Institute for Applied Optics
and Precision Engineering IOF, Albert-Einstein-Str. 7, 07745 Jena, Germany
Frank Setzpfandt  Institute of Applied Physics, Abbe Center of Photonics,
Friedrich Schiller University Jena, Albert-Einstein-Str. 15, 07745 Jena,
Germany  Fraunhofer Institute for Applied Optics and Precision Engineering
IOF, Albert-Einstein-Str. 7, 07745 Jena, Germany  Sina Saravi  Institute of
Applied Physics, Abbe Center of Photonics, Friedrich Schiller University Jena,
Albert-Einstein-Str. 15, 07745 Jena, Germany

###### Abstract

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érot 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érot 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.

††preprint: APS/123-QED

##  I Introduction

The optical process of spontaneous parametric down-conversion (SPDC) in
$\chi^{(2)}$-materials spdc , in which a pump ($p$) photon of wavelength
$\lambda_{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 Flamini2019 ,
quantum cryptography Bennet ; pirandola2020 , quantum imaging FrankRev ;
moreau  and sensing clark , and testing fundamental quantum effects Horodecki
. 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 wang2021 . Entangled photon pairs have been
generated in single nonlinear bulk crystals Kwiat1995 , crossed bulk crystals
Kwiat1999 , or bulk crystals in a Sagnac loop Hentschel . More recently,
several experiments of photon-pair generation in much smaller nonlinear
systems such as single GaAs nanowires saerens2023 , dielectric nanoantennas
Marino2019 , and metasurfaces with subwavelength thicknesses santiago2021 ;
Santiago2022 ; sukho2022 ; Son2023  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) Okoth2019 ; Chek_2021 ;
Chek_2022 , as well as in van der Waals guo2023  or transition metal
dichalcogenide (TMD) crystals weissflog2023tunable ;
trovatelloQuasiphasematchedDownconversionPeriodically2023 . 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 vega2022 ; santos2022 .

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 Krasnok . 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
degrees.

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 poddubny2016 . 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
welsch1996 ; welsch1998 . 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 Sipe . 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 to keep 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. While previous efforts pevrina2011 ; pevrina2014  have sought to
elucidate SPDC in nonlinear layered structures through the quantization of the
electromagnetic field in free space mandel ; howell , our approach utilizing
the GF of the system offers a more intuitive, direct vectorial method, which
facilitates the exploration of entanglement properties. Additionally, the GF
inherently accounts for losses and near-field interactions, which are crucial
considerations in such systems.

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 $\bar{4}3$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 weissflog2024 ; Chek_2022 . 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 Sec. (II) 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 Sec.
(III), 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
Sec. (IV), the classical pump beam is treated. In Sec. (V), we numerically
compute the far-field radiation for a nonlinear source made of GaAs deposited
on a silicon dioxide (SiO2) substrate, and study its far-field properties. In
Sec. (VI), 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 Sec. (VII).

##  II Far-Field Joint Spatial Probability in Fourier Domain

Consider a general nonlinear thin-film photon-pair source, as depicted in Fig.
1(a). It comprises three layers with relative permittivities $\epsilon_{1}$,
$\epsilon_{2}$, and $\epsilon_{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
$\mathbf{E}_{p}(\mathbf{r},t)=\mathbf{E}_{p}(\mathbf{r})\exp{-\mathrm{i}\omega_{p}t}+\text{c.c.}$.
The far-field coincidence detection rate for a pair of signal and idler
photons of different frequencies as a function of the spatial coordinates
follows poddubny2016 ; weissflog2024 :

$\displaystyle\mathcal{R}(\mathbf{r}_{s},\mathbf{r}_{i},\mathbf{e}_{s},\mathbf{e}_{i},\omega_{s},\omega_{p}-\omega_{s})\propto
n_{i}n_{s}r_{i}^{2}r_{s}^{2}\frac{(\omega_{p}-\omega_{s})^{3}\omega_{s}^{3}}{c^{6}}$
$\displaystyle\times\bigg{|}\sum_{\alpha,\beta,\gamma}\sum_{\sigma_{s},\sigma_{i}}e_{s,\sigma_{s}}e_{i,\sigma_{i}}\,\int_{V}\differential\mathbf{r}\,\chi^{(2)}_{\alpha\beta\gamma}(\mathbf{r})$
$\displaystyle\times
E_{p,\gamma}(\mathbf{r})\,G_{\sigma_{s}\alpha}(\mathbf{r}_{s},\mathbf{r},\omega_{s})\,G_{\sigma_{i}\beta}(\mathbf{r}_{i},\mathbf{r},\omega_{p}-\omega_{s})\bigg{|}^{2}\,.$
(1)

The detection rate, $\mathcal{R}\equiv\frac{\differential^{4}N}{\differential
t\differential\Omega_{s}\differential\Omega_{i}\differential\omega_{s}}$, is
the number of photon pairs given per units of solid angles
$\differential\Omega_{i}$ (for the idler photon) and $\differential\Omega_{s}$
(for the signal photon), and per unit of signal-photon angular frequency
$\omega_{s}$. The idler angular frequency is fixed by conservation of energy
to $\omega_{i}=\omega_{p}-\omega_{s}$. The detection rate is also a function
of the idler and signal photons’ detection polarization, which are along the
unitary $\mathbf{e}_{s}$ and $\mathbf{e}_{i}$ directions. The $\alpha$,
$\beta$ $\gamma$, $\sigma_{s}$ and $\sigma_{i}$ indices run over the $x$, $y$,
and $z$ cartesian coordinates. $G_{ij}(\mathbf{r,r^{\prime}},\omega)$ are the
tensor-components of the electric GF,
$\chi^{(2)}_{\alpha\beta\gamma}(\mathbf{r})$ are the components of the second-
order nonlinear tensor and the integral runs over the volume of nonlinearity
$V$.

Figure 1: (a) 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.

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 can be represented as an inverse Fourier transform of the form
lau2007

$\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{G}}(\mathbf{r},\mathbf{r}^{\prime},\omega)=\frac{1}{(2\pi)^{2}}\int_{-\infty}^{\infty}\differential\mathbf{q}\;\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{g}}(\mathbf{q},z,z^{\prime},\omega)\,e^{\mathrm{i}\mathbf{q}\cdot(\mathbf{r}_{\perp}-\mathbf{r}^{\prime}_{\perp})}\,,$
(2)

where $\mathbf{q}=k_{x}\hat{x}+k_{y}\hat{y}$ and
$\mathbf{r}_{\perp}=x\hat{x}+y\hat{y}$ are the two-dimensional real wave-
vector and spatial vector in the $x$-$y$ plane, respectively. We define
$\int_{-\infty}^{\infty}\differential\mathbf{q}\equiv\int_{-\infty}^{\infty}\differential
k_{x}\int_{-\infty}^{\infty}\differential k_{y}$. Similarly, for a classical
pump field treated in the undepleted pump approximation we have,

$\mathbf{E}_{p}(\mathbf{r})_{=}\frac{1}{(2\pi)^{2}}\int_{-\infty}^{\infty}\differential\mathbf{q}\;\mathbf{E}_{p}(\mathbf{q},z)\;e^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}_{\perp}}\,.$
(3)

Then, by introducing Eqs. (2) and (3) into (II), we can separate the spatial
integral into $\differential\mathbf{r}\rightarrow\differential
z\,\differential\mathbf{r}_{\perp}$. After rearranging terms, we can solve for
$\differential\mathbf{r}_{\perp}$ the integral

$\int_{-\infty}^{\infty}\differential\mathbf{r}_{\perp}\;e^{\mathrm{i}(\mathbf{q}_{p}-\mathbf{q}_{s}-\mathbf{q}_{i})\cdot\mathbf{r}_{\perp}}=(2\pi)^{2}\delta(\mathbf{q}_{p}-\mathbf{q}_{s}-\mathbf{q}_{i})\,,$
(4)

which results into the joint detection rate

$\displaystyle\mathcal{R}(\mathbf{r}_{s},\mathbf{r}_{i},\mathbf{e}_{s},\mathbf{e}_{i},\omega_{s},\omega_{p}-\omega_{s})\propto
n_{i}n_{s}r_{i}^{2}r_{s}^{2}\frac{(\omega_{p}-\omega_{s})^{3}\omega_{s}^{3}}{c^{6}}$
$\displaystyle\times\bigg{|}\frac{1}{(2\pi)^{4}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\differential\mathbf{q}_{s}\,\differential\mathbf{q}_{i}\;\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},z_{s},z_{i},\omega_{s},\omega_{p}-\omega_{s})\;$
(5) $\displaystyle\times
e^{\mathrm{i}(\mathbf{q}_{s}\cdot\mathbf{r}_{s\perp}+\mathbf{q}_{i}\cdot\mathbf{r}_{i\perp})}\bigg{|}^{2}\,,$

with

$\displaystyle\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},z_{s},z_{i},\omega_{s},\omega_{p}-\omega_{s})=\sum_{\alpha,\beta,\gamma}\sum_{\sigma_{s},\sigma_{i}}e_{s,\sigma_{s}}e_{i,\sigma_{i}}$
$\displaystyle\times\int_{-a}^{0}\differential
z\,\chi^{(2)}_{\alpha\beta\gamma}(z)\,E_{p,\gamma}(\mathbf{q}_{s}+\mathbf{q}_{i},z)$
$\displaystyle\times
g_{\sigma_{s}\alpha}(\mathbf{q}_{s},z_{s},z,\omega_{s})\,g_{\sigma_{i}\beta}(\mathbf{q}_{i},z_{i},z,\omega_{p}-\omega_{s})\,.$
(6)

Here, equation (II) is understood as a joint angular probability (JAP).
Equation (II) 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, $\mathcal{R}$ in the form of (II) 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 born . We aim to detect both photon pairs in the far field in
medium 1 with refractive index $n_{1}$. Here, $k_{1}=\frac{\omega}{c}n_{1}$ is
the wave number of plane waves in medium 1, and
$k_{z,1}=\sqrt{k_{1}^{2}-k_{x}^{2}-k_{y}^{2}}$ denotes the $z$-component of
its wave vector $\mathbf{k}_{1}=\mathbf{q}+k_{z,1}\hat{z}$. Notably, $k_{x}$
and $k_{y}$ represent 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.,
$\mathbf{r}_{s}\gg\mathbf{r}$ and $\mathbf{r}_{i}\gg\mathbf{r}$. Additionally,
we assume a finite-sized 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. (II) is simplified to
(see appendix A)

$\displaystyle\mathcal{R}^{(\text{far})}(\theta_{s},\theta_{i},\varphi_{s},\varphi_{i},\mathbf{e}_{s},\mathbf{e}_{i},\omega_{s},\omega_{p}-\omega_{s})\propto$
$\displaystyle
n_{i}n_{s}\frac{(\omega_{p}-\omega_{s})^{3}\omega_{s}^{3}}{c^{6}}\,\bigg{|}\tilde{R}(\mathbf{q}_{s},\mathbf{q}_{i},\omega_{s},\omega_{p}-\omega_{s})\;k_{z,1s}k_{z,1i}\bigg{|}^{2},$
(7)

with $k_{x,s}=k_{1s}\sin\theta_{s}\cos\varphi_{s}$ and
$k_{y,s}=k_{1s}\sin\theta_{s}\sin\varphi_{s}$, and similarly for the idler
photon. Thus, Eq. (II) describes the joint detection rate in the far-field as
a function of the detection angles $\\{\theta_{s}$, $\theta_{i}$,
$\varphi_{s}$, $\varphi_{i}\\}$ in spherical coordinates. Note that $k_{z,1s}$
and $k_{z,1i}$ in Eq. (II) correspond to the detection medium 1. If the
detection medium changes, these factors should be adjusted accordingly (see
appendix A). In the far-field, the joint detection rate is proportional to its
Fourier transform multiplied by the signal and idler $z$-components of the
wave vector, $k_{z,1s}$ and $k_{z,1i}$. In the paraxial regime ($q\ll k_{1}$
with $q=\absolutevalue{\mathbf{q}}$), $k_{z,1s}$ and $k_{z,1i}$ tend to
$k_{1s}$ and $k_{1i}$, resulting in the joint detection rate being directly
proportional to its Fourier transform. This calculation extends the concept of
classical Fraunhofer approximation born  to two-photon joint-detection. 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 were also shown in studying the fundamental
resolution limit of quantum imaging schemes vega2022 . Equations (II), (II),
and (II) are the main results of this section, to be numerically solved in
section (V).

##  III 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.
(II), we employ a GF suitable for the multilayered scheme depicted in Fig.
(1). In this scenario, the GFs for signal and idler,
$\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{\mathbf{G}}}(\mathbf{r},\mathbf{r}_{s})$
and
$\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{\mathbf{G}}}(\mathbf{r},\mathbf{r}_{i})$
connect a point $\mathbf{r}$ within the $\chi^{(2)}$-nonlinearity region
(medium 2, see Eq. (II)) to the detection positions $\mathbf{r}_{s}$ and
$\mathbf{r}_{i}$ in medium 1.

If a plane wave of the form
$\mathbf{E}(\mathbf{r})=\mathbf{E}e^{\mathrm{i}\mathbf{k}\cdot\mathbf{r}}$,
with $\sqrt{\mathbf{k}\cdot\mathbf{k}}=k=\frac{\omega}{c}n$ satisfies
Maxwell’s equation in a homogeneous medium of refractive index $n$, two
possible wave vectors $\mathbf{k}$ for a given $\mathbf{q}$ are obtained,
representing upward and downward propagation. These wave vectors take the form
$\mathbf{k}_{+}=\mathbf{q}+k_{z}\hat{z}$ (upward) and
$\mathbf{k}_{-}=\mathbf{q}-k_{z}\hat{z}$ (downward). When $n$ is real and
$q>n\omega/c$, $k_{z}$ becomes imaginary, resulting in an evanescent wave.
Additionally, at a given $\mathbf{q}$, a solution can be either $s$-polarized
(transverse electric) or $p$-polarized (transverse magnetic), defined by unit
polarization vectors $\hat{s}$ and $\hat{p}$ as

$\displaystyle\hat{s}=$
$\displaystyle\frac{\hat{z}\crossproduct\mathbf{k}_{\pm}}{\absolutevalue{\hat{z}\crossproduct\mathbf{k}_{\pm}}}=\frac{1}{q}(-k_{y}\hat{x}+k_{x}\hat{y})\;,$
(8) $\displaystyle\hat{p}_{\pm}=$
$\displaystyle\frac{\mathbf{k}_{\pm}}{\absolutevalue{\mathbf{k}_{\pm}}}\crossproduct\hat{s}=\frac{1}{k}\bigg{(}\mp\frac{k_{z}k_{x}}{q}\hat{x}\mp\frac{k_{z}k_{y}}{q}\hat{y}+q\hat{z}\bigg{)}\,.$
(9)

These vectors are normalized according to $\hat{s}\cdot\hat{s}=1$ and
$\hat{p}_{\pm}\cdot\hat{p}_{\pm}=1$. Note that $\hat{s}$ is real, rests on the
$x$-$y$ plane, and is medium independent, while $\hat{p}$ can be complex and
is medium dependent. Also, $\hat{s}$ does not depend on the direction of
$\mathbf{k}$, while $\hat{p}$ changes for upward- or downward-generated
fields.

Using the $s$\- $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 Sipe ; lau2007

$\displaystyle\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{g}}_{21}(\mathbf{q},z,z^{\prime},\omega)=\frac{-\mathrm{i}}{2k_{z,2}}\bigg{(}T_{21}^{(s)}(q,z^{\prime},\omega)\;\hat{s}\hat{s}+$
$\displaystyle
T_{21}^{(p+)}(q,z^{\prime},\omega)\;\hat{p}_{1+}\hat{p}_{2+}+T_{21}^{(p-)}(q,z^{\prime},\omega)\;\hat{p}_{1+}\hat{p}_{2-}\bigg{)}e^{\mathrm{i}k_{z,1}z}\,,$
(10)

where $T_{21}(q,z^{\prime},\omega)$ is the generalized transmission
coefficient of the $s$\- and $p$-polarized fields generated from the source
plane at $z^{\prime}\in\\{-a,0\\}$ and transmitted to medium 1 (see appendix
(B) for analytical expression). The unit vector $\hat{s}$ is parallel to the
interface and $\hat{p}$ is perpendicular to the wave vector $\mathbf{k}$ and
$\hat{s}$. The subscript numbers in $k_{z}$ indicate the corresponding medium.

The GF in Eq. (III) describes from right to left (apart from the phase factor
$e^{\mathrm{i}k_{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}$.

##  IV Pump Field with Multiple Reflections

The pump field undergoes diffraction and multiple reflections within the slab.
Illustrated in Fig. (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
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

$\displaystyle\mathbf{E}_{p}(\mathbf{q},$ $\displaystyle z=-a)=$
$\displaystyle
E_{p,x}(\mathbf{q},z=-a)\hat{x}-E_{p,x}(\mathbf{q},z=-a)\frac{k_{x}}{k_{z,3p}}\hat{z}\,,$
(11)

where $k_{z,3p}$ 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 $\mathbf{k}_{p}\cdot\mathbf{E}_{p}=0$ (see Appendix (C.2)).

Taking
$\mathbf{E}_{p}(\mathbf{q},z=-a)=U_{p}(\mathbf{q},z=-a)\mathbf{\hat{e}}(\mathbf{q})$
one finds
$\mathbf{\hat{e}}(\mathbf{q})={k_{z,3p}}/{\sqrt{k_{3p}^{2}-k_{y}^{2}}}\;\hat{x}-{k_{x}}/{\sqrt{k_{3p}^{2}-k_{y}^{2}}}\;\hat{z}$,
which allows us to express the pump field in the form

$\displaystyle\mathbf{E}_{p}$ $\displaystyle(\mathbf{q},z=-a)=$ $\displaystyle
U_{p}(\mathbf{q},z=-a)\left[\frac{k_{z,3p}}{\sqrt{k_{3p}^{2}-k_{y}^{2}}}\hat{x}-\frac{k_{x}}{\sqrt{k_{3p}^{2}-k_{y}^{2}}}\hat{z}\right]\,,$
(12)

where we take $\mathbf{\hat{e}}=\hat{x}$ for $k_{x}=k_{z}=0$.
$U_{p}(\mathbf{q},z=-a)$ is the angular spectrum of the pump field at the
interface. We assume it follows a Gaussian distribution,

$U_{p}(\mathbf{q},z=-a)=\begin{cases}Ae^{-\mathbf{q}^{2}/w^{2}}\,,&\quad
q^{2}\leq k_{3p}^{2}\\\ 0\,,&\quad\text{otherwise,}\end{cases}$ (13)

where $w$ is the width of the Gaussian angular spectrum. The condition
$q^{2}\leq k_{3p}^{2}$ 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 $w$ 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 $s$\- and
$p$-polarized fields. After transmission into the slab, the $z$-dependent pump
beam can be calculated as (see Appendix (C.2) for details)

$\displaystyle\mathbf{E}_{p}(\mathbf{q},z)=Q^{(s)}(q,z)E^{(s)}_{p}(\mathbf{q},z=-a)\hat{s}$
$\displaystyle+E^{(p)}_{p}(\mathbf{q},z=-a)\left(Q^{(p+)}(q,z)\hat{p}_{2+}+Q^{(p-)}(q,z)\hat{p}_{2-}\right)\,,$
(14)

where $Q^{(s)}(q,z)$ is the transmission coefficient from medium 3 to medium
2, while $E^{(s)}_{p}(\mathbf{q},z=-a)$ and $E^{(p)}_{p}(\mathbf{q},z=-a)$ are
$s$\- and $p$-components of the pump before the slab. Unit vectors $\hat{s}$,
$\hat{p}_{2+}$, and $\hat{p}_{2-}$ are defined in Eqs. (8, 9) for medium 2.

##  V Far-Field Joint-Radiation Properties

###  V.1 Degenerate photon-pair detection

Now that we have the model, let us consider a Zinc-Blende GaAs crystal (medium
2) with thickness $a$ and orientation $\left\langle 100\right\rangle$. The
crystal axes are aligned with the laboratory coordinates $\\{x,y,z\\}$. The
GaAs is on a SiO2 substrate (medium 3, assumed semi-infinite), with air as
medium 1 for detection. The pump from medium 3 is described by Eq. (IV) with
$\lambda_{p}=500$ nm and Gaussian angular spectrum (Eq. (13)) with
$w=6.6\times 10^{5}$ m-1 (corresponding to a beam waist of $W=3\,\mu$m). At
degenerate wavelengths ($\lambda_{s}=\lambda_{i}=1\,\mu$m), we have
$\epsilon_{2s}=\epsilon_{2i}=12.06$ in GaAs, while the pump permittivity at
$\lambda_{p}=500$ nm in GaAs is $\epsilon_{2p}=17.63+3.83\mathrm{i}$
papatry2021 . In SiO2 substrate, $\epsilon_{3p}=2.14$ and
$\epsilon_{2s}=\epsilon_{2i}=2.10$ Malitson65 , while in the detection medium,
$\epsilon_{1p}=\epsilon_{1s}=\epsilon_{1i}=1$. The non-vanishing elements of
the second-order susceptibility tensor of GaAs are
$\chi^{(2)}_{yzx}=\chi^{(2)}_{zyx}=\chi^{(2)}_{xzy}=\chi^{(2)}_{zxy}=\chi^{(2)}_{xyz}=\chi^{(2)}_{yxz}=\chi^{(2)}_{0}$.

We start by computing the coincidence detection rate
$\mathcal{R}^{(\text{far})}$ in Eq. (II) for unpolarized joint-detection. This
involves summing over the three polarization directions of each signal and
idler detector ($\mathbf{e}=\hat{x}$, $\hat{y}$, $\hat{z}$), resulting in nine
combinations. Due to the relatively narrow angular range of our pump beam,
opposite transverse wave vectors ($\mathbf{q}_{s}\sim-\mathbf{q}_{i}$)
significantly contribute to the generation, as consequence of the transverse
phase-matching condition ($\mathbf{q}_{p}=\mathbf{q}_{s}+\mathbf{q}_{i}$). For
degenerate photon-pairs, this motivates investigating
$\mathcal{R}^{(\text{far})}$ in terms of propagation angles
$\\{\theta_{s},\varphi_{s}\\}$ for the signal photon, where the idler is
detected symmetrically at
$\\{\theta_{i}=\theta_{s},\varphi_{i}=\pi+\varphi_{s}\\}$. This configuration
is shown in Fig. 2(a) and is called $\varphi$-symmetric throughout the
manuscript. The $\varphi$-symmetric setup has already been proven crucial for
entangled state generation in nanoresonators weissflog2024  and will also show
to be essential for thin films in this work.

The coincidence rate in the $\varphi$-symmetric configuration is shown in Fig.
2(b) for an ultra-thin source of thickness $a=0.01\lambda_{p}$, emphasizing
again that both photons are detected in medium 1. Here, the $\chi^{(2)}_{yzx}$
and $\chi^{(2)}_{zyx}$ dominate the pair generation process due to weak $y$\-
and $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. (IV)). Hence,
all components of the GaAs nonlinear tensor participate in the generation
process. Yet, the dominance of $\chi^{(2)}_{yzx}$ and $\chi^{(2)}_{zyx}$
components implies that signal and idler fields within the slab are mainly
generated by $y$\- and $z$-polarized point-like nonlinear sources, resulting
in no radiation in the $y$-$z$ plane, as observed in Fig. 2(b). Additionally,
emission at large angles is significantly reduced, despite the GaAs nonlinear
tensor favoring stronger emission near $\theta_{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 $\theta_{s}=90$° due to total
internal reflection experienced by photons inside the slab.

Figure 2:  (a) $\varphi$-symmetric configuration for joint detection, which is
considered for all the results shown in this figure. (b) far-field joint
detection probability for an ultra-thin film of GaAs of thickness
$a=0.01\lambda_{p}$ in the $\varphi$-symmetric configuration. (c,d) slice cut
through far-field radiation pattern in the $xz$ plane for different
thicknesses (below and above $a=\lambda_{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 $\sim 58.5$°), two in the middle
(green and purple stars at $\sim 45$°) and the valley (red star at $\sim
27.5$°). (f) Intensity of joint detection at $\theta_{s}=45$° which exhibit FP
oscillations. (g) Angle uncertainty $\Delta\theta$ 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 $\Delta a\sim 0.29\lambda_{p}$.

In Fig. 2(c), we display cuts along the $xz$-plane for various subwavelength
thicknesses of the slab, $a=0.01\lambda_{p}$, $a=0.1\lambda_{p}$,
$a=0.5\lambda_{p}$ and $a=0.8\lambda_{p}$. Fig. 2(d) illustrates $xz$-plane
cuts for larger slab thicknesses, $a=3.61\lambda_{p}$, $a=3.67\lambda_{p}$,
$a=3.72\lambda_{p}$ and $a=3.82\lambda_{p}$, which are special points for
understanding the Fabry-Pérot dynamics, marked in Figs. 2(e-g). All graphs are
normalized to one for better observation of angular dependencies. Fig. 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, Fig.
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 Figs. 2(e) the angle of maximum intensity,
(f) the intensity at $\theta_{s}=45$°, and (g) and the angular width of the
emission pattern $\Delta\theta$, respectively, all as a function of the slab
thickness. In Fig. 2(e), the maximum emission angle exhibits a growing
oscillatory pattern with increasing the thickness. The oscillatory behavior
occurs around a mean angle of $\sim 45$°. At the detection angle of
$\theta_{s}=45$°, we plot in Fig. 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 $\Delta a\sim
0.29\lambda_{p}$. This period of oscillation is the same on all three plots of
Figs. 2(e-g). The angular uncertainty of emitted photons exhibits a similar
thickness-dependent variation. We define the angle uncertainty $\Delta\theta$
as the difference between the angles at 50$\%$ of the maximum intensity (full
width at half maximum). In Fig. 2(g), we again observe increasing oscillations
for thicknesses above $a=\lambda_{p}$, while it remains relatively constant
for subwavelength thicknesses, as also seen in Fig. 2(c).

It can be verified, that the peaks and valleys in Fig. 2(f), correspond to
constructive and destructive interference conditions, satisfying the relations
$\Delta=m\lambda_{s}$ and $\Delta=(m+\frac{1}{2})\lambda_{s}$, respectively,
where $m$ is a positive integer. $\Delta=2\,a\,n_{2s}\cos(\theta_{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 lipson .
$\theta_{2s}$ is the angle in medium 2 such that it yields $\theta_{1s}=45$°
in medium 1 after refraction. The period of the resonance peaks is found by
$\Delta a=\frac{\lambda_{s}}{2\,n_{2s}\cos(\theta_{2s})}$, which results in
$\Delta a\approx 0.29\lambda_{p}$, which matches the numerical results. These
points of constructive and destructive interferences for $\theta_{1s}=45$°
also coincide with peaks and valleys of $\Delta\theta$, see in Fig. 2(g). We
point out that the angles of constructive interference oscillate around the
$\theta_{1s}=45$°, as the thickness is increased, which results in an
oscillation of the angles of maximum intensity around $\theta_{1s}=45$° in
Fig. 2(e). This can also can be seen in Fig. 2(d), where at thicknesses of
$a=3.61\lambda_{p}$ or $a=3.72\lambda_{p}$, the directionality plots are
maximized at angles $\sim 58.5$° or $\sim 27.5$°, respectively.

A crucial observation is that due to the inherent losses in GaAs at
$\lambda_{p}=500$ nm, there are constrictions on longitudinal phase-matching
effect. Calculating the pump’s decay length inside the GaAs slab, denoted as
$\alpha_{p}=1/k^{\prime\prime}_{2p}$, where
$k_{2p}=k^{\prime}_{2p}+\mathrm{i}\,k^{\prime\prime}_{2p}$, yields
$\alpha_{p}=0.35\lambda_{p}$. When compared with the coherence length of our
system under normal incidence, $L_{c}=\pi/\Delta k_{\parallel}$, where $\Delta
k_{\parallel}=k^{\prime}_{2p}-k_{2s}\cos\theta_{2s}-k_{2i}\cos\theta_{2i}$, we
obtain $L_{c}=0.6\lambda_{p}$ for joint detection in medium 1 at
$\theta_{1s}=\theta_{1i}=45$°. Notably, the decay length is smaller than the
coherence length ($\alpha_{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 behaviour observed in Figs.
2(e-g), and not the longitudinal phase matching. Hence, the longitudinal phase
matching, only affects the height of the first few peaks in Fig. 2(f), while
subsequent peak intensities stabilize.

###  V.2 Non-degenerate photon-pair detection

Now we consider a frequency non-degenerate scenario for photon-pair detection.
Unlike the degenerate case, the $\varphi$-symmetric configuration is not
expected to dominate in the non-degenerate case. For a normally incident
plane-wave pump, with $\mathbf{q}_{p}=0$, the emitted pairs satisfy the
transversal phase-matching condition of $\mathbf{q}_{s}=-\mathbf{q}_{i}$. This
right away sets $\varphi_{i}=\varphi_{s}+\pi$ for signal and idler photons in
a pair. Then, taking $\theta$ as the emission angle with respect to the
forward $z$-direction we get $k_{x}^{2}+k_{y}^{2}=k^{2}\sin^{2}\theta$ for
both signal and idler. Using the degeneracy factor $r$ defined as
$r\equiv\lambda_{s}/\lambda_{i}$, we get
$|\sin\theta_{s}|=r\absolutevalue{n_{i}/n_{s}}|\sin\theta_{i}|$ 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 Fig. 3(a), we calculate the joint detection probability with a fixed signal
detector at $\theta_{s}=45$°, $\varphi_{s}=180$° for nondegeneracy factors
$r=0.8$, $1$, $1.2$, and $1.5$ in the ultra-thin case ($a=0.01\lambda_{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 $\theta_{i}(r=0.8)=62.11$°, $\theta_{i}(r=1)=45$°,
$\theta_{i}(r=1.2)=36.1$° and $\theta_{i}(r=1.5)=28.12$°, matching Fig.3(a).
For $r=1$ the maximum coincidence rate lie along $\theta_{s}=\theta_{i}$,
justifying the previously used $\varphi$-symmetric configuration.
Additionally, in Fig. 3(a), we observe that as the degeneracy factor $r$
decreases, the uncertainty in idler angles increases, meaning weaker
correlations in emitted photons.

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

Due to transverse phase matching, most pairs are generated under the condition
$\varphi_{s}=\varphi_{i}+180$°, which suggests visualizing the joint-detection
rate as a function of the angles $\\{\theta_{s},\theta_{i}\\}$, in a
particular plane defined by $\varphi$. In Fig. 3(b) we plot $\mathcal{R}$
along the $xz$-plane, where $\varphi_{i}=0$ and $\varphi_{s}=180$°, as a
function of $\theta_{s}$ and $\theta_{i}$ for $r=1.5$ and $a=0.01\lambda_{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 of the process. Additionally, emission at $\sim 45$°, dominant in
the degenerate case, is no longer prevalent in the non-degenerate scenario,
which are now shifted by the transverse phase-matching condition. In the
depicted scenario in Fig. 3(b), these angles are $\theta_{s}=55.35$° and
$\theta_{i}=33.3$°.

##  VI 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 polarization-independent 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 weissflog2024 . There, the density matrix $\hat{\rho}$,
expressed in the $\\{\ket{HH},\ket{HV},\ket{VH},\ket{VV}\\}$ basis, is
retrieved from projective measurements computed for a set of 16 tomographic
states $\\{\ket{\psi_{\nu}}\\}_{\nu\in[[1;16]}$ howell ; weissflog2024 . The
details of the tomographic approach can be found in the appendix (D).

Figure 4: (a) Real and imaginary parts of the elements of the density matrix,
$\hat{\rho}$, for the polarization state of a degenerate ($r=1$) photon pair,
generated from an ultra-thin GaAs slab of thickness $a=0.01\lambda_{p}$. The
detectors are in a $\varphi$-symmetric configuration with $\theta_{s}=45$° and
$\varphi_{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 $\varphi$-symmetric configuration, for two
thicknesses $a=0.01\lambda_{p}$ and $a=0.1\lambda_{p}$

Using this method, we examine the density matrix for several scenarios,
starting with an ultra-thin slab ($a=0.01\lambda_{p}$) in a
$\varphi$-symmetric configuration detection ($\theta_{s}=45$°,
$\varphi_{s}=0$) for $r=1$. The pump is the same as in Sec. (V). In Fig. 4(a),
we display the real and imaginary parts of the density matrix $\hat{\rho}$
representing the polarization state. Remarkably, under these conditions, the
density matrix corresponds to a maximally polarization-entangled state of the
form
$\ket{\psi}=\ket{H}_{s}\ket{V}_{i}-\ket{V}_{s}\ket{H}_{i}=\ket{HV}-\ket{VH}$.
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 unentangled 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\lambda_{p}$ and $a=0.1\lambda_{p}$), including non-
degeneracy, in Fig. 4(b) at the same $\varphi$-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\neq 1$ due to
moving away from the transversal phase-matching condition (as in Fig. 3). This
demonstrates GaAs’s ability to produce polarization-entangled Bell states of
the form $\ket{\psi}=\ket{HV}-\ket{VH}$ across a wide range of wavelengths,
maintaining near-perfect entanglement at $\varphi$-symmetric detection with
minor degradation at non-degenerate wavelengths.

It is interesting to explore the preservation of broad entanglement obtained
in the $\varphi$-symmetric setup across various emission angles, especially in
non-degenerate scenarios where the $\varphi$-symmetric configuration is not
dominant. Fig. 5 displays the Schmidt number ($K$) variation with
$\\{\theta_{i},\varphi_{i}\\}$, with fixed signal detection at
$\theta_{s}=45$° and $\varphi_{s}=0$, for both degenerate (Fig. 5(a, c)) and
non-degenerate cases with $r=1.5$ (Fig. 5(b, d)), with $a=0.01\lambda_{P}$ and
$a=1\lambda_{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 $\varphi$-symmetric angle. Notably,
although entanglement slightly degrades in the non-degenerate
$\varphi$-symmetric case (the values of $K\sim 1.98$ and $1.99$ in Figs. 5(b,
d)), there actually exists a close by angle at which the entanglement is
maximum again at $K\sim 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 non-degenerate 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. weissflog2024 ,
where similar responses were observed in Bell state generation in GaAs
nonlinear nanoresonators.

Figure 5: Schmidt number as a function of the idler angles
$\\{\theta_{i},\varphi_{i}\\}$ with fixed $\theta_{s}=45$° and $\varphi_{s}=0$
for (a) $a=0.01\lambda_{p}$ and $r=1$, (b) $a=1\lambda_{p}$ and $r=1.5$, (c)
$a=1\lambda_{p}$ and $r=1$, (c) $a=1\lambda_{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.

##  VII 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 down-conversion (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 far-
field properties of photon pairs detected within one of the media surrounding
the slab. Using a Zinc-Blende $\left\langle 100\right\rangle$ GaAs crystal as
a specific example, deposited on a SiO2 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
$\ket{\psi}=\ket{HV}-\ket{VH}$, within a $\varphi$-symmetric configuration by
analyzing its density matrix $\hat{\rho}$. 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 phase-matching 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.

##  VIII acknowledgements

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) and by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) – Project number
505897284.

##  Appendix A Far-Field Approximation of the Joint-Detection Probability

Let us consider a 2-dimensional Fourier expansion integral of the form

$\displaystyle I(x,y,z)$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{2}}\iint\differential k_{x}\,\differential
k_{y}\;U(k_{x},k_{y})\,\frac{e^{\mathrm{i}k_{z}z}}{k_{z}}\,e^{\mathrm{i}(k_{x}x+k_{y}y)}\,,$
(15)
$\displaystyle=\mathcal{F}^{-1}\left[U(k_{x},k_{y})\,\frac{e^{\mathrm{i}k_{z}z}}{k_{z}}\right]\,,$

that describes the amplitude propagation of a field with angular spectrum
$U(k_{x},k_{y})/k_{z}$ at the initial position
$\mathbf{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime}=0)$ to the final
position $\mathbf{r}=(x,y,z)$. With wave vector component in the $z$ direction
$k_{z}=\sqrt{k^{2}-k_{x}^{2}-k_{y}^{2}}$.

To do a far-field approximation of this propagation integral, we first use the
convolution theorem to write Eq. (A) as a convolution of the form

$\displaystyle I(x,y,z)$
$\displaystyle\propto\mathcal{F}^{-1}\left[U(k_{x},k_{y})\right]*\mathcal{F}^{-1}\left[\frac{e^{\mathrm{i}k_{z}z}}{k_{z}}\right]\,,$
(16) $\displaystyle\propto u(x,y)*\frac{e^{\mathrm{i}kr}}{r}\,,$ (17)
$\displaystyle\propto\iint\differential x^{\prime}\,\differential
y^{\prime}\,u(x^{\prime},y^{\prime})\,\frac{e^{\mathrm{i}kR}}{R}\,,$ (18)

where $*$ denotes the convolution operation,
$U(k_{x},k_{y})=\mathcal{F}[u(x,y)]$ and we have used the Weyl expansion of an
outgoing spherical wave as a superposition of plane waves,
$\frac{e^{\mathrm{i}kr}}{r}=\frac{\mathrm{i}}{2\pi}\iint\differential
k_{x}\,\differential
k_{y}\;\frac{e^{\mathrm{i}k_{z}z}}{k_{z}}\,e^{\mathrm{i}(k_{x}x+k_{y}y)}$.
Where $r=\sqrt{x^{2}+y^{2}+z^{2}}$ is the magnitude of the position vector
$\mathbf{r}$ and $R=\sqrt{(x-x^{\prime})^{2}+(y-y^{\prime})^{2}+z^{2}}$ is the
magnitude of the vector $\mathbf{R}=\mathbf{r}-\mathbf{r}^{\prime}$.

Eq. (18) is interpreted as the beam propagation of an initial field
$u(x^{\prime},y^{\prime})$ to a distance $z$ with an outgoing spherical wave
as the response function $h={e^{\mathrm{i}kr}}/{r}$. In the spherical
coordinates representation with $x=r\sin\theta\cos\varphi$,
$y=r\sin\theta\sin\varphi$ and $z=r\cos\theta$, $R$ takes the form

$\displaystyle R=$ $\displaystyle
r\bigg{(}1-2\sin\theta\cos\varphi\frac{x^{\prime}}{r}-2\sin\theta\sin\varphi\frac{y^{\prime}}{r}+$
(19) $\displaystyle\frac{x^{\prime 2}}{r^{2}}+\frac{y^{\prime
2}}{r^{2}}\bigg{)}^{1/2}\,.$

By taking the condition that $\mathbf{r}\gg\mathbf{r^{\prime}}$ the last to
terms are approximated to be zero, and using a first order binomial expansion
of $R$, we arrive at

$\displaystyle R\approx
r\left(1-\sin\theta\cos\varphi\frac{x^{\prime}}{r}-\sin\theta\sin\varphi\frac{y^{\prime}}{r}\right)\,.$
(20)

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 is a
good enough approximation for the denominator. Thus, we get the far-field-
approximated expression for beam propagation,

$\displaystyle I(r,\theta,\varphi)$
$\displaystyle\propto\frac{e^{\mathrm{i}kr}}{r}\iint\differential
x^{\prime}\,\differential y^{\prime}\,u(x^{\prime},y^{\prime})\,$ (21)
$\displaystyle\times
e^{-\mathrm{i}(k\sin\theta\cos\varphi\,x^{\prime}+k\sin\theta\sin\varphi\,y^{\prime})}\,,$
$\displaystyle\propto\frac{e^{\mathrm{i}kr}}{r}\,U(k\sin\theta\cos\varphi,k\sin\theta\sin\varphi)\,,$
(22)

which tells that, under the far-field approximation, the amplitude propagation
to the far-field of an initial beam field at $\mathbf{r}^{\prime}$ with
angular spectrum $U(k_{x},k_{y})/k_{z}$ is proportional to $U(k_{x},k_{y})$
with $k_{x}=k\sin\theta\cos\varphi$ and $k_{y}=k\sin\theta\sin\varphi$.

This approximation, which is similar to the Fraunhofer approximation for the
diffraction of classical beams born , 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. (II), and notice that the integral
inside can be written as

$\displaystyle I_{\mathcal{R}}=$
$\displaystyle\int_{-\infty}^{\infty}\differential\mathbf{q}_{s}\;e^{\mathrm{i}k_{z,1s}z_{s}}$
$\displaystyle\times\int_{-\infty}^{\infty}\,\differential\mathbf{q}_{i}\;\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},\omega_{s},\omega_{p}-\omega_{s})e^{\mathrm{i}k_{z,1i}z_{i}}\;$
$\displaystyle\times
e^{\mathrm{i}(\mathbf{q}_{s}\cdot\mathbf{r}_{s\perp}+\mathbf{q}_{i}\cdot\mathbf{r}_{i\perp})}\,,$
(23)

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^{\mathrm{i}k_{z,1s}z_{s}}$
and $e^{\mathrm{i}k_{z,1i}z_{i}}$ come from the definition of the GF in Eq.
(III), which allows us to write

$\displaystyle\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},$
$\displaystyle z_{s},z_{i},\omega_{s},\omega_{p}-\omega_{s})=$
$\displaystyle\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},\omega_{s},\omega_{p}-\omega_{s})e^{\mathrm{i}k_{z,1s}z_{s}}e^{\mathrm{i}k_{z,1i}z_{i}}\,.$
(24)

The integral in Eq. (A) is a two-photon analogue of the propagation integral
in Eq. (A) and can be interpreted as the propagation of the probability
amplitude of the generated signal and idler photon pair with joint angular
probability (JAP)
$\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},\omega_{s},\omega_{p}-\omega_{s})$
at $\mathbf{r}^{\prime}$ to the detection position at $\mathbf{r}$.

Therefore, we can perform a similar two-photon-like far-field approximation
where we identify $\tilde{\mathcal{R}}$ as $U/k_{z}$ in Eq. (A), and by using
the identity in (22) twice (for signal and idler propagation), we obtain

$I_{\mathcal{R}}\propto\frac{e^{\mathrm{i}k_{s}r_{s}}}{r_{s}}\frac{e^{\mathrm{i}k_{i}r_{i}}}{r_{i}}\tilde{\mathcal{R}}(\mathbf{q}_{s},\mathbf{q}_{i},\omega_{s},\omega_{p}-\omega_{s})\,k_{z,1s}k_{z,1i}\,,$
(25)

with $k_{x,s}=k_{1s}\sin\theta_{s}\cos\varphi_{s}$ and
$k_{y,s}=k_{1s}\sin\theta_{s}\sin\varphi_{s}$ and similarly for the idler
photon. After inserting (25) in (II) we arrive at the expression in (II) which
describes a far-field-approximated joint-probability rate of detecting
coincidences in the far-field as a function of the detection angles
$\\{\theta_{s},\theta_{i},\varphi_{s},\varphi_{i}\\}$.

##  Appendix B Green’s Function of a Dielectric Slab

The GFs in Fourier domain for the multilayered system depicted in Fig. (1),
where a field is generated in medium 2 and propagates outside the slab can be
constructed as

$\displaystyle\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{g}}_{21}(\mathbf{q},z,z^{\prime},\omega)=\frac{-\mathrm{i}}{2k_{z2}}\bigg{(}T_{21}^{(s)}(q,z^{\prime},\omega)\;\hat{s}\hat{s}+$
$\displaystyle
T_{21}^{(p+)}(q,z^{\prime},\omega)\;\hat{p}_{1+}\hat{p}_{2+}+T_{21}^{(p-)}(q,z^{\prime},\omega)\;\hat{p}_{1+}\hat{p}_{2-}\bigg{)}e^{\mathrm{i}k_{z1}z}\,,$
(26)

and

$\displaystyle\overset{\text{\tiny$\leftrightarrow$}}{\mathbf{g}}_{23}(\mathbf{q},z,z^{\prime},\omega)=\frac{-\mathrm{i}}{2k_{z2}}\big{[}T_{23}^{(s)}(q,z^{\prime},\omega)\;\hat{s}\hat{s}+$
$\displaystyle
T_{23}^{(p-)}(q,z^{\prime},\omega)\;\hat{p}_{3-}\hat{p}_{2-}+T_{23}^{(p+)}(q,z^{\prime},\omega)\;\hat{p}_{3-}\hat{p}_{2+}\big{]}e^{-\mathrm{i}k_{z3}z}\,,$
(27)

where we have also added the GF connecting medium 2 and 3 for completeness.

Using the definitions for the $s$\- and $p$-polarized fields in Eqs. (8) and
(9) we can compute the dyadics $\hat{s}\hat{s}$, $\hat{p}_{1+}\hat{p}_{2\pm}$
and $\hat{p}_{3-}\hat{p}_{2\pm}$ as

$\displaystyle\hat{s}\hat{s}$
$\displaystyle=\frac{1}{q^{2}}\left(\begin{matrix}k_{y}^{2}&-k_{x}k_{y}&0\\\
-k_{x}k_{y}&k_{x}^{2}&0\\\ 0&0&0\end{matrix}\right)\;,$ (28)
$\displaystyle\hat{p}_{1+}\hat{p}_{2\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{k_{1}k_{2}}\left(\begin{matrix}\pm\frac{k_{x}^{2}k_{z1}k_{z2}}{q^{2}}&\pm\frac{k_{x}k_{y}k_{z1}k_{z2}}{q^{2}}&-k_{x}k_{z1}\\\
\pm\frac{k_{x}k_{y}k_{z1}k_{z2}}{q^{2}}&\pm\frac{k_{y}^{2}k_{z1}k_{z2}}{q^{2}}&-k_{y}k_{z1}\\\
\mp k_{x}k_{z2}&\mp k_{y}k_{z2}&q^{2}\end{matrix}\right)\,,$ (29)
$\displaystyle\hat{p}_{3-}\hat{p}_{2\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{k_{2}k_{3}}\left(\begin{matrix}\mp\frac{k_{x}^{2}k_{z2}k_{z3}}{q^{2}}&\mp\frac{k_{x}k_{y}k_{z2}k_{z3}}{q^{2}}&k_{x}k_{z3}\\\
\mp\frac{k_{x}k_{y}k_{z2}k_{z3}}{q^{2}}&\mp\frac{k_{y}^{2}k_{z2}k_{z3}}{q^{2}}&k_{y}k_{z3}\\\
\mp k_{x}k_{z2}&\mp k_{y}k_{z2}&q^{2}\end{matrix}\right)\,.$ (30)

The quantities $T^{(s)}$ and $T^{(p)}$ correspond to the generalized
transmission coefficients for the $s$\- and $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$,

$\displaystyle
r_{ij}^{(p)}=\frac{k_{zi}\epsilon_{j}-k_{zj}\epsilon_{i}}{k_{zi}\epsilon_{j}+k_{zj}\epsilon_{i}}\;,\quad
r_{ij}^{(s)}=\frac{k_{zi}-k_{zj}}{k_{zi}+k_{zj}}\;,$ (31) $\displaystyle
t^{(p)}_{ij}=\frac{2n_{i}n_{j}k_{zi}}{k_{zi}\epsilon_{j}+k_{zj}\epsilon_{i}}\;,\qquad
t^{(s)}_{ij}=\frac{2k_{zi}}{k_{zi}+k_{zj}}\,.$ (32)

Thus, $T^{(s,p)}$ can be calculated as Sipe ; lau2007

$T_{21}(q,z^{\prime})=\frac{\left[e^{-\mathrm{i}k_{z2}z^{\prime}}+r_{23}e^{\mathrm{i}k_{z2}(2a+z^{\prime})}\right]}{1-r_{23}r_{21}e^{\mathrm{i}k_{z2}(2a)}}\times
t_{21}\,,$ (33)

and

$T_{23}(q,z^{\prime})=\frac{\left[e^{\mathrm{i}k_{z2}(a+z^{\prime})}+r_{21}e^{\mathrm{i}k_{z2}(a-z^{\prime})}\right]}{1-r_{23}r_{21}e^{\mathrm{i}k_{z2}(2a)}}\times
t_{23}\,,$ (34)

where the corresponding Fresnel coefficient for the $s$\- and $p$-polarized
vectors must be used.

As explained in Sec. (III), 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
$p$-polarized field one can write,

$T^{(p)}_{21}(q,z^{\prime})=T^{(p+)}_{21}(q,z^{\prime})+T^{(p-)}_{21}(q,z^{\prime})\,,$
(35)

which splits the transmission into an upward field generated in medium 2 that
transmits to medium 1, namely $T^{(p+)}_{21}$, and a downward wave generated
in medium 2 that also transmits to medium 1, namely $T^{(p-)}_{21}$. From Eq.
(33) one can distinguish,

$\displaystyle
T^{(p+)}_{21}(q,z^{\prime})=\frac{e^{-\mathrm{i}k_{z2}z^{\prime}}}{1-r^{(p)}_{23}r^{(p)}_{21}e^{\mathrm{i}k_{z2}(2a)}}\times
t^{(p)}_{21}\,,$ (36) $\displaystyle
T^{(p-)}_{21}(q,z^{\prime})=\frac{r^{(p)}_{23}e^{\mathrm{i}k_{z2}(2a+z^{\prime})}}{1-r^{(p)}_{23}r^{(p)}_{21}e^{\mathrm{i}k_{z2}(2a)}}\times
t^{(p)}_{21}\,.$ (37)

This distinction becomes important when we want to construct the GF of our
system and it was not considered in Ref. lau2007 .

For completeness,the same analysis is applied to the GF in Eq. (B) for waves
generated in medium 2 and transmitting to medium 3. Here,

$\displaystyle
T_{23}^{(p+)}(q,z^{\prime})=\frac{r_{21}^{(p)}e^{\mathrm{i}k_{z2}(a-z^{\prime})}}{1-r^{(p)}_{23}r^{(p)}_{21}e^{\mathrm{i}k_{z2}(2a)}}\times
t^{(p)}_{23}\,,$ (38) $\displaystyle
T_{23}^{(p-)}(q,z^{\prime})=\frac{e^{\mathrm{i}k_{z2}(a+z^{\prime})}}{1-r^{(p)}_{23}r^{(p)}_{21}e^{\mathrm{i}k_{z2}(2a)}}\times
t^{(p)}_{23}\,,$ (39)

where similarly, as depicted in Fig. 1(c), $T^{(p+)}_{23}$ describes an upward
$p$-wave generated in medium 2 that transmits to medium 3 after multiple
reflections and $T^{(p-)}_{23}$ 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 (IV).
Using $\lambda_{p}=500$ nm as the pump wavelength, we can observe in Fig. 6
the intensity of the spectrum, $x$\- and $z$-component of the field in Fourier
domain for two spectral pump widths $w$. For a width $w=6.6\times 10^{5}$ m-1
Fig. 6(a) shows that the spectrum of spatial frequencies are much smaller than
the wave-vector ($q<<k_{p}=2\pi/\lambda_{p}$), this is the so-called
nonparaxial regime (weakly focused beam in position space). Here,
$U_{p}(\mathbf{q})$ corresponds to the angular spectrum of a Gaussian beam of
width $W=3\mu$m defined as $u_{p}(x,y)\propto 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
Fig.6(a) where the maximum Max$[E^{2}_{P,z}(\mathbf{q})]$ is about $\sim
10^{4}$ weaker than Max$[E^{2}_{P,x}(\mathbf{q})]$. 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 Fig. 6(b) we show the same
quantities for a width $w=2\times 10^{7}$ m-1. This is a nonparaxial scenario
where the angular spectrum saturates to values close to $\sim 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.

Figure 6: Intensity of the spectrum, $x$\- and $z$-polarized components of a
pump field of wavelength $\lambda_{p}=500$ nm in (a) for $w=6.6\times 10^{5}$
m-1 (paraxial regime) and (b) for a width $w=2\times 10^{7}$ m-1 (nonparaxial
regime).

###  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,
$\mathbf{E}_{p}(\mathbf{q},z=-a)=(E_{p,x},0,E_{p,z})$, where its longitudinal
component $E_{p,z}$ is calculated from from $E_{p,x}$ using the transverse
condition $\mathbf{k}_{p}\cdot\mathbf{E}_{p}=0$. This results in

$E_{p,z}(\mathbf{q},z=-a)=-\frac{k_{x}}{k_{z,3p}}E_{p,x}(\mathbf{q},z=-a)\,.$
(40)

We express the pump field at the interface in terms of the $s$\- and
$p$-polarized fields as

$\displaystyle\mathbf{E}_{p}$ $\displaystyle(\mathbf{q},z=-a)$
$\displaystyle={E}_{p,x}(\mathbf{q},z=-a)\hat{x}+{E}_{p,z}(\mathbf{q},z=-a)\hat{z}$
$\displaystyle={E}^{(s)}_{p}(\mathbf{q},z=-a)\hat{s}+{E}^{(p)}_{p}(\mathbf{q},z=-a)\hat{p}_{3+}\,,$
(41)

where we used the definitions of $\hat{s}$ and $\hat{p}$ in (8) and (9) for
medium 3. After solving we find that

$\displaystyle
E_{p}^{(s)}(\mathbf{q},z=-a)=\frac{k_{z,3p}k_{y}}{qk_{x}}{E}_{p,z}(\mathbf{q},z=-a)\,,$
(42) $\displaystyle
E_{p}^{(p)}(\mathbf{q},z=-a)=\frac{k_{3p}}{q}{E}_{p,z}(\mathbf{q},z=-a)\,,$
(43)

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 $x$-polarized, that is
$\mathbf{E}_{P}(\mathbf{q}=0,z=-a)={E}_{P,x}(\mathbf{q}=0,z=-a)\hat{x}$.

To describe the transmission of the field to the slab, one can calculate a
corresponding coefficient associated to the $s$\- and $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
take the form

$Q_{32}(q,z)=\frac{e^{\mathrm{i}k_{z,2p}(a+z)}(1+r_{21}e^{-\mathrm{i}k_{z,2p}(2z)})}{1-r_{23}r_{21}e^{\mathrm{i}k_{z,2p}(2a)}}\times
t_{32}\,,$ (44)

where the first term

$Q^{(p+)}_{32}(q,z)=\frac{e^{\mathrm{i}k_{z,2p}(a+z)}}{1-r_{23}r_{21}e^{\mathrm{i}k_{z,2p}(2a)}}\times
t_{32}\,,$ (45)

describes upward waves and the second term

$Q^{(p-)}_{32}(q,z)=\frac{r_{21}e^{\mathrm{i}k_{z,2p}(a-z)}}{1-r_{23}r_{21}e^{\mathrm{i}k_{z,2p}(2a)}}\times
t_{32}\,,$ (46)

downward waves inside the slab. This distinction is again only important for
the $p$-polarized fields, since a $p$-field changes direction when the field
propagates upward and downward (see Eq. (9)), while the $s$-field remains
unchanged. Notice that the factors are $z$-dependent and consider the multiple
reflections of the pump field in the slab.

Thus, we arrive at Eq. (IV), where the pump field inside the slab is
$\mathbf{E}_{p}(x,y,z)=\frac{1}{(2\pi)^{2}}\int\differential\mathbf{q}\;\mathbf{E}_{p}(\mathbf{q},z)\,e^{\mathrm{i}\mathbf{q\cdot
r}_{\perp}}\,,$ with

$\displaystyle\mathbf{E}_{p}(\mathbf{q},z)=Q^{(s)}(q,z)E^{(s)}_{p}(\mathbf{q},z=-a)\hat{s}$
$\displaystyle+E^{(p)}_{p}(\mathbf{q},z=-a)\left(Q^{(p+)}(q,z)\hat{p}_{2+}+Q^{(p-)}(q,z)\hat{p}_{2-}\right)\,,$
(47)

where again the $s$\- and $p$\- directions for each plane wave in the
expansion are defined as in Eqs. (8,9) for medium 2. Note that, for a plane
inside the slab located at $z\in\\{-a,0\\}$, the term
$Q^{(s)}(q,z)E^{(s)}_{p}(\mathbf{q},z=-a)\hat{s}$ in equation (IV) contains
the sum of all the direction-independent $s$-polarized fields, while the
direction-dependent term
$E^{(p)}_{p}(\mathbf{q},z=-a)\left(Q^{(p+)}(q,z)\hat{p}_{2+}+Q^{(p-)}(q,z)\hat{p}_{2-}\right)$
takes all $p_{+}$ (upward) and $p_{-}$ (downward) fields that add up all after
transmission to the slab and multiple reflections.

Finally, using the definitions in Eqs. (8,9), we can insert the computed Eqs.
(42) and (43) into Eq. (IV) to express the pump field inside the slab in
cartesian coordinates.

###  C.3 Generalized transmission coefficient for the pump field inside the
slab

A field coming from region 3 transmits into region 2 at the interface $z=-a$,
then it propagates to $z$ gaining a phase factor of
$\alpha=e^{ik_{z}(a-\absolutevalue{z})}$. The field then propagates to $z=0$
and reflects back to $z$ gaining an extra phase of
$\beta=e^{ik_{z}\absolutevalue{z}}r_{21}e^{ik_{z}\absolutevalue{z}}$. The
filed now propagates to $z=-a$ and reflects again to $z$ with an extra phase
of
$\gamma=e^{ik_{z}(a-\absolutevalue{z})}r_{23}e^{ik_{z}(a-\absolutevalue{z})}$.
Finally this process is repeated indefinitely due to multiple reflections. The
field inside the slab at $z$ can then be describe through the recurrence

$\displaystyle
Q^{(s,p)}(q,z)=(\alpha+\beta\alpha+\gamma\beta\alpha+\gamma\beta^{2}\alpha+\gamma^{2}\beta^{2}\alpha$
$\displaystyle+\gamma^{2}\beta^{3}\alpha+\gamma^{3}\beta^{3}\alpha+...)t_{32}\,,$
(48)

which converges to

$Q^{(s,p)}(q,z)=\frac{\alpha(1+\beta)}{1-\beta\gamma}t_{32}\,,$ (49)

leading to Eq. (44). If we split the recurrence into the upward and downward
waves at $z$ we find $Q^{(p+)}(q,z)=\frac{\alpha}{1-\beta\gamma}t_{32}\,,$ and
$Q^{(p-)}(q,z)=\frac{\alpha\beta}{1-\beta\gamma}t_{32}\,,$

##  Appendix D Quantum Polarization Tomography

Following Ref. James , the classical Stokes parameters can be generalized to
measure the state of multiple photon beams. We use the
$\\{\ket{HH},\ket{HV},\ket{VH},\ket{VV}\\}$ basis to express the polarization
density matrix $\hat{\rho}$ which is reconstructed by 16 projections into a
set of states $\\{\ket{\psi_{\nu}}\\}_{\nu\in[\\![1;16]\\!]}$. We employ the
same basis states $\ket{\psi_{\nu}}$ as in Ref. James . The orientation of the
orthogonal basis vectors $\ket{H}$ and $\ket{V}$ can in principle be chosen
freely, however will influence the shape of the density matrix (but
importantly not the degree of entanglement, which is independent of the
basis). We choose here a convention for $\ket{H}$ and $\ket{V}$ that relates
to an experimental scenario, where photon pairs are first collimated and then
projected into different measurement bases by polarizers fixed in the
laboratory coordinate frame $\hat{x}=\ket{H}$, $\hat{y}=\ket{V}$. Following
this, we define a new set of vectors $\\{\hat{x}^{\prime},\hat{y}^{\prime}\\}$
in the plane $\\{\hat{\theta},\hat{\varphi}\\}$, which for a propagation
direction $\varphi$ follow from

$\begin{cases}\hat{x}^{\prime}&=\cos{\varphi}\,\hat{\theta}-\sin{\varphi}\,\hat{\varphi}\\\
\hat{y}^{\prime}&=\sin{\varphi}\,\hat{\theta}+\cos{\varphi}\,\hat{\varphi}\end{cases}.$
(50)

When collected and collimated as in an experimental scheme, these vectors
coincide with the lab coordinate frame and therefore
$\hat{x}^{\prime}=\ket{H}$, $\hat{y}^{\prime}=\ket{V}$.

The Schmidt entanglement parameter K (Schmidt number) quantifying the degree
of entanglement of the pair is given by grobe ; eberly :

$K=\frac{1}{Tr\left(\hat{\rho}_{i}\right)^{2}}=\frac{1}{Tr\left(\hat{\rho}_{s}\right)^{2}}\,,$
(51)

where $\hat{\rho}_{s}$ and $\hat{\rho}_{i}$ are the reduced density matrices
corresponding to the single particle wavefunctions which are obtained from the
partial trace of the biphoton density matrix:
$\hat{\rho}_{s,i}=Tr_{s,i}(\hat{\rho})$.

## References

  * (1) C. K. Hong and L. Mandel, “Theory of parametric frequency down conversion of light,” Phys. Rev. A, vol. 31, pp. 2409–2418, Apr 1985.
  * (2) F. Flamini, N. Spagnolo, and F. Sciarrino, “Photonic quantum information processing: a review,” Reports on Progress in Physics, vol. 82, p. 016001, nov 2018.
  * (3) C. H. Bennett and G. Brassard, “Quantum cryptography: Public key distribution and coin tossing,” Theoretical Computer Science, vol. 560, pp. 7–11, 2014\.  Theoretical Aspects of Quantum Cryptography – celebrating 30 years of BB84.
  * (4) S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, D. Bunandar, R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani, J. L. Pereira, M. Razavi, J. S. Shaari, M. Tomamichel, V. C. Usenko, G. Vallone, P. Villoresi, and P. Wallden, “Advances in quantum cryptography,” Adv. Opt. Photon., vol. 12, pp. 1012–1236, Dec 2020.
  * (5) M. Gilaberte Basset, F. Setzpfandt, F. Steinlechner, E. Beckert, T. Pertsch, and M. Gräfe, “Perspectives for applications of quantum imaging,” Laser & Photonics Reviews, vol. 13, no. 10, p. 1900097, 2019.
  * (6) P.-A. Moreau, E. Toninelli, T. Gregory, and M. J. Padgett, “Imaging with quantum states of light,” Nature Reviews Physics, vol. 1, no. 6, pp. 367–380, 2019.
  * (7) A. S. Clark, M. Chekhova, J. C. F. Matthews, J. G. Rarity, and R. F. Oulton, “Special Topic: Quantum sensing with correlated light sources,” Applied Physics Letters, vol. 118, p. 060401, 02 2021.
  * (8) R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, “Quantum entanglement,” Rev. Mod. Phys., vol. 81, pp. 865–942, Jun 2009.
  * (9) Y. Wang, K. D. Jöns, and Z. Sun, “Integrated photon-pair sources with nonlinear optics,” Applied Physics Reviews, vol. 8, p. 011314, 03 2021\.
  * (10) P. G. Kwiat, K. Mattle, H. Weinfurter, A. Zeilinger, A. V. Sergienko, and Y. Shih, “New high-intensity source of polarization-entangled photon pairs,” Phys. Rev. Lett., vol. 75, pp. 4337–4341, Dec 1995.
  * (11) P. G. Kwiat, E. Waks, A. G. White, I. Appelbaum, and P. H. Eberhard, “Ultrabright source of polarization-entangled photons,” Phys. Rev. A, vol. 60, pp. R773–R776, Aug 1999.
  * (12) M. Hentschel, H. Hübel, A. Poppe, and A. Zeilinger, “Three-color sagnac source of polarization-entangled photon pairs,” Opt. Express, vol. 17, pp. 23153–23159, Dec 2009.
  * (13) G. Saerens, T. Dursap, I. Hesner, N. M. H. Duong, A. S. Solntsev, A. Morandi, A. Maeder, A. Karvounis, P. Regreny, R. J. Chapman, A. Danescu, N. Chauvin, J. Penuelas, and R. Grange, “Background-free near-infrared biphoton emission from single gaas nanowires,” Nano Letters, vol. 23, no. 8, pp. 3245–3250, 2023.  PMID: 37057961.
  * (14) G. Marino, A. S. Solntsev, L. Xu, V. F. Gili, L. Carletti, A. N. Poddubny, M. Rahmani, D. A. Smirnova, H. Chen, A. Lemaître, G. Zhang, A. V. Zayats, C. D. Angelis, G. Leo, A. A. Sukhorukov, and D. N. Neshev, “Spontaneous photon-pair generation from a dielectric nanoantenna,” Optica, vol. 6, pp. 1416–1422, Nov 2019.
  * (15) T. Santiago-Cruz, A. Fedotova, V. Sultanov, M. A. Weissflog, D. Arslan, M. Younesi, T. Pertsch, I. Staude, F. Setzpfandt, and M. Chekhova, “Photon pairs from resonant metasurfaces,” Nano Letters, vol. 21, no. 10, pp. 4423–4429, 2021.  PMID: 33971095.
  * (16) T. Santiago-Cruz, S. D. Gennaro, O. Mitrofanov, S. Addamane, J. Reno, I. Brener, and M. V. Chekhova, “Resonant metasurfaces for generating complex quantum states,” Science, vol. 377, no. 6609, pp. 991–995, 2022.
  * (17) J. Zhang, J. Ma, M. Parry, M. Cai, R. Camacho-Morales, L. Xu, D. N. Neshev, and A. A. Sukhorukov, “Spatially entangled photon pairs from lithium niobate nonlocal metasurfaces,” Science Advances, vol. 8, no. 30, p. eabq4240, 2022\.
  * (18) C. Son, V. Sultanov, T. Santiago-Cruz, A. P. Anthur, H. Zhang, R. Paniagua-Dominguez, L. Krivitsky, A. I. Kuznetsov, and M. V. Chekhova, “Photon pairs bi-directionally emitted from a resonant metasurface,” Nanoscale, vol. 15, pp. 2567–2572, 2023.
  * (19) C. Okoth, A. Cavanna, T. Santiago-Cruz, and M. V. Chekhova, “Microscale generation of entangled photons without momentum conservation,” Phys. Rev. Lett., vol. 123, p. 263602, Dec 2019.
  * (20) T. Santiago-Cruz, V. Sultanov, H. Zhang, L. A. Krivitsky, and M. V. Chekhova, “Entangled photons from subwavelength nonlinear films,” Opt. Lett., vol. 46, pp. 653–656, Feb 2021.
  * (21) V. Sultanov, T. Santiago-Cruz, and M. V. Chekhova, “Flat-optics generation of broadband photon pairs with tunable polarization entanglement,” Opt. Lett., vol. 47, pp. 3872–3875, Aug 2022.
  * (22) Q. Guo, X.-Z. Qi, L. Zhang, M. Gao, S. Hu, W. Zhou, W. Zang, X. Zhao, J. Wang, B. Yan, et al., “Ultrathin quantum light source with van der waals nbocl2 crystal,” Nature, vol. 613, no. 7942, pp. 53–59, 2023.
  * (23) M. A. Weissflog, A. Fedotova, Y. Tang, E. A. Santos, B. Laudert, S. Shinde, F. Abtahi, M. Afsharnia, I. P. Pérez, S. Ritter, et al., “A tunable transition metal dichalcogenide entangled photon-pair source,” arXiv preprint arXiv:2311.16036, 2023.
  * (24) C. Trovatello, C. Ferrante, B. Yang, J. Bajo, B. Braun, X. Xu, Z. H. Peng, P. K. Jenke, A. Ye, M. Delor, D. N. Basov, J. Park, P. Walther, L. A. Rozema, C. Dean, A. Marini, G. Cerullo, and P. J. Schuck, “Quasi-phase-matched up- and down-conversion in periodically poled layered semiconductors,” arXiv preprint arXiv:2312.05444, 2023.
  * (25) A. Vega, E. A. Santos, J. Fuenzalida, M. Gilaberte Basset, T. Pertsch, M. Gräfe, S. Saravi, and F. Setzpfandt, “Fundamental resolution limit of quantum imaging with undetected photons,” Phys. Rev. Res., vol. 4, p. 033252, Sep 2022.
  * (26) E. A. Santos, T. Pertsch, F. Setzpfandt, and S. Saravi, “Subdiffraction quantum imaging with undetected photons,” Phys. Rev. Lett., vol. 128, p. 173601, Apr 2022.
  * (27) A. Krasnok, S. Lepeshov, and A. Alú, “Nanophotonics with 2d transition metal dichalcogenides,” Opt. Express, vol. 26, pp. 15972–15994, Jun 2018\.
  * (28) A. N. Poddubny, I. V. Iorsh, and A. A. Sukhorukov, “Generation of photon-plasmon quantum states in nonlinear hyperbolic metamaterials,” Phys. Rev. Lett., vol. 117, p. 123901, Sep 2016.
  * (29) T. Gruner and D.-G. Welsch, “Green-function approach to the radiation-field quantization for homogeneous and inhomogeneous kramers-kronig dielectrics,” Phys. Rev. A, vol. 53, pp. 1818–1829, Mar 1996.
  * (30) H. T. Dung, L. Knöll, and D.-G. Welsch, “Three-dimensional quantization of the electromagnetic field in dispersive and absorbing inhomogeneous dielectrics,” Phys. Rev. A, vol. 57, pp. 3931–3942, May 1998.
  * (31) J. E. Sipe, “New green-function formalism for surface optics,” J. Opt. Soc. Am. B, vol. 4, pp. 481–489, Apr 1987.
  * (32) J. Peřina, “Chapter 2 - spontaneous parametric down-conversion in nonlinear layered structures,” in Progress in Optics (E. Wolf, ed.), vol. 59 of Progress in Optics, pp. 89–158, Elsevier, 2014.
  * (33) J. Peřina, “Spatial properties of entangled photon pairs generated in nonlinear layered structures,” Phys. Rev. A, vol. 84, p. 053840, Nov 2011.
  * (34) C. K. Hong and L. Mandel, “Theory of parametric frequency down conversion of light,” Phys. Rev. A, vol. 31, pp. 2409–2418, Apr 1985.
  * (35) J. Schneeloch and J. C. Howell, “Introduction to the transverse spatial correlations in spontaneous parametric down-conversion through the biphoton birth zone,” Journal of Optics, vol. 18, no. 5, p. 053501, 2016.
  * (36) M. A. Weissflog, R. Dezert, V. Vinel, C. Gigli, G. Leo, T. Pertsch, F. Setzpfandt, A. Borne, and S. Saravi, “Nonlinear nanoresonators for Bell state generation,” Applied Physics Reviews, vol. 11, p. 011403, 01 2024\.
  * (37) W. T. Lau, J.-T. Shen, G. Veronis, and S. Fan, “Spatial coherence of the thermal electromagnetic field in the vicinity of a dielectric slab,” Phys. Rev. E, vol. 76, p. 016601, Jul 2007.
  * (38) M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light.  Elsevier Science, 2013.
  * (39) K. Papatryfonos, T. Angelova, A. Brimont, B. Reid, S. Guldin, P. R. Smith, M. Tang, K. Li, A. J. Seeds, H. Liu, and D. R. Selviah, “Refractive indices of mbe-grown alxga(1$-$x)as ternary alloys in the transparent wavelength region,” AIP Advances, vol. 11, p. 025327, 02 2021.
  * (40) I. H. Malitson, “Interspecimen comparison of the refractive index of fused silica$\ast$,†,” J. Opt. Soc. Am., vol. 55, pp. 1205–1209, Oct 1965\.
  * (41) S. Lipson, H. Lipson, and D. Tannhauser, Optical Physics.  Optical Physics, Cambridge University Press, 1995.
  * (42) D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White, “Measurement of qubits,” Phys. Rev. A, vol. 64, p. 052312, Oct 2001.
  * (43) R. Grobe, K. Rzazewski, and J. H. Eberly, “Measure of electron-electron correlation in atomic physics,” Journal of Physics B: Atomic, Molecular and Optical Physics, vol. 27, p. L503, aug 1994.
  * (44) J. Eberly, “Schmidt analysis of pure-state entanglement,” Laser physics, vol. 16, pp. 921–926, 2006.