# Spectro-temporally tailored Non-Gaussian Quantum Operations in Thin-Film Waveguides

Peter Namdar, Patrick Folge, Carlos Lopetegui, Silia Babel, Benjamin Brecht, Christine Silberhorn, Valentina Parigi

► To cite this version:

Peter Namdar, Patrick Folge, Carlos Lopetegui, Silia Babel, Benjamin Brecht, et al. Spectro-temporally tailored Non-Gaussian Quantum Operations in Thin-Film Waveguides. 2025. hal-05362338

HAL Id: hal-05362338

https://hal.science/hal-05362338v1

Preprint submitted on 13 Nov 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Spectro-temporally tailored Non-Gaussian Quantum Operations in
Thin-Film Waveguides

Peter Namdar¹, Patrick Folge², Carlos E. Lopetegui¹, Silia Babel²,
Benjamin Brecht², Christine Silberhorn², Valentina Parigi¹

¹Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, 4 place Jussieu, F-75252 Paris, France

²Integrated Quantum Optics, Institute for Photonic Quantum Systems (PhoQS), Paderborn University, Warburgerstr. 100, Paderborn 33098, Germany

*Corresponding author: peter.namdar@lkb.upmc.fr

Advancements in photonic platforms have enabled the precise control of light’s spectral and temporal degrees of freedom, a capability crucial for the development of scalable quantum information systems. In this work, we address the challenge of implementing spectro-temporal mode-selective non-Gaussian quantum operations, specifically single-photon subtraction (SPS) and addition (SPA), in the telecom wavelength regime. Building on prior experimental demonstrations of mode-selective near-infrared SPS, we present the first design framework for achieving mode-selective SPA and SPS using thin-film lithium niobate nonlinear waveguide platforms. We introduce an inverse-design optimization scheme by modeling the quantum-optical response via the Joint Spectral Amplitude (JSA) and Transfer Function (TF), in order to identify optimal waveguide and pump parameters that maximize mode selectivity and state purity. This approach is first tested on a metallic waveguide design. We then exploit the dispersion engineering capabilities of thin-film waveguides, which offer enhanced nonlinear interactions through tighter light confinement. Our findings demonstrate that tailored nonlinear processes, particularly parametric down-conversion and frequency up-conversion, can support high-fidelity non-Gaussian operations essential for next-generation quantum photonic networks.

# 1. Introduction

Photonics platforms provide powerful and versatile tools for encoding and transmitting information in both classical and quantum protocols. To fully harness the potential of light multiplexing—encoding across the complete set of optical degrees of freedom—comprehensive spectro-temporal control is essential for advancing quantum applications [1, 2]. In particular, multimode spectro-temporal encoding for quantum information processing necessitates precise manipulation of quantum operations in both the spectral and temporal domains. While significant progress has been made in generating spectrally and temporally multiplexed quantum states with Gaussian field statistics [3–6] operations exploiting non-Gaussian quantum field statistics remain costly and challenging to control [7].

Achieving spectro-temporal mode-selective operations at the single-photon level requires highly tailored nonlinear optical processes. Thin-film waveguides now serve as an effective platform for these processes [8–10]. However, optimizing fabrication parameters to fine-tune the nonlinear response introduces a complex inverse-design problem. In this work, we identify a set of optimal parameters for designing spectro-temporal mode-selective single-photon subtraction (SPS) and addition (SPA) operations at telecom wavelengths via an optimization scheme. Previous experiments have successfully demonstrated SPA and SPS using nonlinear media [11–14], typically performed directly on the whole input fields to analyze their resulting (non-)classicality. While mode-selective SPS was experimentally realized in the near-infrared by Ra et al. in 2017 [15], our work extends this concept into the telecom regime and additionally proposes the first implementation of mode-selective SPA in general. While waveguides with micrometer-scale cross-sections have been implemented in various quantum optical experiments [16, 17], recent technological advancements have introduced an additional trend toward so-called thin-film waveguides. These structures confine light to significantly smaller cross-sections, often in the sub-micron scale, thereby greatly enhancing nonlinear interactions [8, 18]. Recent studies have shown that by carefully tailoring the waveguide geometry, the spectral properties of the interacting fields can be precisely engineered to achieve optimal phase-matching conditions, maximize mode
overlap, and enhance conversion efficiency [19–21]. Advances in manufacturing and calibration tech-
niques enable correction of fabrication imperfections, thereby making it feasible to realize theoretical
configurations, such as those proposed in this study [22]. Additionally, significant progress has been
made in developing low-loss thin-film waveguides in the telecom-wavelength regime, further enhancing
the potential for high-performance integrated photonic systems [23–27].

Starting from a fully analytical model that has been shown to be a sufficient guideline, we could later
apply our schematics to thin-film lithium niobate as described in Sec. 2.4.b. For a quantum-optical
description of the proposed experimental scenarios, we express the so-called Joint Spectral Amplitude
(JSA) or Transfer Function (TF) of the nonlinear scheme, which is determined by the waveguide pa-
rameters and the incoming classical light field (see Sec. 2).

To optimize these waveguide and pump-defining parameters, we introduce an algorithmic optimiza-
tion approach designed to efficiently identify optimal configurations within the parameter space (see
Methods). In this context, we establish the necessary conditions for valid solutions and analyze the
resulting purity of the non-Gaussian states, providing both a description and evaluation. We focus
entirely on the well-studied regime of second-order nonlinear processes ($\chi^{(2)}$-processes), specifically
parametric down-conversion (PDC) and frequency up-conversion. To model the mode-selective sce-
narios, we adopt established approaches, evolving the incoming states either through a beam splitter
or via a nonlinear PDC Hamiltonian [28–30].

Given the established infrastructure of classical communication networks and the recent demonstration
of multimode cluster states in the telecom regime, we focus our approach on applying non-Gaussian pro-
cesses to light at a central wavelength of $\lambda_{center} = 1560$ nm [31]. However, this scheme is wavelength-
independent and can, in principle, be applied at any wavelength.

For clarification, in this work we always sketch a non-collinear configuration involving three interact-
ing beams. In reality, these beams are assumed to be well-confined within a single spatial mode of a
waveguide, where they can be distinguished by either wavelength or polarization.

## 2. Non-Gaussian processes on multimode quantum states

The two processes of SPA and SPS rely on nonlinear interactions in a $\chi^{(2)}$-medium, where an
input quantum state interacts with a strong classical beam, facilitating controlled photon addition or
subtraction in specific spectral-temporal modes. In the following, we provide a detailed description of
each process, highlighting their theoretical foundations and experimental feasibility.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Schematic representations of the experimental setups discussed in this work. For clarity, the optical paths are illustrated in a non-collinear configuration, although in well-confined waveguides, they are collinear and can be separated by wavelength and/or polarization. The two processes are heralded by photon detection on one of the paths.</figcaption></figure>
## 2.1. Mode-Selective Single-Photon Addition

Fig. 1a illustrates a setup for performing single-photon addition on an incoming quantum state $\hat{\rho}$. In this configuration, the two incident fields — the signal $\hat{\rho}$ and a strong coherent pump — interact within a $\chi^{(2)}$-medium, generating a third beam via parametric down-conversion. The signal state is associated with quantum mode $|a\rangle$, with the creation and annihilation operators $\{\hat{a}^\dagger, \hat{a}\}$, while the up-converted beam is generated in the idler mode $|b\rangle$, associated with $\{\hat{b}^\dagger, \hat{b}\}$, which is initially prepared in the vacuum state $|0\rangle$. The system's nonlinear optical dynamics are governed by the Hamiltonian:

$$ \hat{H} = \iint d\omega_s d\omega_i J(\omega_s, \omega_i) \hat{a}^\dagger(\omega_s) \hat{b}^\dagger(\omega_i) + \text{h.c.} $$

where $J(\omega_s, \omega_i)$ is the Joint Spectral Amplitude (JSA), describing the spectral correlations between the two emerging quantum fields. $J(\omega_s, \omega_i)$ can be expressed as the product of the pump beam's spectral shape $\alpha(\omega_s + \omega_i)$ and the crystal's phase-matching function $\phi(\omega_s, \omega_i)$, where $\omega_s$ and $\omega_i$ are the frequencies of the signal and idler photons [31]. Since the JSA encodes the spectral correlations between these fields, analyzing its structure is crucial for understanding the mode composition of the generated quantum state and further the mode-selectivity of the process. A natural approach to achieve this is through Schmidt decomposition, which is mathematically equivalent to a singular-value decomposition (SVD) of the JSA. This decomposition provides eigenmodes and corresponding eigenvalues $\lambda_l$, quantifying the contributions of specific mode pairs $\{h_l(\omega_s), g_l(\omega_i)\}$ to the overall state. To diagonalize the JSA to obtain the eigenstates and their contractub, in 2000 C.K. Law extended the Schmidt decomposition framework to the spectro-temporal domain [32].

