# General simulation method for spontaneous parametric down- and parametric
up-conversion experiments

Felix Riexinger<EMAIL_ADDRESS>Fraunhofer Institute for
Industrial Mathematics ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany
Department of Physics and Research Center OPTIMAS, Technische Universität
Kaiserslautern (TUK), 67663 Kaiserslautern, Germany  Mirco Kutas  Fraunhofer
Institute for Industrial Mathematics ITWM, Fraunhofer-Platz 1, 67663
Kaiserslautern, Germany  Department of Physics and Research Center OPTIMAS,
Technische Universität Kaiserslautern (TUK), 67663 Kaiserslautern, Germany
Björn Haase  Fraunhofer Institute for Industrial Mathematics ITWM, Fraunhofer-
Platz 1, 67663 Kaiserslautern, Germany  Department of Physics and Research
Center OPTIMAS, Technische Universität Kaiserslautern (TUK), 67663
Kaiserslautern, Germany  Patricia Bickert  Fraunhofer Institute for Industrial
Mathematics ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany  Daniel
Molter  Fraunhofer Institute for Industrial Mathematics ITWM, Fraunhofer-Platz
1, 67663 Kaiserslautern, Germany  Michael Bortz  Fraunhofer Institute for
Industrial Mathematics ITWM, Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany
Georg von Freymann  Fraunhofer Institute for Industrial Mathematics ITWM,
Fraunhofer-Platz 1, 67663 Kaiserslautern, Germany  Department of Physics and
Research Center OPTIMAS, Technische Universität Kaiserslautern (TUK), 67663
Kaiserslautern, Germany

###### Abstract

Spontaneous parametric down-conversion (SPDC) sources are an important
technology for quantum sensing and imaging. We demonstrate a general
simulation method, based on modeling from first principles, reproducing the
spectrally and spatially resolved absolute counts of a SPDC experiment. By
simulating additional processes and effects we accomplish good agreement with
the experimental results. This method is broadly applicable and allows for the
separation of contributing processes, virtual characterization of SPDC
sources, and enables the simulation of many quantum based applications.

††preprint: APS/123-QED

Entangled and correlated photon pairs have become the basis for many
applications in quantum optics. They are used in various quantum based schemes
such as ghost imaging [1, 2], optical coherence tomography [3], spectroscopy
[4, 5], quantum sensing [6], and imaging with undetected photons [7, 8, 9].
One of the most prominent methods for generating entangled photon pairs is
spontaneous parametric down-conversion (SPDC), where a pump laser photon
decays into two lower-energy photons in a medium with a second-order
nonlinearity. The entanglement can exist in many of the photon properties,
such as polarization, wavelength, or momentum, making SPDC a versatile source
for entangled photons. Additionally, the process is easy to implement and well
understood experimentally. However, an accurate simulation method for SPDC
sources is required as the foundation for simulations of many applications
with entangled photons. The availability of detailed simulations can help
resolve problems such as finding limitations to the resolution or visibility
of quantum imaging [8, 9]. This issue is becoming more relevant as
applications for SPDC sources are on the verge of a breakthrough, but their
limiting factors need to be understood better to exploit the full potential of
such quantum based applications.

The theory of SPDC is well developed and multiple approaches to simulating the
properties of the created photons are available [10, 11]. However, these SPDC
simulations are often limited to the paraxial regime [11], or do not predict
absolute photon conversion rates [12]. In this letter we propose and
demonstrate a novel simulation method for SPDC sources and the subsequent
measurement setup. The sparse use of approximations makes the underlying model
applicable to a wide range of SPDC sources from the ultraviolet to the
terahertz regime.

Our method reproduces the spectrally and spatially resolved absolute photon
count rates. This is demonstrated on an experiment with idler photons in the
terahertz range and signal photons in the visible range. The extreme
wavelength spread between signal and idler leads to a setup that covers a
large range in frequency and emission directions and further has multiple
quasi-phasematching (QPM) orders overlapping in the same wavelength range. In
the terahertz range additional processes such as parametric up- and
(nonspontaneous) down-conversion occur parallel to SPDC. In order to
adequately match the experimental data we include these processes in our
simulation.

The high qualitative and quantitative agreement with experimental results
demonstrates the capabilities of our simulation method even for complex SPDC
sources.

Our model is based on the second-order nonlinear interaction of
electromagnetic fields together with a first-order perturbation theory
approximation. We start with the Hamiltonian [13]

$H_{\mathrm{NL}}(t)=\frac{1}{3}\varepsilon_{0}\int\mathrm{d}\mathbf{r}\zeta_{jkl}^{(2)}(\mathbf{r})\mathbf{D}_{j}(\mathbf{r},t)\mathbf{D}_{k}(\mathbf{r},t)\mathbf{D}_{l}(\mathbf{r},t),$
(1)

where $\zeta_{jkl}^{(2)}$ is the second-order inverse susceptibility tensor
and $\mathbf{D}$ are the displacement fields. The formulation with the
displacement fields instead of the electric fields is necessary to ensure
consistency after quantization [14, 15].

We describe the pump beam as a classical monochromatic Gaussian beam with
linear polarization. We assume that the pump is undepleted, meaning its power
is large enough such that it is not significantly reduced due to nonlinear
conversion events. In addition, we use the approximation of a collimated beam
such that the curvature and the Gouy phase can be neglected. The pump
propagates along the z-axis, which we define parallel to the optical axis of
the system. The signal and idler fields are described using a plane-wave
decomposition separated into positive and negative frequency components
$\hat{\mathbf{D}}^{+}$ and $\hat{\mathbf{D}}^{-}$ with

$\hat{\mathbf{D}}^{+}(\mathbf{r},t)=\sum_{\mathbf{k},\sigma}i\sqrt{\frac{\varepsilon_{0}n_{\mathbf{k}}^{2}\hbar\omega_{\mathbf{k}}}{2V}}\hat{a}_{\mathbf{k},\sigma}\mathbf{\epsilon}_{\mathbf{k},\sigma}e^{i(\mathbf{k}\cdot\mathbf{r}+\omega_{\mathbf{k}}t)},$
(2)

and $\hat{\mathbf{D}}^{-}$ being the hermitian conjugate of
$\hat{\mathbf{D}}^{+}$. Here, $V$ is the quantization volume,
$\hat{a}_{\mathbf{k},\sigma}$ is the annihilation operator for a photon with
momentum $\mathbf{k}$, and $\mathbf{\epsilon}_{\mathbf{k},\sigma}$ is the
direction of the displacement field vector indexed with the polarization
$\sigma$.

We then calculate the two-photon state $\ket{\psi(t)}$ using first-order
perturbation theory:

$\displaystyle\ket{\psi(t)}\approx\mathopen{}\mathclose{{}\left(1-\frac{i}{\hbar}\int_{0}^{t}\mathrm{d}\tau
H_{\mathrm{NL}}(\tau)}\right)\ket{\psi(0)},$ (3)

where the initial state $\ket{\psi(0)}$ is a vacuum state. Under the
assumption that the quantization volume is large, we can make a transition
from sums to integrals in Eq. (2). From this we obtain the signal count rate
density:

$\displaystyle\Gamma_{\mathrm{d}}(\mathbf{k}_{\mathrm{s}})$
$\displaystyle=\frac{1}{T_{\mathrm{I}}}\bra{\psi(T_{\mathrm{I}})}\hat{a}^{\dagger}(\mathbf{k}_{\mathrm{s}})\hat{a}(\mathbf{k}_{\mathrm{s}})\ket{\psi(T_{\mathrm{I}})}$
(4) $\displaystyle=Z\sum_{\sigma_{\mathrm{s}},\sigma_{\mathrm{i}}}\sum_{m\
\mathrm{odd}}\int\mathrm{d}k_{\mathrm{i}}^{3}\mathopen{}\mathclose{{}\left\|A(\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})}\right\|^{2},$
(5)

with

$\displaystyle A($
$\displaystyle\mathbf{k}_{\mathrm{s}},\mathbf{k}_{\mathrm{i}})=\frac{\chi^{(2)}_{\mathrm{eff}}}{m}\sqrt{\frac{\omega_{\mathrm{s}}\omega_{\mathrm{i}}}{n_{\mathrm{s}}^{2}n_{\mathrm{i}}^{2}}}\,\operatorname{sinc}\mathopen{}\mathclose{{}\left[\frac{1}{2}\Delta
k_{z}L}\right]$
$\displaystyle\times\exp\mathopen{}\mathclose{{}\left[-\frac{1}{4}(\Delta
k_{x}^{2}+\Delta
k_{y}^{2})w_{\mathrm{p}}^{2}}\right]\,\operatorname{sinc}\mathopen{}\mathclose{{}\left[\frac{1}{2}\Delta\omega
T_{\mathrm{I}}}\right]$ (6)