$$ J(\omega_s, \omega_i) = \alpha(\omega_s + \omega_i)\phi(\omega_s, \omega_i) = \sum_l \sqrt{\lambda_l} h_l(\omega_s) g_l(\omega_i) \quad (1) $$

The role of the phase-matching function $\phi(\omega_s, \omega_i)$ is discussed in sec. 2.4, while the impact of the pump spectral shape $\alpha(\omega_s + \omega_i)$ is detailed in sec. 2.2.5.

## 2.2. Mode-selective Single-Photon Subtraction

Fig. 1b illustrates a setup for implementing mode-selective single-photon subtraction from an incoming multimode state $\hat{\rho}$. Similar to the SPA scheme from sec. 2.1, two incident light beams overlap in a $\chi^{(2)}$-medium, resulting in the generation of a third beam. In this configuration, the strong semi-classical beam is labelled as the gate field. The signal state is associated with mode $|a\rangle$, while the up-converted field is associated with mode $|c\rangle$, each characterized by the ladder operators $\{\hat{a}^\dagger, \hat{a}\}$ and $\{\hat{c}^\dagger, \hat{c}\}$, respectively. As photons are created in the up-converted beam, photons are simultaneously annihilated from the signal beam. This evolution can be described by a beam splitter Hamiltonian [7]:

$$ \hat{H} = \iint d\omega_s d\omega_{\text{up}} T(\omega_s, \omega_{\text{up}}) \hat{a}(\omega_s) \hat{c}^\dagger(\omega_{\text{up}}) + \text{h.c.}, \quad (2) $$

where $T(\omega_s, \omega_{\text{up}})$ is the transfer function, analogous to the JSA in the single-photon addition scheme [28]. Therefore, $T(\omega_s, \omega_{\text{up}})$ describes the frequency up-conversion mapping from signal frequencies to output frequencies:

$$ \begin{aligned} T(\omega_s, \omega_{\text{up}}) &= \alpha_g(\omega_{\text{up}} - \omega_s)\phi(\omega_s, \omega_{\text{up}}) \\ &= \sum_l \sqrt{\lambda_l} m_l(\omega_s) n_l(\omega_{\text{up}}) \end{aligned} \quad (3) $$

This can be expressed as the product of the gate beam's spectral shape $\alpha_g(\omega_{\text{up}} - \omega_s)$ and the phase-matching function $\phi(\omega_s, \omega_{\text{up}})$, where $\omega_s$ and $\omega_{\text{up}}$ denote the frequencies of the signal and up-converted
photons, respectively. This decomposition naturally reveals the mode structure of the process, where the eigenmodes {$m_l(\omega_s), n_l(\omega_{up})$} and their corresponding eigenvalues $\lambda_l$ quantify the spectral contributions of individual eigenmode pairs.

Notably, the transfer function and the coupling between the signal and up-converted fields are directly analogous to the quantum pulse gate, first investigated by C. Silberhorn's group [33].

## 2.3. Goal: Desired JSA and TF geometry

To provide a clear schematic representation of the intended processes, we illustrate a target JSA for the SPA process in Fig. 2 and a TF for the SPS process in Fig. 3, each illustrating two different incoming classical modes. In both cases, the advantage of a horizontally oriented phase-matching function $\phi(\omega_i, \omega_k)$ becomes evident. When a strong, coherent beam — the pump in SPA or the gate in SPS — is present, its spectral shape $\alpha(\omega)$ is transferred onto the signal mode, enabling mode-selective heralded operations. This highlights the potential of materials such as thin-film lithium niobate (TFLN), which enable dispersion engineering through geometrical tailoring. When combined with algorithmic optimization, these platforms provide a powerful route toward designing waveguides with desired spectral and quantum-optical properties. In all configurations shown, the idler mode remains confined to a narrow, well-defined spectral region, facilitating a straightforward heralding scheme. However, due to the non-orthogonal geometric relation between the phase-matching function and the pump/gate envelope, the tilted overlap region induces a small deviation from perfect mode purity. This effect remains within an acceptable range, as discussed in Sec. 3.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Schematic JSAs for mode-selective SPA. Here the incoming pump mode (dashed green line) gets transferred onto the signal axis (red signal mode) due to the intersection with the narrow, horizontal phase-matching function (blue area). There, the narrow and horizontal shape of the phase-matching function results (dashed blue line) in a narrow idler mode (red mode on the idler axis), suitable for heralding.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3: Schematic Transfer Function for mode-selective SPS. Here the incoming gate mode gets transferred onto the signal axis due to the intersection with the narrow, horizontal phase-matching function. There, the narrow and horizontal shape of the phase-matching function results in a narrow up-converted mode, suitable for heralding.</figcaption></figure>

## 2.4. Waveguide Simulator

### a. Metallic Waveguide Approximation

To put our optimization approach to the test, we begin by modeling the system using the metallic waveguide approximation described in this section. This model is fully analytically derivable and has been shown to yield good agreement between theoretical predictions and experimental results for dif-fused waveguides [34]. We derive the refractive index and, subsequently, the phase-matching function used to construct the Joint Spectral Amplitudes (JSAs) and Transfer Functions (TFs). This approximation assumes perfectly rectangular waveguides with a discrete refractive index profile, confined strictly to the rectangular geometry of the waveguide (see Fig. 4a). The phase-matching function is approximated as

$$ \phi(\omega) = \text{sinc}\left(\frac{L}{2}\Delta k\right), \qquad (4) $$

where $\Delta k = k_1 + k_2 - k_3 - \frac{2\pi}{\Lambda}$, and $k(\omega) = \frac{n(\omega)\omega}{c}$ is the wave vector, with $n(\omega)$ representing the refractive index as a function of frequency [35]. $k_{1,2,3}$ therefore denote the wave vectors of each field involved in the nonlinear process. $\Lambda$ is the poling period of the magnetic domains in the medium, which can be set during the fabrication process. To get an optimal efficiency, we aim for so-called perfect phase-matching, so $\Delta k = 0$. The poling period is typically chosen to satisfy this condition at the central frequencies of the interacting fields. The frequency-dependent refractive index $n(\omega)$ is calculated using Sellmeier equations specific to the nonlinear material. For the metallic waveguide approximation, we use coefficients for periodically poled lithium niobate (pp:LN) and potassium titanyl phosphate (pp:KTP) as reported in [36, 37].

### b. Thin-Film Model

Secondly, we investigate thin-film lithium niobate (TFLN) waveguides, where the effective refractive indices $n_{eff}(\omega)$ of the guided modes strongly depend on the waveguide geometry. This high sensitivity to structural parameters provides a promising foundation for tailoring waveguides to our desired quantum processes by optimizing their geometrical structure [38, 39]. Here, we consider trapezoidal waveguides as depicted in Fig. 4b, which are characterized by the thin film thickness h, the etching depth g, the waveguide width D and the side wall angle $\phi$. To obtain $n_{eff}$ we solve the wave equation utilizing the
<figure><img src="image_5.png" /><figcaption>Fig. 4: Comparison of two waveguide geometries: (a) a metallic waveguide and (b) a thin-film waveguide. In (a), the cross-section is characterized by the height h and the width w. In (b), the cross-section is characterized by the sidewall angle φ, the waveguide width D, the etching depth g, and the thin-film thickness h. (c) Example values for n<sub>eff</sub> obtained for different TFLN waveguide geometries (besides D the parameters are equal to the ones reported in Table 2a).</figcaption></figure>

commercial software MODE: Waveguide Simulator [40], where we use the refractive indices provided in [41] for lithium niobate and [42] for the SiO2-Substrate.

In the simulation results we identify the fundamental TE and TM modes, by searching for the modes with the highest n<sub>eff</sub> with TE- and TM-fractions above 90%. These n<sub>eff</sub>(ω) profiles of the fundamental modes we then use to obtain the JSAs and TFs. Note, that we ignore geometries in which we identify mode crossings of the fundamental modes (TE-fractions below 90%).

Some example profiles for n<sub>eff</sub> which we obtain for different TFLN waveguide geometries are depicted in Fig. 4c. These highlight the strong tunability of the dispersion properties via changes in the waveguides, which in the following will allow us to design SPA and SPS processes as discussed in Sec. 2.3.

## 2.5. Purity and Mode-Selectivity