and

$Z=\frac{16Pw_{\mathrm{p}}^{2}L^{2}T_{I}}{(2\pi)^{7}\varepsilon_{0}n_{\mathrm{p}}c}.$
(7)

Here $T_{\mathrm{I}}$ is the interaction time for a pump photon with the
nonlinear medium, $m$ denotes the odd QPM orders, and we sum over permutations
of the indices $j$, $k$, $l$ and substitute new indices to denote pump (p),
signal (s) and idler (i). $\Delta
k_{z}=\mathbf{k}_{\mathrm{p}z}-\mathbf{k}_{\mathrm{s}z}-\mathbf{k}_{\mathrm{i}z}+k_{\Lambda}$
and $\Delta k_{j}=\mathbf{k}_{\mathrm{s}j}+\mathbf{k}_{\mathrm{i}j}$ with
$j=x,y$ are the longitudinal and transversal phase mismatches, and
$\Delta\omega=\omega_{\mathrm{p}}-\omega_{\mathrm{s}}-\omega_{\mathrm{i}}$
corresponds to the energy mismatch. The width of the transverse part is
determined by the waist radius $w_{\mathrm{p}}$ of the pump beam. The periodic
poling offset $k_{\mathrm{\Lambda}}=2m\pi/\Lambda$ depends on the poling
period of the crystal $\Lambda$ and the QPM order $m$. $P$ denotes the power
of the pump beam and $L$ the length of the crystal. The refractive indices $n$
and $\chi^{(2)}_{\mathrm{eff}}$ are functions of the wave vectors
$\mathbf{k}_{\mathrm{p}},\,\mathbf{k}_{\mathrm{s}},\,\mathbf{k}_{\mathrm{i}}$
and their corresponding polarizations. The spectral dependence of the
nonlinear coefficient is modeled with Miller’s rule [16]. The spatial
variation of $\chi^{(2)}_{\mathrm{eff}}$ is considered by calculating the
effective value from the tensor components and displacement field directions
[17]. The spectral dependence of the refractive indices is modeled with
Sellmeier equations in the visible range and polynomial fits to experimental
data from literature in the terahertz range. Spatial variation in the
refractive indices is considered following [18]. The spectral dependencies are
especially relevant since we simulate a large spectral range in the terahertz
regime, from $\sim 0.1\ \mathrm{THz}$ to $3.6\ \mathrm{THz}$. The spatial
variation cannot be neglected here, because the transverse momentum
conservation dictates a large emission angle range for the idler photons.

With the idler in the terahertz range, at room temperature, thermal photons at
the idler wavelength have to be taken into account [19].

These thermal photons interact with the pump laser as well. Through parametric
down-conversion, additional photons at the signal wavelength are created. We
derive this process analogously to SPDC, where instead of an initial vacuum
state we have a thermal state for the idler. This leads to the same
phasematching properties, but instead of a $"1"$ from the vacuum fluctuations
we obtain the thermal fluctuations

$N_{\mathrm{th}}=\frac{1}{\text{exp}(\hbar\omega_{\mathrm{i}}/k_{\mathrm{B}}T_{\mathrm{c}})-1},$
(8)

where $T_{\mathrm{c}}$ is the temperature of the crystal. The influence of
thermal fluctuations for the signal can be neglected.

The emitted spectrum is propagated from the crystal onto a detector through an
optical system. The number of counts for a single detector pixel is obtained
by integrating over all signal wave vectors that are propagated to this pixel.
However, not all photons created in the crystal arrive at the detector. We
consider the losses from internal reflection in the crystal, during
propagation through the optical setup and the efficiency of the detector
summarized into a single factor $\eta$. With this we obtain the photon counts
for the detector pixel with indices $(i,j)$

$R_{\mathrm{d}}^{(i,j)}=\eta
T(1+N_{\mathrm{th}})\int_{\Omega(i,j)}\mathrm{d}\mathbf{k}_{\mathrm{s}}\Gamma_{\mathrm{d}}(\mathbf{k}_{\mathrm{s}}),$
(9)