For the sake of simplicity we consider the photon subtraction and addition processes on pure Gaussian multimode input states. Such states can always be written in some mode basis as $\hat{\rho} = |\phi_1\rangle\langle\phi_1| \otimes |\phi_2\rangle\langle\phi_2| \otimes \dots \otimes |\phi_n\rangle\langle\phi_n|$. In temporal-spectral multimode Gaussian states, the modes in which the state is factorized is usually a Hermite Gauss basis [34].

We are interested in subtracting or adding photons in a specific mode of the system, which, for simplicity can be just one of the Hermite Gauss modes, on which all the modes are disentangled. For this, it is essential to reduce the number of modes effectively involved in the non-linear interaction as close to 1 as possible. Since the JSA/TF of the process can be decomposed into a set of principal modes, the process itself can be characterized by the number of effective modes, similar to the Schmidt-number [43]:

$$K = \frac{(\sum_n \lambda_n)^2}{\sum_n \lambda_n^2}.$$

This metric is directly related to the output state's purity, where K = 1 corresponds to a pure photon addition or subtraction from a pure input. It should be noted that when we refer to the Schmidt number and to the purity, we refer to the purity of the process itself, and therefore relate to an outcome multimode-state from a pure incoming state.

In case of multimode processes, we can express the output state of the single photon subtraction process as

$$\hat{\rho}_s^{\text{out}} = \sum_l \tilde{\lambda}_l \hat{A}_l \hat{\rho}_s^{\text{in}} \hat{A}_l^\dagger \quad (5)$$
where $\tilde{\lambda}_l = \lambda_l/P^-$ are the normalized squared eigenvalues of the decomposed TF, weighting the contribution of the specific mode with ladder-operators $\{\hat{A}_l, \hat{A}_l^\dagger\}$. The normalization constant, $P^-= \sum_{l,m} \lambda_l O_{lm} n_{lm}$, defines the probability of success of the heralding detection. Here $n_{lm}$ defines the average photon number in the $m^{th}$ input mode, and $O_{lm}$ are the coefficients in the decomposition of the $l^{th}$ eigenmode of the process $(\hat{A}_l)$, in terms of the input modes $(\hat{a}_m)$, $\hat{A}_l = \sum_m O_{lm} \hat{a}_m$. The purity of the output state, which, given that the input state is pure, is directly linked to the mode selectivity of the photon subtraction process is thus given by [see Appendix]:

$$ \mu^- = \mathrm{Tr}(\hat{\rho}_{\mathrm{out}}^2) = \sum_{l,k} \tilde{\lambda}_l \tilde{\lambda}_k \left( \sum_m O_{lm} O_{km} \bar{n}_m \right)^2 . $$

In the case in which $K=1$ and only a single mode is involved in the non-linear process the expression above evaluates exactly to one. It is also possible to notice that if only one of the input modes has non zero average photon number then the purity is also one, as all other modes do not contribute, with the caveat that the output mode gets projected into this single mode. Similarly, the output for a SPA-process can be written as:

$$ \hat{\rho}_s^{\mathrm{out}} = \sum_l \tilde{\lambda}_l \hat{A}_l^\dagger \hat{\rho}_s^{\mathrm{in}} \hat{A}_l \qquad (6) $$

Using the same procedure as above, one can arrive at:

$$ \mu^+ = \sum_{l,k} \tilde{\lambda}_l \tilde{\lambda}_k \left| \sum_m O_{lm} O_{km} (1 + \bar{n}_m) \right|^2, \qquad (7) $$

where $\bar{n}_k$ is the mean photon number in the input mode k. Differently to the case of photon subtraction, here all input modes, including those in vacuum, contribute to the final state. This implies that in the case of photon addition, only by having exactly $K=1$ the process can produce a pure outcome state.

These results provide a framework to understand the connection between the Schmidt eigenvalues and the mode selectivity of the process. The purity $\mu=1$ is attainable for both processes only when the Schmidt decomposition reduces to a single non-zero eigenvalue $\lambda$, corresponding to a single involved mode in the signal and idler fields. This condition requires symmetry along the signal and idler axes in the Joint Spectral Amplitude (JSA), as outlined in [44]. However, in our experimental setup, such symmetry is inherently constrained by the interplay between the waveguide's phase-matching function and the pump spectral profile [see Fig. 7, 8, 11]. Despite this limitation, approximating this symmetry remains desirable for our objectives. While achieving a Schmidt number of $K=1$ is theoretically unattainable, a value close to unity remains our focus. Therefore, we direct our search within the parameter space towards configurations approaching $K \approx 1$. Our primary goal is to design a system that supports a narrow idler-frequency bandwidth for the heralding scheme, while fully exploiting the pump bandwidth in the signal state. To achieve this, we make use of the narrow phase-matching width along the idler axis within our target frequency range, as discussed in Sec. 2.3.

An additional method that can be straightforwardly applied to any outcome of our optimization should be further mentioned — namely, spectral filtering along the idler axis. This can be experimentally realized using a narrowband optical filter or, when fine spectral selectivity is required, a dispersive element, e.g., a diffraction grating, combined with a slit, to isolate a specific spectral window. The filtering is implemented between the photon generation and the heralding process. While it reduces the overall heralding rate by discarding undesired frequency components, it does not degrade the quantum properties of the heralded state or the purity of the process. With this technique, the process can asymptotically approach unit purity while decreasing the event rate of the process [45].
## 2.6. Optimization and Numerical Search

Since we have demonstrated in Sec. 2.5 that achieving a Schmidt number of $K \approx 1$, together with the appropriate shape of the phase-matching function in the LN and KTP waveguides, is sufficient to configure the idler as a narrow-bandwidth heralding mode while transferring the pump bandwidth to the broader signal field, we proceed with a numerical optimization approach, described in Methods. To apply our algorithm to a fully analytical testbed, we model the nonlinear waveguide using the metallic-waveguide approximation, incorporate periodic poling for quasi-phase matching, and determine the refractive index via the Sellmeier equations introduced in Sec. 2.4.a. The phase-matching function is computed based on the waveguide geometry and the type of nonlinear process (Type 0, I, or II), all treated as variable parameters. In addition, we explore the influence of the Hermite-Gauss pump mode order and bandwidth on the resulting JSAs and transfer functions. Using a genetic algorithm, we perform a parameter-space search to minimize the Schmidt number toward $K \approx 1$, thereby identifying waveguide configurations suitable for mode-selective single-photon addition and subtraction.

## 2.7. Testbed Results - Metallic Waveguide approximation

As mentioned in Sec. 2.4.a, we begin by implementing our scheme within the metallic waveguide approximation. Within this framework, we present results that already demonstrate systems capable of performing mode-selective non-Gaussian operations. In cases where differences in mode selectivity were minimal, we opted for the marginally larger pump width to ensure easier experimental implementation in terms of pulse shaping.

### a. Mode-Selective Single-Photon Addition

We start by showing the results for the SPA scheme. We identified an optimized waveguide-pump configuration, with the corresponding parameters listed in Table 1a. These parameters yield the desired JSAs and spectral modes, shown in Figs. 5 and 6. Fig. 5 displays the JSAs for Hermite-Gaussian (HG) pump modes of orders 0, 1, and 2 - all for the same waveguide. Correspondingly, Fig. 6 illustrates the signal and idler modes obtained from the JSA decomposition. The signal modes span over approximately 100 nm, while the idler mode serves as the desired narrow-band mode, enabling an efficient heralding detection scheme. Since the idler mode remains nearly unchanged across all HG orders, we present only the zero-order mode as a representative example for all idler modes.

<figure><img src="image_6.png" /><figcaption>Fig. 5: Joint Spectral Amplitudes (JSAs) corresponding to different pump modes, calculated using the waveguide and pump parameters listed in Table 1a. The involved frequency ranges of the signal and idler modes are directly indicated. Via a diagonalization similar to the Schmidt decomposition, the involved eigenmodes characterizing the purity of the process can be obtained and are shown in Fig. 6.</figcaption></figure>
<table><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Length</td><td>7 mm</td></tr><tr><td>Width</td><td>2.8 µm</td></tr><tr><td>Height</td><td>2.3 µm</td></tr><tr><td>Material</td><td>KTP</td></tr><tr><td>Pump width</td><td>7 nm</td></tr><tr><td>Pump mode</td><td>HG 0/1/2</td></tr><tr><td>Type</td><td>II</td></tr><tr><td>Polarization (signal)</td><td>H</td></tr><tr><td>Polarization (idler)</td><td>V</td></tr><tr><td>Polarization (pump)</td><td>H</td></tr></tbody></table>

(a) Waveguide and pump parameters used to generate the JSAs for the SPA process and the associated modes shown in Figs. 5 and 6.