where $T$ is the illumination time of the detector and $\Omega(i,j)$ is
defined such that every ray with $\mathbf{k}_{\Omega}\in\Omega(i,j)$ at the
crystal exit is propagated through the optical system onto the detector pixel
$(i,j)$. Along the same lines we obtain the rate for the parametric up-
conversion process, where a thermal photon and a laser photon merge into a
signal photon. The only changes are the signs of the terahertz frequency
$\omega_{\mathrm{i}}$ in the terms $\Delta\omega$ and the wave vector
$\mathbf{k}_{\mathrm{i}j}$ in the respective $\Delta k_{j}$ terms with
$j=x,y,z$. Up-conversion is only caused by thermal fluctuations, such that the
up-conversion rate is given by:

$R_{\mathrm{u}}^{(i,j)}=\eta
TN_{\mathrm{th}}\int_{\Omega(i,j)}\mathrm{d}\mathbf{k}_{\mathrm{s}}\Gamma_{\mathrm{u}}(\mathbf{k}_{\mathrm{s}}).$
(10)

Figure 1: Layout of the experimental setup. Only the components used in the
simulations are included. f: lens. PPLN: periodically poled MgO-doped
LiNbO3-crystal. M: mirror. Length of optical paths are not to scale. (a)
represents a cross section of the spectrum after the slit. (b) sketches the
separation of spectral components by the grating. The ratio
$d_{\mathrm{a}}/d_{\mathrm{b}}$ quantifies the asymmetry introduced in this
step. Figure 2: Experimental (a) and simulated (b) frequency-angular spectrum.
The average dark count rate is subtracted from the experimental spectrum. The
dashed lines in (a) denote the location of the cuts presented in Fig. 3 and 4.
The numbers $m$ in (b) indicate the contributions of the different QPM orders.
The tails on the left correspond to up-conversion, the ones on the right to
down-conversion.

A sketch of the model for the experimental setup is shown in Fig. 1, where
only the simulated components are shown. The setup consists of a narrow
bandwidth laser, a nonlinear crystal and an imaging system. The spatial and
spectral components of the created signal photons are separated by the imaging
system and imaged onto the detector. The imaging system contains a slit to
limit the transmission of rays with large $\mathbf{k}_{x}$, shown in Fig. 1
(a), resulting in a sharper image on the detector. This effect is shown in
Fig. 1 (b), where the remaining parts of the ellipses are separated along the
$\theta_{\mathrm{s}}$ axis. A transmission grating separates the spectral
components of the photons. Photons with a fixed wavelength are emitted in a
cone shape, shown in Fig. 1 (a). This leads to an imperfect separation of
spatial and spectral components, which is shown by the full and partial
ellipses in Fig. 1 (b). The ratio $d_{\mathrm{a}}/d_{\mathrm{b}}$ is a measure
of how much an ellipse is squeezed. As it approaches zero, the ellipses are
imaged as vertical lines, making the spread in $\theta_{\mathrm{s}}$ large and
the one in $\lambda_{\mathrm{s}}$ small. The system further contains several
filters and Bragg gratings to suppress the pump radiation that are not
considered in this model. A detailed description of the experimental
realization is given by Haase et al. [20]. The laser is modeled with a
wavelength of $\lambda_{\mathrm{p}}=659.58\ \mathrm{nm}$ and a beamwaist
$w_{\mathrm{p}}=43\ \mathrm{\mu m}$. The nonlinear crystal is a periodically
poled MgO-doped lithium niobate crystal with dimensions $5\times 1\times 10\
\mathrm{mm}^{3}$ $(H\times W\times L)$ and a poling period of $170\
\mathrm{\mu m}$. Due to symmetries in the nonlinear susceptibility only the
values of the $\chi^{(2)}_{333}$ and $\chi^{(2)}_{311}$ components need to be
considered. The $\chi^{(2)}_{222}$ component is not relevant to any processes
in our setup. We use a value of $\chi^{(2)}_{333}=327\ \mathrm{pm}/\mathrm{V}$
at $661\ \mathrm{nm}$ and $0.75\ \mathrm{THz}$ which we obtained from a fit to
the experimental results. The fitted $\chi^{(2)}_{333}$ value agrees with a
scaled value from the infrared range [21] but is $\sim 1.5$ times larger than
other scaled values from the terahertz range [22, 23]. The difference can be
explained by different factors for the Hamiltonian ($\frac{1}{2}$ and
$\frac{1}{3}$) used in the classical and quantum approaches. Compared to our
previous estimation [20] this value is two times larger. There are two reasons
for this: First, we did not use Miller’s rule in the previous method, but
assumed a constant value. Second, we added the full optical model to the
simulation, which changes the shape of the spectrum. The value of
$\chi^{(2)}_{311}=49\ \mathrm{pm}/\mathrm{V}$ at $661\ \mathrm{nm}$ and $0.75\
\mathrm{THz}$ is a scaled value from [21]. The contributions of this parameter
are small such that a reliable fit is not possible. The refractive indices,
for $\sim 5\ \mathrm{mol.}\%$ MgO-doped LiNbO3, are taken from [24, 25].

The propagation through the optical system is modeled by paraxial ray optics
assuming ideal optical components. The paraxial approximation is reasonable as
only the signal photons are detected. The idler photons which are emitted at
larger angles are not considered. Since the longitudinal positions of the
optical components are not exactly known, they are estimated from the imaging
properties of the setup. This is done using an optimization algorithm to
minimize the squared difference of three experimentally measured values: the
transformation of wavelength into a position on the x-axis, the relation
between emission angle and position on the y-axis, and the squeezing of a
monochromatic circular beam, defined by the ratio
$d_{\mathrm{a}}/d_{\mathrm{b}}$. We also penalize deviations from the measured
positions. This procedure of estimating the parameters of the setup mimics the
alignment process performed in the experiment, where the components are moved
around their nominal positions to obtain a sharper image. The calculated
values deviate from the nominal distances of the setup, but reproduce the
imaging properties of the actual experiment. Further sources of errors in the
optical system such as misalignment transverse to the optical axis or
derivations from nominal values are not included in our model.

We employ a Monte-Carlo integration scheme to evaluate the integrals in Eq.
(4), approximating each squared sinc functions with the sum of three scaled
Gaussian functions. This allows for efficient sampling while maintaining a
small error in the approximated function. Other options such as a single
Gaussian or a Cosine Gauss approximation [26] are either less accurate or less
efficient. The sum in Eq. (4) is evaluated up to the seventh QPM order. Since
positive and negative orders are possible, a total of eight summands are
evaluated. Higher orders contribute less than $70$ counts per pixel to the
spectrum in the observed range and are therefore neglected. We consider both
ordinary and extraordinary polarization for signal and idler. Type 0
parametric conversion contributes most of the spectrum in the investigated
range, type I contributes less than $165$ counts per pixel in the region above
$663.7\ \mathrm{nm}$ and less than $2$ counts per pixel elsewhere. Type II
contributes less than $1$ count per pixel in the whole range.

The experimental and numerical results for the full spectrum are shown in Fig.
2. The experimental spectrum is a single image as recorded by the sCMOS camera
with the average dark count rate subtracted. Simulated and experimental
spectrum show four distinct tails for up- as well as down-conversion. The
contributing QPM orders for each tail are given in panel (b) of the figure.
The tails corresponding to the fifth QPM order are barely visible, while the
contributions of the seventh order cannot be distinguished at all due to the
lower conversion efficiency of higher orders. The experiment is limited to a
scattering angle of around $\pm 2.3^{\circ}$ by the apertures of the optical
components. The simulated spectrum extends beyond this angle as the limiting
apertures are not included in the model.

The simulation shows more counts than the experiment at wavelengths lower than
$655\ \mathrm{nm}$ and higher than $663.5\ \mathrm{nm}$ for up- and down-
conversion, respectively. Potential reasons are limiting apertures in the
experiment or an incorrect model for the terahertz refractive indices. The
idler photons in this range have a frequency of over $2.5\ \mathrm{THz}$ which
is beyond the measured range of our reference, covering frequencies from $0.3\
\mathrm{THz}$ to $1.9\ \mathrm{THz}$. The value of $\chi^{(2)}_{\mathrm{eff}}$
depends on the refractive indices in our model such that an overestimation of
$n$ increases the counts.

Figure 3: Horizontal cut through the experimental and simulated spectra at
$\theta_{\mathrm{s}}=0.2^{\circ}$. The positive and negative sign for
$\theta_{\mathrm{s}}$ correspond to cuts through the upper and lower half of
Fig. 2 respectively. The simulated spectrum is symmetric in
$\theta_{\mathrm{s}}$, thus only one line is shown.  Figure 4: Vertical cuts
through the experimental and simulated spectra. Panel (a) shows a cut in the
up-conversion regime, panel (b) in the down-conversion regime. The color
coding is the same as in Fig. 3.