<table><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Length</td><td>2.0 mm</td></tr><tr><td>Width</td><td>2.5 µm</td></tr><tr><td>Height</td><td>1.8 µm</td></tr><tr><td>Material</td><td>KTP</td></tr><tr><td>Pump width</td><td>10 nm</td></tr><tr><td>Pump mode</td><td>HG 0/1/2</td></tr><tr><td>Type</td><td>0</td></tr><tr><td>Polarization (signal)</td><td>V</td></tr><tr><td>Polarization (idler)</td><td>V</td></tr><tr><td>Polarization (pump)</td><td>V</td></tr></tbody></table>

(b) Waveguide and pump parameters used to generate the TFs for the SPS process and the associated modes shown in Figs. 7 and 8.

<figure><img src="image_7.png" /><figcaption>Fig. 6: Signal and idler modes corresponding to the JSAs shown in Fig. 5, obtained using the parameters listed in Table 1a.</figcaption></figure>

b. Mode-Selective Single-Photon Subtraction

We further could identify a sufficient waveguide-pump-configuration for the SPS-process. The parameters are listed in Table 1b and yield the desired transfer functions (TFs) and spectral modes, shown in Figs. 7 and 8. Fig. 7 displays the TFs for Hermite-Gaussian gate modes of orders 0, 1, and 2. Correspondingly, Fig. 8 illustrates the signal and up-converted modes obtained from the TF decomposition. The signal modes span approximately 110 nm, while the up-converted mode serves as the desired narrow-band mode, enabling an efficient heralding detection scheme. Since the gate mode remains nearly unchanged across all HG orders, we present only the zero-order mode as a representative example for all up-converted modes.
<figure><img src="image_8.png" /><figcaption>Fig. 7: Transfer functions (TFs) corresponding to the pump modes HG₀, HG₁, and HG₂, calculated using the parameters listed in Table 1b. The involved frequency ranges of the signal and up-converted modes are directly indicated.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Signal and gate modes corresponding to the TFs shown in Fig. 7, obtained using the parameters listed in Table 1b.</figcaption></figure>

### 3. Results - Thin-Film Lithium Niobate waveguides

Since we could apply our optimization scheme to the analytical model, we now present our implementation to the TFLN-model (see Sec. 2.4.b), showing how finely tuned the dispersion engineering can be done to achieve the desired non-gaussian processes.

#### 3.1. Single-Photon Addition

Starting from the thin-film model, we identified an optimized waveguide-pump configuration for mode-selective single-photon addition. The respective waveguide- and pump-parameters are listed
<table><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Length</td><td>7 mm</td></tr><tr><td>Width / D</td><td>608 nm</td></tr><tr><td>Height / h</td><td>600 nm</td></tr><tr><td>Etching angle / φ</td><td>78.2°</td></tr><tr><td>Etching depth / g</td><td>441 nm</td></tr><tr><td>Material</td><td>LN</td></tr><tr><td>Pump width</td><td>6 mm</td></tr><tr><td>Pump mode</td><td>HG 0/1/2</td></tr><tr><td>Type</td><td>II</td></tr><tr><td>Polarization (signal)</td><td>V</td></tr><tr><td>Polarization (idler)</td><td>H</td></tr><tr><td>Polarization (pump)</td><td>V</td></tr></tbody></table>

(a) Waveguide and pump parameters used to generate the JSAs for the SPA process and the associated modes shown in Figs. 9 and 10.

<table><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Length</td><td>7 mm</td></tr><tr><td>Width / D</td><td>883 nm</td></tr><tr><td>Height / h</td><td>788 nm</td></tr><tr><td>Etching angle / φ</td><td>75.2°</td></tr><tr><td>Etching depth / g</td><td>774 nm</td></tr><tr><td>Material</td><td>LN</td></tr><tr><td>Pump width</td><td>7 mm</td></tr><tr><td>Pump mode</td><td>HG 0/1/2</td></tr><tr><td>Type</td><td>II</td></tr><tr><td>Polarization (signal)</td><td>H</td></tr><tr><td>Polarization (idler)</td><td>V</td></tr><tr><td>Polarization (pump)</td><td>V</td></tr></tbody></table>

(b) Waveguide and pump parameters used to generate the TFs for the SPS process and the associated modes shown in Figs. 11 and 12.

<caption>Tab. 2: Comparison of waveguide and pump parameters for thin-film configurations enabling single-photon subtraction and addition. The parameters width (D), height (h), etching angle (φ), and etching depth (g) are illustrated in Fig. 4b.</caption>

in Table 2a. The obtained Joint Spectral Amplitudes (JSAs) for HG pump modes of orders 0, 1, and 2 are shown in Fig. 9, while the corresponding signal and idler modes, extracted from the JSA decomposition, are presented in Fig. 10. Since the idler mode remains nearly unchanged across all HG orders, we present only the zero-order signal mode as a representative example.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Joint-Spectral Amplitudes (JSA) for HG₀, HG₁, and HG₂, corresponding to the parameters in Table 2a. As in the metallic approximation, the involved quantum modes can be obtained by a Schmidt-decomposition and are presented in Fig. 10.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) Obtained signal mode with pump: HG₀</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Obtained signal mode with pump: HG₁</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Obtained signal mode with pump: HG₂</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Obtained idler mode for HG₀-pump. Representative for all gate modes.</figcaption></figure>

Fig. 10: Signal and idler modes corresponding to the JSAs shown in Fig. 9, obtained using the parameters listed in Table 2a.

## 3.2. Single-Photon Subtraction

Following the same approach, we now present the results for single-photon subtraction. The optimized waveguide and pump parameters yield the desired transfer functions (TFs) for HG pump modes of orders 0, 1, and 2, as shown in Fig. 11. The corresponding signal and up-converted modes are presented in Fig. 12. Since the up-converted mode remains nearly unchanged across all incoming HG orders, we present only the zero-order signal mode as a representative example.

<figure><img src="image_15.png" /><figcaption>(a) Pump mode: HG₀</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Pump mode: HG₁</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Pump mode: HG₂</figcaption></figure>

Fig. 11: Transfer functions (TFs) for HG₀, HG₁, and HG₂, corresponding to the parameters in Table 2b. The decomposed signal and up-converted modes are shown in Fig. 12.
<figure><img src="image_18.png" /><figcaption>Fig. 12: Signal and up-converted modes corresponding to the TFs shown in Fig. 11, obtained using the parameters listed in Table 2b.</figcaption></figure>

4. Conclusion

In this work, we investigated mode-selective single-photon addition and subtraction using two complementary waveguide modeling approaches. Starting from an analytically derivable metallic waveguide approximation—representative of indiffused waveguides on the micrometer scale—we extended our optimization scheme to a self-developed model of subwavelength-scale thin-film lithium niobate (TFLN) waveguides. The strong potential for dispersion engineering through geometric design in this platform enabled us to directly propose configurations suitable for implementing SPA and SPS processes.

By optimizing the waveguide–pump parameters in each case, we identified parameter sets yielding the desired Joint Spectral Amplitudes and Transfer Functions, thus enabling mode-selective single-photon addition and subtraction and facilitating the generation of non-Gaussian quantum states. For both processes, we found configurations that can serve as standalone experimental proposals. In the case of the fundamental Hermite-Gauss pump mode ($HG_0$), Schmidt numbers below $K = 1.07$ were found for SPA using both the metallic and TFLN waveguide models. For SPS, the optimization yielded Schmidt numbers of approximately $K = 1.04$ for the metallic model and $K = 1.03$ for TFLN. In all cases, higher-order pump modes ($HG_1, HG_2$) led to similarly satisfactory results, confirming the robustness of the identified solutions with respect to mode-selectivity.

These findings demonstrate that, beyond the well-established metallic waveguide framework, the TFLN platform also constitutes a highly promising source for non-Gaussian state generation. Both platforms enable the selective generation of non-Gaussian states in a specific mode within a multimode input, using a single, static device configuration. The implementation is based on heralding via photon detection, which, to date, is the most efficient procedure for generating complex states in continuous variables [46–56]. Here, we control the spectro-temporal mode where such an operation is applied [15, 28–30, 57], which is particularly relevant since full scalability requires involving a large number of spectro-temporal modes of the field [58]. Moreover, this specific implementation preserves coherence among multiple modes and enables the exploitation of the interplay between non-Gaussianity and quantum entanglement [59, 60].
## Funding

This work was supported by the European Research Council under the Consolidator Grant COQ-COoN (Grant No. 820079), by Agence Nationale de la Recherche (OQuLus, ANR-22-PETQ-0013) and by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – SFB-Geschäftszeichen TRR142/3-2022 – Projektnummer 231447078. P. N. was supported by QuantEdu-France (ANR-22-CMAS-0001) in the framework of France 2030.