Figures 3 and 4 show horizontal and vertical cuts of the spectrum. The
positive and negative values for $\theta_{\mathrm{s}}$ correspond to the sign
given in Fig. 2. As the simulated spectrum is symmetric in
$\theta_{\mathrm{s}}$ only one line is shown. Both cuts show good agreement in
position, width and height of the peaks. The experimental spectrum shows some
remaining pump light around the laser wavelength at the center of Fig. 3. This
is more pronounced for the $+\theta_{\mathrm{s}}$ case. In the model we assume
the pump to be blocked completely, therefore the simulated spectrum does not
show the pump rest.

The tails of the experimental spectrum are slightly narrower and more
pronounced, which is due to the experiment being adjusted for a maximally
sharp image, while the simulated setup is not. It is optimized to reproduce
three imaging characteristics of the experiment. Observed differences between
simulation and experiment are of the same order of magnitude as the
differences between the two experimental results. Peaks found in the cuts from
larger angles or larger wavelengths match in shape and position. The peaks of
the simulated spectrum are significantly higher for this region. Note that the
simulated spectrum exhibits some minor modulation along the
$\theta_{\mathrm{s}}$-axis which are numerical artifacts caused by the
employed sampling method.

The demonstrated simulation of parametric conversion spectra shows good
agreement with the experiment in spectral and angular distributions as well as
absolute photon counts. Qualitative and quantitative features of the
experimentally obtained spectrum can be reproduced. Various effects such as
SPDC, the influence of thermal photons and parametric up-conversion were
simulated. This shows the potential of applying our model to predict accurate
characteristics of photon sources. And the possibility of identifying the
contributions of different processes. The model can be simplified for faster
computation times. An important step to improve the simulation results is to
use better estimates for the crystal characteristics in the investigated
frequency range. The method provides significant benefits over traditional
methods such as the calculation of phasematching curves as no information
about the width or intensity of the spectrum is provided there. Our numerical
method also allows for reconstruction of the spectrum at the crystal face
which allows to investigate spatial, spectral and correlation properties
without the need of building optical setups for measuring them. Compared to
previous results [20], the propagation of the spectrum through the measurement
setup improves the simulation results. Further research is needed to separate
influences of model errors in the up- and down-conversion model from those in
the optical measurement setup. Nonetheless, the presented model provides the
necessary basis for the simulation of many quantum optic applications such as
quantum imaging.

## References

  * Pittman  _et al._ [1995] T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, Optical imaging by means of two-photon quantum entanglement, Phys. Rev. A 52, R3429 (1995).
  * Padgett and Boyd [2017] M. J. Padgett and R. W. Boyd, An introduction to ghost imaging: quantum and classical, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20160233 (2017).
  * Vanselow  _et al._ [2020] A. Vanselow, P. Kaufmann, I. Zorin, B. Heise, H. M. Chrzanowski, and S. Ramelow, Frequency-domain optical coherence tomography with undetected mid-infrared photons, Optica 7, 1729 (2020).
  * Lindner  _et al._ [2021] C. Lindner, J. Kunz, S. J. Herr, S. Wolf, J. Kießling, and F. Kühnemann, Nonlinear interferometer for fourier-transform mid-infrared gas spectroscopy using near-infrared detection, Opt. Express 29, 4035 (2021).
  * Kutas  _et al._ [2021] M. Kutas, B. Haase, J. Klier, D. Molter, and G. von Freymann, Quantum-inspired terahertz spectroscopy with visible photons, Optica 8, 438 (2021).
  * Kutas  _et al._ [2020] M. Kutas, B. Haase, P. Bickert, F. Riexinger, D. Molter, and G. von Freymann, Terahertz quantum sensing, Science Advances 6, eaaz8065 (2020).
  * Lemos  _et al._ [2014] G. B. Lemos, V. Borish, G. D. Cole, S. Ramelow, R. Lapkiewicz, and A. Zeilinger, Quantum imaging with undetected photons, Nature 512, 409 (2014).
  * Fuenzalida  _et al._ [2020] J. Fuenzalida, A. Hochrainer, G. B. Lemos, E. Ortega, R. Lapkiewicz, M. Lahiri, and A. Zeilinger, Resolution of quantum imaging with undetected photons, arXiv preprint arXiv:2010.07712 (2020).
  * Gilaberte Basset  _et al._ [2021] M. Gilaberte Basset, A. Hochrainer, S. Töpfer, F. Riexinger, P. Bickert, J. R. León-Torres, F. Steinlechner, and M. Gräfe, Video-rate imaging with undetected photons, Laser & Photonics Rev. 15, 2000327 (2021).
  * Schneeloch and Howell [2016] J. Schneeloch and J. C. Howell, Introduction to the transverse spatial correlations in spontaneous parametric down-conversion through the biphoton birth zone, J. Opt. 18, 053501 (2016).
  * Trajtenberg-Mills  _et al._ [2020] S. Trajtenberg-Mills, A. Karnieli, N. Voloch-Bloch, E. Megidish, H. S. Eisenberg, and A. Arie, Simulating correlations of structured spontaneously down-converted photon pairs, Laser & Photonics Rev. 14, 1900321 (2020).
  * Ramírez-Alarcón  _et al._ [2013] R. Ramírez-Alarcón, H. Cruz-Ramírez, and A. B. U’Ren, Effects of crystal length on the angular spectrum of spontaneous parametric downconversion photon pairs, Laser Phys. 23, 055204 (2013).
  * Schneeloch  _et al._ [2019] J. Schneeloch, S. H. Knarr, D. F. Bogorin, M. L. Levangie, C. C. Tison, R. Frank, G. A. Howland, M. L. Fanto, and P. M. Alsing, Introduction to the absolute brightness and number statistics in spontaneous parametric down-conversion, J. Opt. 21, 043501 (2019).
  * Hillery and Mlodinow [1984] M. Hillery and L. D. Mlodinow, Quantization of electrodynamics in nonlinear dielectric media, Phys. Rev. A 30, 1860 (1984).
  * Quesada and Sipe [2017] N. Quesada and J. E. Sipe, Why you should not use the electric field to quantize in nonlinear optics, Opt. Lett. 42, 3443 (2017).
  * Miller [1964] R. C. Miller, Optical second harmonic generation in piezoelectric crystals, Appl. Phys. Lett. 5, 17 (1964).
  * Boyd [2008] R. W. Boyd,  _Nonlinear Optics_ , 3rd ed. (Academic Press, Inc., USA, 2008) pp. 40–46.
  * Born and Wolf [2019] M. Born and E. Wolf,  _Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light_, 60th ed. (Cambridge University Press, 2019) pp. 806–807.
  * Kitaeva  _et al._ [2011] G. K. Kitaeva, S. Kovalev, A. Penin, A. Tuchak, and P. Yakunin, A method of calibration of terahertz wave brightness under nonlinear-optical detection, J. Infrared Millim. Terahertz Waves 32, 1144 (2011).
  * Haase  _et al._ [2019] B. Haase, M. Kutas, F. Riexinger, P. Bickert, A. Keil, D. Molter, M. Bortz, and G. von Freymann, Spontaneous parametric down-conversion of photons at 660 nm to the terahertz and sub-terahertz frequency range, Opt. Express 27, 7458 (2019).
  * Nikogosyan [2006] D. N. Nikogosyan,  _Nonlinear optical crystals: a complete survey_ (Springer Science & Business Media, 2006) pp. 41–44.
  * Hebling  _et al._ [2004] J. Hebling, A. G. Stepanov, G. Almási, B. Bartal, and J. Kuhl, Tunable THz pulse generation by optical rectification of ultrashort laser pulses with tilted pulse fronts, Appl. Phys. B 78, 593 (2004).
  * Vodopyanov [2008] K. Vodopyanov, Optical THz-wave generation with periodically-inverted GaAs, Laser & Photonics Rev. 2, 11 (2008).
  * Wu  _et al._ [2015] X. Wu, C. Zhou, W. R. Huang, F. Ahr, and F. X. Kärtner, Temperature dependent refractive index and absorption coefficient of congruent lithium niobate crystals in the terahertz range, Opt. Express 23, 29729 (2015).
  * Gayer  _et al._ [2008] O. Gayer, Z. Sacks, E. Galun, and A. Arie, Temperature and wavelength dependent refractive index equations for MgO-doped congruent and stoichiometric LiNbO3, Appl. Phys. B 91, 343 (2008).
  * Kolenderski  _et al._ [2009] P. Kolenderski, W. Wasilewski, and K. Banaszek, Modeling and optimization of photon pair sources based on spontaneous parametric down-conversion, Phys. Rev. A 80, 013811 (2009).