## Disclosure

The authors declare no conflicts of interest.

## Data availabilty

Data and code underlying the results presented in this work are not publicly available at this time but may be obtained from the authors upon reasonable request.

[1] C. Fabre and N. Treps, Rev. Mod. Phys. 92, 035005 (2020), URL https://link.aps.org/doi/10.1103/RevModPhys.92.035005.

[2] B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, Phys. Rev. X 5, 041017 (2015), URL https://link.aps.org/doi/10.1103/PhysRevX.5.041017.

[3] M. V. Larsen, X. Guo, C. R. Breum, J. S. Neergaard-Nielsen, and U. L. Andersen, Science 366, 369 (2019), ISSN 0036-8075, URL https://science.sciencemag.org/content/366/6463/369.

[4] A. Inoue, T. Kashiwazaki, T. Yamashima, N. Takanashi, T. Kazama, K. Enbutsu, K. Watanabe, T. Umeki, M. Endo, and A. Furusawa, Applied Physics Letters 122, 104001 (2023), https://doi.org/10.1063/5.0137641, URL https://doi.org/10.1063/5.0137641.

[5] T. Kashiwazaki, T. Yamashima, K. Enbutsu, T. Kazama, A. Inoue, K. Fukui, M. Endo, T. Umeki, and A. Furusawa, Applied Physics Letters 122, 234003 (2023), ISSN 0003-6951, https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/5.0144385/17972007/234003_1_5.0144385.pdf, URL https://doi.org/10.1063/5.0144385.

[6] R. Nehra, R. Sekine, L. Ledezma, Q. Guo, R. M. Gray, A. Roy, and A. Marandi, Science 377, 1333 (2022), https://www.science.org/doi/pdf/10.1126/science.abo6213, URL https://www.science.org/doi/abs/10.1126/science.abo6213.

[7] M. Walschaers, PRX Quantum 2, 030204 (2021), URL https://link.aps.org/doi/10.1103/PRXQuantum.2.030204.

[8] D. Zhu, L. Shao, M. Yu, R. Cheng, B. Desiatov, C. J. Xin, Y. Hu, J. Holzgrafe, S. Ghosh, A. Shams-Ansari, et al., Adv. Opt. Photon. 13, 242 (2021), URL https://opg.optica.org/aop/abstract.cfm?URI=aop-13-2-242.

[9] A. Boes, L. Chang, C. Langrock, M. Yu, M. Zhang, Q. Lin, M. Lončar, M. Fejer, J. Bowers, and A. Mitchell, Science 379, eabj4396 (2023).

[10] L. Wang, H. Du, X. Zhang, and F. Chen, Journal of Physics D: Applied Physics 58, 023001 (2024), URL https://dx.doi.org/10.1088/1361-6463/ad7ff7.

[11] A. Zavatta, S. Viciani, and M. Bellini, Science 306, 660 (2004).

[12] A. Zavatta, S. Viciani, and M. Bellini, Phys. Rev. A 72, 023820 (2005), URL https://link.aps.org/doi/10.1103/PhysRevA.72.023820.

[13] A. Zavatta, V. Parigi, and M. Bellini, Phys. Rev. A 75, 052106 (2007), URL https://link.aps.org/doi/10.1103/PhysRevA.75.052106.

[14] V. Parigi, A. Zavatta, M. Kim, and M. Bellini, Science 317, 1890 (2007), ISSN 00368075, 10959203, URL http://www.jstor.org/stable/20048471.

[15] Y.-S. Ra, C. Jacquard, A. Dufour, C. Fabre, and N. Treps, Phys. Rev. X 7, 031012 (2017), URL https://link.aps.org/doi/10.1103/PhysRevX.7.031012.

[16] J. O'Brien, A. Furusawa, and J. Vučković, Nature Photonics 3, 687 (2009), ISSN 1749-4885.
[17] Y. Qi and Y. Li, Nanophotonics 9, 1287 (2020), URL https://doi.org/10.1515/nanoph-2020-0013.

[18] M. G. Vazimali and S. Fathpour, Advanced Photonics 4, 034001 (2022), URL https://doi.org/10.1117/1.AP.4.3.034001.

[19] Y. Koyaz, C. Lafforgue, H. Zarebidaki, O. Hefti, D. Grassani, H. Sattari, and C.-S. Brès, Opt. Express 32, 46776 (2024), URL https://opg.optica.org/oe/abstract.cfm?URI=oe-32-26-46776.

[20] M. Arora, J. Arya, P. Chokkara, and J. Lugani, arXiv preprint arXiv:2411.17369 (2024), accessed: 2025-02-21, URL https://arxiv.org/abs/2411.17369.

[21] R. Luo, Y. He, H. Liang, M. Li, and Q. Lin, Optica 5, 1006 (2018), URL https://opg.optica.org/optica/abstract.cfm?URI=optica-5-8-1006.

[22] C. J. Xin, S. Lu, J. Yang, A. Shams-Ansari, B. Desiatov, L. S. Magalhães, S. S. Ghosh, E. McGee, D. Renaud, N. Achuthan, et al., Communications Physics 8, 136 (2025), URL https://doi.org/10.1038/s42005-025-02068-3.

[23] I. Krasnokutska, J.-L. J. Tambasco, X. Li, and A. Peruzzo, Opt. Express 26, 897 (2018), URL https://opg.optica.org/oe/abstract.cfm?URI=oe-26-2-897.

[24] A. Shams-Ansari, G. Huang, L. He, Z. Li, J. Holzgrafe, M. Jankowski, M. Churaev, P. Kharel, R. Cheng, D. Zhu, et al., APL Photonics 7, 081301 (2022), ISSN 2378-0967, https://pubs.aip.org/aip/app/article-pdf/doi/10.1063/5.0095146/16492670/081301_1_online.pdf, URL https://doi.org/10.1063/5.0095146.

[25] X. Zhu, Y. Hu, S. Lu, H. K. Warner, X. Li, Y. Song, L. Magalhães, A. Shams-Ansari, A. Cordaro, N. Sinclair, et al., Photonics Research 12, A63 (2024), ISSN 2327-9125, publisher: Optica Publishing Group, URL https://opg.optica.org/abstract.cfm?URI=prj-12-8-A63.

[26] J. Kellner, A. Sabatti, T. Kuttner, R. J. Chapman, and R. Grange, *Counter-propagating spontaneous parametric down-conversion source in lithium niobate on insulator* (2025), 2506.21396, URL https://arxiv.org/abs/2506.21396.

[27] A. Sabatti, J. Kellner, R. J. Chapman, and R. Grange, *Nanodomain poling unlocking backward nonlinear light generation in thin film lithium niobate* (2025), 2507.13004, URL https://arxiv.org/abs/2507.13004.

[28] V. A. Averchenko, V. Thiel, and N. Treps, Phys. Rev. A 89, 063808 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.89.063808.

[29] V. Averchenko, C. Jacquard, V. Thiel, C. Fabre, and N. Treps, New Journal of Physics 18, 083042 (2016), URL https://dx.doi.org/10.1088/1367-2630/18/8/083042.

[30] G. Roeland, S. Kaali, V. R. Rodriguez, N. Treps, and V. Parigi, New Journal of Physics 24, 043031 (2022), URL https://dx.doi.org/10.1088/1367-2630/ac5f85.

[31] V. Roman-Rodriguez, B. Brecht, S. K, C. Silberhorn, N. Treps, E. Diamanti, and V. Parigi, New Journal of Physics 23, 043012 (2021), ISSN 1367-2630, URL http://dx.doi.org/10.1088/1367-2630/abef96.

[32] C. K. Law, I. A. Walmsley, and J. H. Eberly, Phys. Rev. Lett. 84, 5304 (2000), URL https://link.aps.org/doi/10.1103/PhysRevLett.84.5304.

[33] A. Eckstein, B. Brecht, and C. Silberhorn, Opt. Express 19, 13770 (2011), URL https://opg.optica.org/oe/abstract.cfm?URI=oe-19-15-13770.

[34] V. Roman-Rodriguez, D. Fainsin, G. L. Zanin, N. Treps, E. Diamanti, and V. Parigi, Physical Review Research 6, 043113 (2024).

[35] R. W. Boyd, *Nonlinear Optics* (Academic Press, San Diego, 2008), 3rd ed., ISBN 978-0-12-369470-6.

[36] G. J. Edwards and M. Lawrence, Opt. Quantum Electron. 16, 373 (1984).

[37] D. H. Jundt, Opt. Lett. 22, 1553 (1997), URL https://opg.optica.org/ol/abstract.cfm?URI=ol-22-20-1553.

[38] M. Jankowski, J. Mishra, and M. M. Fejer, Journal of Physics: Photonics 3, 042005 (2021), ISSN 2515-7647, publisher: IOP Publishing, URL https://iopscience.iop.org/article/10.1088/2515-7647/ac1729.

[39] C. J. Xin, J. Mishra, C. Chen, D. Zhu, A. Shams-Ansari, C. Langrock, N. Sinclair, F. N. C. Wong, M. M. Fejer, and M. Lončar, Optics Letters 47, 2830 (2022), ISSN 0146-9592, 1539-4794, publisher: Optica Publishing Group, URL https://opg.optica.org/abstract.cfm?URI=ol-47-11-2830.

[40] Lumerical Inc., *Lumerical: Photonic design software* (2023), URL https://www.lumerical.com/.

[41] O. Gayer, Z. Sacks, E. Galun, and A. Arie, Applied Physics B 91, 343 (2008).

[42] E. D. Palik, ed., *Handbook of Optical Constants of Solids* (Academic Press, 1985).

[43] A. Ekert and P. L. Knight, American Journal of Physics 63, 415 (1995).

[44] W. P. Grice, A. B. U'Ren, and I. A. Walmsley, Phys. Rev. A 64, 063815 (2001), URL https://link.aps.org/doi/10.1103/PhysRevA.64.063815.

[45] M. G. Raymer, J. Noh, K. Banaszek, and I. A. Walmsley, Physical Review A 72, 023825 (2005), quant-ph/0504062, URL https://journals.aps.org/pra/abstract/10.1103/PhysRevA.72.023825.

[46] M. V. Larsen, J. E. Bourassa, S. Kocsis, J. F. Tasker, R. S. Chadwick, C. González-Arciniegas, J. Hastrup,
C. E. Lopetegui-González, F. M. Miatto, A. Motamedi, et al., Nature (2025), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-025-09044-5.

[47] M. S. Winnel, J. J. Guanzon, D. Singh, and T. C. Ralph, Phys. Rev. Lett. 132, 230602 (2024), URL https://link.aps.org/doi/10.1103/PhysRevLett.132.230602.

[48] S. Konno, W. Asavanant, F. Hanamura, H. Nagayoshi, K. Fukui, A. Sakaguchi, R. Ide, F. China, M. Yabuno, S. Miki, et al., Science 383, 289 (2024), https://www.science.org/doi/pdf/10.1126/science.adk7560, URL https://www.science.org/doi/ abs/10.1126/science.adk7560.

[49] A. Anteneh, L. Brunel, and O. Pfister, Optica Quantum 2, 296 (2024), URL https://opg.optica.org/ opticaq/abstract.cfm?URI=opticaq-2-4-296.

[50] H. Simon, L. Caron, R. Journet, V. Cotte, and R. Tualle-Brouri, Phys. Rev. Lett. 133, 173603 (2024), URL https://link.aps.org/doi/10.1103/PhysRevLett.133.173603.

[51] D. V. Sychev17, A. E. Ulanov, A. A. Pushkina, M. W. Richards, I. A. Fedorov, and A. I. Lvovsky, Nature Photonics 11, 379 (2017), URL https://doi.org/10.1038/nphoton.2017.57.

[52] J. Etesse, M. Bouillard, B. Kanseri, and R. Tualle-Brouri, Phys. Rev. Lett. 114, 193602 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.114.193602.

[53] T. Gerrits, S. Glancy, T. S. Clement, B. Calkins, A. E. Lita, A. J. Miller, A. L. Migdall, S. W. Nam, R. P. Mirin, and E. Knill, Phys. Rev. A 82, 031802 (2010), URL https://link.aps.org/doi/10.1103/ PhysRevA.82.031802.

[54] H. Takahashi, K. Wakui, S. Suzuki, M. Takeoka, K. Hayasaka, A. Furusawa, and M. Sasaki, Phys. Rev. Lett. 101, 233605 (2008), URL https://link.aps.org/doi/10.1103/PhysRevLett.101.233605.

[55] A. Ourjoumtsev, H. Jeong, R. Tualle-Brouri, and P. Grangier, Nature 448, 784 (2007), URL https://doi.org/10.1038/nature06054.

[56] J. S. Neergaard-Nielsen, B. M. Nielsen, C. Hettich, K. Mølmer, and E. S. Polzik, Phys. Rev. Lett. 97, 083604 (2006), URL https://link.aps.org/doi/10.1103/PhysRevLett.97.083604.

[57] N. Biagi, S. Francesconi, M. Gessner, M. Bellini, and A. Zavatta, Advanced Quantum Technologies 5, 2200039 (2022).

[58] J. E. Bourassa, R. N. Alexander, M. Vasmer, A. Patil, I. Tzitrin, T. Matsuura, D. Su, B. Q. Baragiola, S. Guha, G. Dauphinais, et al., Quantum 5, 392 (2021), ISSN 2521-327X, URL https://doi.org/10. 22331/q-2021-02-04-392.

[59] Y.-S. Ra, A. Dufour, M. Walschaers, C. Jacquard, T. Michel, C. Fabre, and N. Treps, Nature Physics 16, 144 (2020), URL https://doi.org/10.1038/s41567-019-0726-y.

[60] C. E. Lopetegui, M. Isoard, N. Treps, and M. Walschaers, Optica Quantum 3, 312 (2025), URL https://opg.optica.org/opticaq/abstract.cfm?URI=opticaq-3-4-312.

APPENDIX

Analytics of multimode single-photon subtraction and addition

Following the approaches developed in the works of V. Averchenko [28, 29] and G. Roeland [30], we present a unified description of the mathematical frameworks to describe the processes of single-photon subtraction (SPS) and single-photon addition (SPA). This formalism provides a consistent basis for analyzing these quantum operations in a beam-splitter-based scheme as well as a PDC scheme.

Modeling Single-Photon Subtraction (SPS) as a Weak Beam Splitter

To begin, we consider a pure, single mode, SPS process, to then treat the more general case. Such a process can be modelled by an input state undergoing a weak beam-splitter operation. This can be described by the evolution $\hat{U}$:

$$ \hat{U} = \exp(i\theta(\hat{A}\hat{B}^{\dagger} - \hat{A}^{\dagger}\hat{B})) \approx 1 + i\theta(\hat{B}^{\dagger}\hat{A} - \hat{B}\hat{A}^{\dagger}), $$

where a weak coupling regime is assumed for the right hand side approximation, implying that only single photon exchanges between signal and idler take place. The initial state in our SPS-scheme is given by $\hat{\rho}_{in} = \hat{\rho}_{s} \otimes |0\rangle_{B}\langle0|$, where $\hat{\rho}_{s}$ is the state of the input signal, while the idler mode is on the
vacuum. To perform the single photon subtraction from the signal mode, post-selection on a single photon going to the idler is required. The detection event is described by the POVM $\Pi_B = |1\rangle\langle 1|_B$, where $|1\rangle_B$, is the single photon Fock state in the idler mode. Ideally, a photon number resolving detector should be used, but in practice click detectors suffice in the weak coupling regime that we assume. The appropriate POVM for such a detector is $\Pi_B^{\text{click}} = 1 - |0\rangle\langle 0|_B$, which, given that only single photon events play a significant role, effectively acts as $\Pi_B^{\text{click}} = \Pi_B$.

The output state after a successful heralding of in the up-converted mode is given by

$$ \hat{\rho}_{out} = \frac{1}{P} \text{Tr}_B \left( \underbrace{\hat{U} \hat{\rho}_{in} \hat{U}^{\dagger} \Pi_B}_{\rho'_{out}} \right), $$

where $P$ describes the probability of successfully measuring a heralding photon, and thus subtracting it from the signal. The state before heralding is given by $\rho'_{out}$:

$$ \hat{\rho}'_{out} = \hat{\rho}_s \otimes |0\rangle_B \langle 0| + i\theta(\hat{A}\hat{\rho}_s \otimes |1\rangle_B \langle 0|) - i\theta(\hat{\rho}_s \hat{A}^{\dagger} \otimes |0\rangle_B \langle 1|) + \theta^2 \hat{A}\hat{\rho}_s \hat{A}^{\dagger} \otimes |1\rangle_B \langle 1|. $$

Therefore, for a single-photon subtracted state we obtain:

$$ \hat{\rho}_{out} = \frac{1}{P} \text{Tr}_B(\hat{\rho}_{out} \Pi_B) = \frac{\theta^2}{P} \hat{A} \hat{\rho}_s \hat{A}^{\dagger} $$

with $P = \theta^2 \text{Tr}[\hat{A}^{\dagger}\hat{A}\hat{\rho}_s] = \theta^2 \bar{n}_s$, the probability of the process to detect a photon in the heralding path, and $\bar{n}_s$ the average photon number of the input state. Notice that there are two ways to boost the probability of success: increasing the average photon number of the input (equivalently, the amount of squeezing in the target mode), or having a stronger coupling (higher beam splitter transmissivity, or stronger pump and non-linearity in the case of up conversion process). Nevertheless, this cannot be done arbitrarily as otherwise higher photon contributions may appear. The non-linear term is effectively not in the equation for the output state as it appears as a multiplicative factor, meaning we can write the state as $\hat{\rho}_{out} = \frac{1}{N_A} \hat{A}\hat{\rho}_s\hat{A}^{\dagger}$, where $N_A = \bar{n}_s$.

## Purity of the SPS-process

The experimental implementation of single photon subtraction is, on the other hand, a non-pure process. While the beam splitter approach provides an appealing way of deriving the effect of single photon subtraction, the non-linear process through which we implement this process is more subtle. Given that the photon detection scheme is not mode selective, and several signal modes are coupled through the non-linear process to the idler, a successful heralding event is effectively a statistical mixture of subtracted states. In what follows we characterize the purity of these states.

The purity $\mu^-$ of a state is defined as $\mu^- = \text{Tr}[\hat{\rho}_{out}^2]$. Since in this work we assume pure multimode states as an input, we can characterize the purity of the output states by simultaneously describing the purity in mode-selectivity of the nonlinear process per se.

We start by considering an input multimode state $\hat{\rho}_{in} = |\phi\rangle\langle\phi| = \bigotimes_{i=1}^n |\phi_i\rangle\langle\phi_i|$. In all simulated scenarios the Schmidt decomposition of the JSAs/TFs results in several non-zero eigenvalues, which means that the process not only involves our target mode, but also other modes or thogonal to it. The Hamiltonian of the non-linear process can be written as

$$ \hat{H} = \iint d\omega_s d\omega_{up} \hat{a}(\omega_s) \hat{b}^{\dagger}(\omega_{up}) + \text{h.c.} $$

After performing a Singular value decomposition on the transfer function,

$$ T(\omega_s, \omega_{up}) = \sum_l \sqrt{\lambda_l} m_l(\omega_s) n_l(\omega_{up}), $$
the process can be written as

$$ \hat{H} = \sum_{l} \sqrt{\lambda_{l}} (\hat{A}_{l} \hat{B}_{l}^{\dagger} + \text{h.c.}), $$

where $\hat{A}_l = \int d\omega_s m_l(\omega_s) \hat{a}(\omega_s)$, and $\hat{B}_l = \int d\omega_{up} n_l(\omega_{up}) \hat{b}(\omega_{up})$, are the annihilation operators of the eigenmodes of the process $\{m_l(\omega_s)\}, \{n_l(\omega_{up})\}$. Following a calculation similar to the one performed in the single mode case, we obtain the output state

$$ \hat{\rho}_{out} = \frac{1}{P} \sum_{l} \lambda_{l} A_{l} \hat{\rho}_{s} A_{l}^{\dagger}, $$

where $P = \sum_l \lambda_l \tilde{n}_l$, is the success probability of the photon subtraction. Notice that, as discussed on the introduction of this section, the output state is given by a mixture of photon subtractions in different modes, due to the unavoidable multimode nature of the process and the non mode-selective photon counting. Nevertheless, for most cases of interest, we obtain a value of $\lambda_1$ significantly higher than the rest, meaning that the contribution from subtractions in undesired modes is minimal. We recall here that we consider the input state to be a pure state $\hat{\rho}_s = |\phi_s\rangle\langle\phi_s|$. This allows us to use, as a figure of merit for the mode selectivity of the process, the purity of the output state, which is given by

$$ \mathrm{Tr}[\hat{\rho}_{out}^2] = \frac{1}{P^2} \left( \sum_l \lambda_l^2 \tilde{n}_l^2 + \sum_{l \neq k} \lambda_l \lambda_k |\langle\phi_s|\hat{A}_l^\dagger \hat{A}_k|\phi_s\rangle|^2 \right). $$

Defining the normalized coefficients:

$$ \tilde{\lambda}_i = \frac{\lambda_i}{P}, $$

we obtain:

$$ \mu^{-} = \mathrm{Tr}(\hat{\rho}_{out}^{2}) = \sum_{l,k} \tilde{\lambda}_{l} \tilde{\lambda}_{k} |\langle\phi_{s}|\hat{A}_{l}^{\dagger}\hat{A}_{k}|\phi_{s}\rangle|^{2}. $$

To figure out a neat form for this evaluation we recall that the input state, $|\phi_s\rangle$, is given by

$$ |\phi_s\rangle = \bigotimes_i |\phi_i\rangle, $$

where

$$ |\phi_i\rangle = \frac{1}{\sqrt{\cosh r_i}} \sum_n \frac{\sqrt{(2n)!}}{2^n n!} (e^{-i\xi_i} \tanh(r_i))^n |2n\rangle, $$

where $r_i, \xi_i$ are the squeezing level and phase respectively. The basis in which the input state is written as a tensor product $\{\phi_i(\omega_s)\}$, is not necessarily the same as the SVD basis of the photon subtraction process, $\{m_l(\omega_s)\}$. Because of this, contributions from crossed terms on the purity cannot be dismissed. The eigenmodes of the photon subtraction process are related to the squeezed modes by $O_{ij} = \int d\omega_s m_i(\omega_s) \phi_j(\omega_s)$, such that $\hat{A}_l = \sum_n O_{ln} \hat{a}_n$, where $\hat{a}_n$ is the annihilation operator corresponding to the $n^{th}$ squeezed state. As a consequence,

$$ \mu^{-} = \sum_{l,k} \tilde{\lambda}_{l} \tilde{\lambda}_{k} \left| \sum_{m,n} O_{lm} O_{kn} \langle \phi_{s} | \hat{a}_{m}^{\dagger} \hat{a}_{n} | \phi_{s} \rangle \right|^{2}, $$
which, taking into account taht for vacuum squeezed input states $\langle\phi_i|\hat{a}_i|\phi_i\rangle = 0$, implies that

$$ \mu^{-} = \sum_{l,k} \tilde{\lambda}_{l} \tilde{\lambda}_{k} \left( \sum_{m} O_{lm} O_{km} \bar{n}_{m} \right)^{2}. $$

On the other hand, if the states are bright squeezed state, i.e. if they have non zero mean field, then the overlap $\langle\phi_i|\hat{a}_i|\phi_i\rangle \neq 0$ and more contributions have to be considered.

## Modeling Single-Photon Addition (SPA) as a PDC-process

While the SPS process is modeled as a weak beam-splitter interaction, the SPA process can be described by a parametric down-conversion (PDC) process, as sketched in Fig. 1a. The Hamiltonian of the process is given by

$$ \hat{H} = \iint d\omega_s d\omega_i J(\omega_s, \omega_i) \hat{a}^{\dagger}(\omega_s) \hat{b}^{\dagger}(\omega_i) + \text{h.c.} $$

Doing SVD on the JSA $J(\omega_s, \omega_i) = \sum_l \sqrt{\lambda_l} h_l(\omega_s) g_l(\omega_i)$, we obtain the eigenmode decomposition of the process

$$ \hat{H} = \sum_{l} \sqrt{\lambda_{l}} \hat{A}_{l}^{\dagger} \hat{B}_{l}^{\dagger} + \text{h.c.}, $$

where $\hat{A}_l^\dagger = \int d\omega_s h_l(\omega_s) \hat{a}(\omega_s)$, and $\hat{B}_l^\dagger = \int d\omega_s g_l(\omega_s) \hat{b}(\omega_s)$, are the creation operator of the eigenmodes of the process. Under similar assumptions, and following a similar procedure as for the single photon subtraction case, we can show that the output state is given by

$$ \hat{\rho}_{\text{out}}^{+} = \frac{1}{P} \sum_{l} \lambda_{l} \hat{A}_{l}^{\dagger} \hat{\rho} \hat{A}_{l}, $$

where P is the probability of detecting a photon in the idler mode:

$$ P = \text{Tr}(\sum_{l} \lambda_{l} \hat{A}_{l}^{\dagger} \hat{\rho} \hat{A}_{l}) = \sum_{l} \lambda_{l} (1 + \langle \phi_{s} | \hat{A}_{l}^{\dagger} \hat{A}_{l} | \phi_{s} \rangle), $$

which, considering the expansion of the eigenmodes of the process in terms of the signal (Hermite Gauss) modes $\hat{A}_l = \sum_m O_{lm} \hat{a}_m$, and the fact that each of the signal modes is in a squeezed vacuum state, leads to

$$ P = \sum_{l} \lambda_{l} + \sum_{l,m} \lambda_{l} O_{lm}^{2} \bar{n}_{m}. $$

where $\bar{n}_m$ denotes the mean photon number in the $m^{th}$ signal (Hermite Gauss) mode.

The purity of the outcoming state $\mu^+ = \text{Tr}[\hat{\rho}_{out}^2]$, is thus given by:

$$ \mu^{+} = \sum_{l,k} \tilde{\lambda}_{l} \tilde{\lambda}_{k} \left| \sum_{m} O_{lm} O_{km} (1 + \bar{n}_{m}) \right|^{2}, $$

where $\tilde{\lambda}_i = \frac{\lambda_i}{P}$, are normalized coefficients.
Purity Evaluation of the Obtained Processes

The previously derived expressions allow us to directly quantify the intrinsic purity of the obtained processes. One should keep in mind, as discussed in Sec. 2.5, that additional spectral filtering on the idler/up-converted field can increase the purity further towards unity. As highlighted in Sec. 2.7, the spectral width of the pump was chosen to strike a "sweet spot" between spectral shapability and a process Schmidt number of $K \approx 1$. For us, a sufficient configuration was found to correspond to a pump width of 7 nm for the metallic waveguide approximation and 10 nm for the TFLN model in the SPS process. For the SPA process, 7 nm was chosen for the metallic waveguide approximation and 6 nm for the TFLN model.

<table><caption>Tab. 3: SPS process</caption><thead><tr><th></th><th>μ<sup>-</sup><sub>met</sub></th><th>μ<sup>-</sup><sub>TFLN</sub></th></tr></thead><tbody><tr><td>HG<sub>0</sub></td><td>0.969</td><td>0.976</td></tr><tr><td>HG<sub>1</sub></td><td>0.942</td><td>0.941</td></tr><tr><td>HG<sub>2</sub></td><td>0.919</td><td>0.909</td></tr></tbody></table>

<table><caption>Tab. 4: SPA process</caption><thead><tr><th></th><th>μ<sup>+</sup><sub>met</sub></th><th>μ<sup>+</sup><sub>TFLN</sub></th></tr></thead><tbody><tr><td>HG<sub>0</sub></td><td>0.952</td><td>0.954</td></tr><tr><td>HG<sub>1</sub></td><td>0.884</td><td>0.895</td></tr><tr><td>HG<sub>2</sub></td><td>0.826</td><td>0.822</td></tr></tbody></table>

Tab. 5: Process purities for the first three Hermite-Gauss pump/gate modes, containing the metallic waveguide approximation and the TFLN model for (a) the SPS process and (b) the SPA process.

METHODS

Functionality of the Genetic Algorithm

We implement a genetic algorithm, which optimizes a problem by iteratively improving a population of candidate solutions through the operations of selection, crossover, and mutation. Each solution, referred to as an individual, is represented by a vector $x \in \mathbb{R}^L$, where $L$ denotes the number of parameters in the genome. Since the objective is to find a set of parameters that enable mode-selective single-photon addition/subtraction, each vector $x$ encodes a specific configuration of tunable waveguide and pump properties.

Initialization of the Population

A population $\mathcal{P} = \{x^{(i)}\}_{i=1}^N$ of $N$ individuals is initialized. Each individual

$$x^{(i)} = (x_{1}^{(i)}, x_{2}^{(i)}, \dots, x_{L}^{(i)}) \in \mathbb{R}^{L}$$

represents a candidate configuration, such as geometrical waveguide parameters and the spectral pump width. The components $x_j^{(i)}$ are drawn independently from normal distributions $N(\mu_j, \sigma_j)$, subject to a hard constraint $x_j^{(i)} \ge x_{j,\min}$, ensuring manufacturing feasibility:

$$x_j^{(i)} \sim N(\mu_j, \sigma_j), \quad \text{with} \quad x_j^{(i)} \ge x_{j,\min}.$$

Parent Selection

The goal of the optimization process is to iteratively move the population toward individuals with smaller fitness values, i.e., those closer to the target configuration. The fitness value, as described below, quantifies how close an individual is to the target configuration.

To this end, two parent individuals $x^{(p1)}, x^{(p2)} \in \mathcal{P}$ are selected based on their fitness values. The selection is carried out via the tournament selection method. In each selection round, a random subset $T \subset \mathcal{P}$ of the current population is drawn, where $|T| = N/2$. The individual within the subset that
exhibits the best fitness value $f(x)$ is selected as a parent:

$$x_{\text{best}} = \arg \min_{x \in T} f(x).$$

This process is repeated independently to obtain both parents. The selection mechanism thus introduces a degree of stochasticity while preserving selective pressure toward fitter individuals.

### Crossover

To generate new individuals, we apply a linear crossover between two selected parents. Given two parent vectors $x^{(p1)}$ and $x^{(p2)} \in \mathbb{R}^L$, two offspring $x^{(c1)}$ and $x^{(c2)}$ are created by forming weighted linear combinations of the parents. Specifically, each component $i \in \{1, \dots, L\}$ of the offspring is computed as

$$x_i^{(c1)} = \beta_1 \cdot x_i^{(p1)} + (1 - \beta_1) \cdot x_i^{(p2)},$$

$$x_i^{(c2)} = \beta_2 \cdot x_i^{(p1)} + (1 - \beta_2) \cdot x_i^{(p2)},$$

where the coefficients $\beta_1, \beta_2$ are drawn independently and uniformly from the interval $[-d, d+1]$, where $d$ is a hyperparameter, that we set to, for example, $d = 0.5$. This mechanism allows both interpolation ($\beta_i \in [0, 1]$) and extrapolation ($\beta_i < 0$ or $\beta_i > 1$) around the parent vectors, depending on the values of $\beta_k$.

This method promotes genetic diversity by enabling exploration of a broader parameter space. To ensure physical validity, the procedure is repeated until both offspring satisfy the constraint $x_i^{(c)} \ge \min_i$ for all components $i \in \{1, \dots, L\}$.

### Mutation

To maintain genetic diversity and avoid convergence to local minima, random mutations are applied to individual genomes. With a fixed probability defined by the mutation rate, a given individual is selected for mutation. In that case, each component of the genome is perturbed by adding Gaussian noise with zero mean and a standard deviation matching the corresponding initialization value. Specifically, for each gene $i \in \{1, \dots, L\}$, the mutated value is sampled as

$$x_i^{(\text{mut})} = x_i + \mathcal{N}(0, \sigma_i),$$

with the additional constraint $x_i^{(\text{mut})} \ge \min_i$. If any component fails to meet this constraint, the sampling is repeated until a valid configuration is obtained.

### Fitness Evaluation and Simulation

Each individual $x^{(i)} \in \mathcal{P}$ is evaluated by simulating the nonlinear process defined by its parameters. These parameters are passed to the simulation to compute the resulting joint spectral amplitude (JSA) or transfer function (TF). From there, a fitness value is derived, quantifying how well the configuration meets the desired physical behavior — such as process purity. The quality of a given configuration is quantified using the Schmidt number $K$, which serves as an indicator of mode selectivity. The ideal value corresponds to $K_s = 1$, representing perfect single-mode behavior. The fitness value is defined as the absolute deviation from this ideal:

$$f(x^{(i)}) = |K_s - K(x^{(i)})|.$$

To ensure that the solution performs robustly across multiple relevant spectral profiles, the simulation was also carried out using the first three Hermite-Gauss pump modes as input. The fitness score in this case is computed as the average deviation across all test cases:

$$f(x^{(i)}) = \frac{1}{3} \sum_{k=1}^{3} |K_s - K^{(k)}(x^{(i)})|,$$
where $K^{(k)}(x^{(i)})$ denotes the Schmidt number obtained when using the $k$-th Hermite-Gauss pump mode. It should be noted that both approaches—evaluating a single mode or averaging over multiple modes—led qualitatively to similar results, as the orientation of the phase-matching function was found to be the primary limiting factor.

## Iterative Population Update

The evolutionary process is carried out over $G$ generations. In each generation, a subset of individuals is evaluated based on their fitness, and parent individuals are selected via tournament selection. Crossover and mutation are then applied to generate new offspring. The resulting population fully replaces the previous one, forming the basis for the next generation. Over successive generations, the optimization gradually converges toward configurations exhibiting strong mode selectivity, as indicated by consistently low fitness values in the best-performing individuals.