Design of the coherent diffusive photon gun for generating non-classical states

TL;DR

This paper proposes a practical, integrated quantum light source using coherent diffusive photonics in waveguides, capable of generating non-classical states efficiently and feasibly with simplified models.

Contribution

It introduces a novel design for a deterministic non-classical light source based on dissipatively coupled waveguides, with a comprehensive analysis of model validity and feasibility.

Findings

Feasibility of generating bright sub-Poissonian light states from coherent input.

Generation of non-classical states occurs early in the dynamics, simplifying modeling.

Linear models can accurately describe the initial non-classical state generation.

Abstract

We suggest and discuss a concept of deterministic integrated source of non-classical light based on the coherent diffusive photonics, a coherent light flow in a system of dissipatively coupled waveguides. We show how this practical quantum device can be realized with a system of single-mode waveguides laser-inscribed in nonlinear glass. We describe a hierarchy of models, from the complete multi-mode model of the waveguide network to the single mode coupled to a bath, analyze the conditions for validity of the simplest single-mode model and demonstrate feasibility of the generation of bright sub-Poissonian light states merely from a coherent input. Notably, the generation of non-classical states occurs at the initial stages of the dynamics, and can be accounted for in the linear model that allows us to circumvent the prohibiting computational complexity of the exact full quantum representation.

Equations

$$\frac{d}{dt}\rho=-i[H,\rho]+\gamma_{1}\Bigl{(}\mathcal{L}(a)+\mathcal{L}(b)+\sum\limits_{j=0}^{N}\mathcal{L}(c_{j})\Bigr{)}\rho,$$

$$H^{\mathrm{int}}=g_a{a}^{†}{c}_0+g_b{b}^{†}{c}_0+\sum _{j=1}^N{g}_j{c}_{j-1}^{†}{c}_j+h.c.,$$

$$H^{\mathrm{kerr}}=\frac{U}{2}\sum _{\forall x}(x^{†}{)}^2{x}^2,x=a,b,c_j,j=0,1\dots N;$$

$$\frac{d}{dt}\rho_{3}=-i[H_{3},\rho_{3}]+\Bigl{(}\gamma_{1}\mathcal{L}(a)+\gamma_{1}\mathcal{L}(b)+\gamma\mathcal{L}(c_{0})\Bigr{)}\rho_{3},$$

$$s_{+}=\frac{1}{G}(g_{a}a+g_{b}b),s_{-}=\frac{1}{G}(g_{a}b-g_{b}a)$$

$$\frac{d}{dt}{\rho }_2=-i[H_2,{\rho }_2]+\left({\gamma }_1\mathcal{L}(s_{-})+(\Gamma +{\gamma }_1)\mathcal{L}(s_{+})\right){\rho }_2,$$

$$H_2^{\mathrm{self}}={\varsigma }_1(n_{+}^2+n_{-}^2)+{\varsigma }_2{n}_{+}{n}_{-}+{\varsigma }_3(n_{+}+n_{-}),$$

$$H_2^{\mathrm{int}}={\varsigma }_4(s_{+}^{†}{s}_{-}{)}^2+{\varsigma }_5{s}_{+}^{†}{s}_{-}(n_{-}-n_{+}-1)+h.c.$$

$${\varsigma }_1=\frac{U}{2G^4}(g_a^4+g_b^4),{\varsigma }_2=\frac{4U}{G^4}(g_a{g}_b{)}^2,{\varsigma }_3={\varsigma }_2/4-U/2,$$

$${\varsigma }_4={\varsigma }_2/4,{\varsigma }_5=\frac{U}{G^4}{g}_a{g}_b(g_a^2-g_b^2).$$

$$\frac{d}{dt}{\rho }_1=\left({\gamma }_1\mathcal{L}(s_{-})+{\gamma }_2\mathcal{L}(s_{-}^2)+{\gamma }_3\mathcal{L}(n_{-}{s}_{-})\right){\rho }_1+$$

$$-i[{\varsigma }_1{n}_{-}^2+{\varsigma }_3{n}_{-},{\rho }_1],$$

$${\gamma }_2=\frac{4U^2(g_a{g}_b{)}^4}{G^8(\Gamma +{\gamma }_1)},{\gamma }_3=\frac{4U^2(g_a{g}_b{)}^2}{G^8(\Gamma +{\gamma }_1)}(g_a^2-g_b^2{)}^2.$$

$$Q=\frac{⟨(s_{-}^{†}{)}^2{s}_{-}^2⟩}{⟨n_{-}⟩}-⟨n_{-}⟩.$$

$$\frac{d}{dt}{\rho }_n=-({\gamma }_{1}n+{\gamma }_{2}n(n-1)+{\gamma }_{3}n(n-1{)}^2){\rho }_n+$$

$$({\gamma }_1(n+1)+{\gamma }_3(n+1)n^2){\rho }_{n+1}+$$

$${\gamma }_2(n+1)(n+2){\rho }_{n+2},$$

$${\left(g_b/g_a\right)}_{\mathrm{opt}}=\sqrt{2}-1.$$

$${\gamma }_3=\frac{U^2}{4({\gamma }_1+\Gamma )}$$

$$X={\gamma }_3(⟨n_{-}(0)⟩{)}^2{t}_{fix},Y=\frac{{\gamma }_1}{{\gamma }_3⟨n_{-}(0){⟩}^2},$$

$${\rho }_n(t)=\frac{1}{\sqrt{2\pi }\sigma (t)}\exp \left\{-\frac{(n-⟨n_{-}(t)⟩{)}^2}{2\sigma (t{)}^2}\right\}$$

$$⟨n_{-}(t{)}^m⟩=\underset{0}{\overset{+\infty }{\int }}dn{\rho }_n{n}^m.$$

$$\frac{d}{dt}Q\approx -{\gamma }_{1}Q-2{\gamma }_2⟨n_{-}⟩(3Q+1)-{\gamma }_3⟨n_{-}{⟩}^2(5Q+4).$$

$$Q(X)\approx -\frac{4}{5}\left(\frac{1}{(1+2X{)}^{5/2}}-1\right).$$

$$U=2\hslash \omega \frac{\omega }{V_{\mathrm{eff}}}\frac{n_2}{n_{\mathrm{eff}}},$$

$$A_{\mathrm{eff}}^{-1}=\int \int d^2\vec{s}|E(\vec{s}){|}^4,$$

$$\int \int d^2\vec{s}|E(\vec{s}){|}^2=1.$$

$${\bar{\gamma }}^{\mathrm{NL}}=\frac{\omega }{c}\frac{n_2}{A_{\mathrm{eff}}}\Rightarrow U=2\frac{\hslash \omega }{T_{\mathrm{eff}}{n}_{\mathrm{eff}}}{\bar{\gamma }}^{\mathrm{NL}},$$

$${\partial }_t⟨s_{-}⟩=c_1⟨s_{-}⟩+c_2⟨s_{-}^{†}{s}_{-}^2⟩-\frac{{\gamma }_3}{2}⟨s_{-}^{†2}{s}_{-}^3⟩$$

$${\partial }_t⟨s_{-}^{†}⟩=c_1^{\ast }⟨s_{-}^{†}⟩+c_2^{\ast }⟨s_{-}^{†2}{s}_{-}⟩-\frac{{\gamma }_3}{2}⟨s_{-}^{†3}{s}_{-}^2⟩$$

$${\partial }_t⟨s_{-}^2⟩=c_3⟨s_{-}^2⟩+c_4⟨s_{-}^{†}{s}_{-}^3⟩-{\gamma }_3⟨s_{-}^{†2}{s}_{-}^4⟩$$

$${\partial }_t⟨s_{-}^{†2}⟩=c_3^{\ast }⟨s_{-}^{†2}⟩+c_4^{\ast }⟨s_{-}^{†3}{s}_{-}⟩-{\gamma }_3⟨s_{-}^{†4}{s}_{-}^2⟩$$

$${\partial }_t⟨s_{-}^{†}{s}_{-}⟩=-{\gamma }_1⟨s_{-}^{†}{s}_{-}⟩+c_5⟨s_{-}^{†2}{s}_{-}^2⟩-{\gamma }_3⟨s_{-}^{†3}{s}_{-}^3⟩$$

$$H_2^{int}={\varsigma }_4(s_{+}^{†}{s}_{-}{)}^2+h.c.$$

$$|s_{-}{⟩}^{\mathrm{out}}|s_{+}{⟩}^{\mathrm{out}}=|2{⟩}_{-}|0{⟩}_{+}=$$

Full text

I Introduction

Engineered loss has already turned into a powerful and intensively researched tool for quantum state manipulation. With the help of engineered reservoirs, it is possible to drive a target system into a desired state, to generate and protect entanglement, to implement computation, and to transfer quantum states carmichael ; cirac ; metel . Engineered reservoirs have become a popular and fruitful research direction for QED clark , cold atoms and trapped ions zoller ; blatt , Rydberg atoms carr ; rog , and for super-conducting circuits science2015 ; kimchi ; kapit due to the large spectra of possibilities to design systems with necessary coupling between components and strong effective nonlinearities. For example, by devising strong self-Kerr and cross-Kerr interactions, and combining them with two-photon conversion and classical driving, it is possible to produce two-photon loss, and in this way create the Schrödinger-cat states and safeguard them against linear loss science2015 .

It is considerably harder to engineer losses in photonic circuits. Optical nonlinearities are usually quite small and accompanied by rather strong conventional linear loss. However, even in this case there are possibilities to exploit effects of nonlocal loss and devise photonic circuits with quite unusual and potentially very useful features, even in passive optical structures. These possibilities are opened up by a new level of precision and controllability reached in systems of coupled waveguides. For example, localization of the light in the small part of a perfect Lieb lattice and diffraction-less propagation was experimentally demonstrated in a system of single-mode waveguides written in a bulk glass seba1 ; vicencio .

These new experimental advances gave rise to recently introduced concept of “coherent diffusive photonics” (CDP): light-processing in a specifically designed system of single-mode waveguides coupled by common loss reservoirs, as realized in practice in natcom . Even with linear glass, dissipative coupling enables, for example, light equalization in a waveguide chain: any input state tends to a completely symmetrized state over all the modes. In contrast to the common unitary coupling schemes, CDP schemes are quite robust with respect to variations in the coupling length and strength. Further, CDP allows one to realize an optical router directing light in different arms of the structure by selective excitation of the central control nodes, and CDP can even allow light to be localized in a perfect lattice of dissipatively coupled modes natcom .

In this paper we extend the developed CDP schemes to the nonlinear regime, and suggest a new integrated CDP-based device exhibiting features of a “photon gun” PhoG: a deterministic generator of nonclassical states, in particular, sub-Poissonian states (Fig. 1).

It is worth noticing that a true photon gun is still a challenge (for review see SPSreview1 ; SPSreview2 ; SPSreview3 ). Due to high technical overhead, the available single photon sources are not readily implementable, mostly do not have well-defined output mode and many of them cannot simultaneously satisfy the crucial figures of merit: purity, indistinguishability (often traded for brightness), and efficiency. Only very recently top-end quantum dot-based sources could demonstrate simultaneously single-photon purity of $\geq 99\%$, photon indistinguishability of $95-99\%$ and extraction efficiency of up to $65\%$ qdot . Other alternative, heralded photon sources based on entangled photon pairs production are probabilistic and deliver high indistinguishability at the cost of low brightness. Therefore in many quantum technology applications, experimentalists revert to attenuated coherent states instead of single photons. This is true for most if not all quantum cryptography implementations and many other, including quantum simulations and implementations of quantum gates. Dim coherent states are ‘cheap’, easy to handle and deterministic, however their Poissonian photon number statistics represents a serious limitation, resulting in unwanted contribution from multi-photon components. Much better performance can be achieved using light with sub-Poissonian photon statistics as a quasi-single photon source. For example, the key rates for decoy-state quantum key distribution (QKD) have been explicitely calculated for different sources, and it was demonstrated that such quasi-single-photon source can drastically raise the key rate in the decoy-state QKD wang ; wang1 .

We demonstrate that the photon gun device suggested in this paper is potentially able to act either as a deterministic or, upon a slight variation of the scheme, probabilistic quasi-single photon source. We analyze the possibility of realization of these devices with current technology. Whereas approaching deterministic generation of quasi-single photon- or few-photon strongly sub-Poissoinian states still seems challenging, we show that bright sub-Poissonian photon gun is completely feasible and can be quite robust with respect to experimental imperfections and noise.

Bright Sub-Poissonian photon gun (PhoG) suggested in this paper carries important advantage of producing non-classical states from merely coherent input, on demand, and in well-defined temporal and spatial modes. Further, PhoG is a versatile quantum source and can be modified to produce correlated photons and other entangled states at the output. Knowingly, non-classical photon number correlations and sub-Poissonian photon statistics can be exploited for quantum sensing, imaging and other metrology tasks bio ; imaging .

The basic idea underlying our scheme is to create a photonic circuit exhibiting nonlinear coherent loss (NCL), which is a specific kind of engineered single-photon loss mog2010 ; mog2011 . The principal circuit, Fig. 1(a), comprises two signal modes asymmetrically coupled to a common loss reservoir and undergoing NCL. The required light propagation regime can be engineered in different ways. We suggest a nonlinear waveguide network with a chain waveguide structure working as a common reservoir, Fig. 1(b), which allows the NCL to be realized in the superposition of two waveguide modes. We analyze its feasibility under realistic Kerr nonlinearity of nonlinear glass and for realistic waveguide parameters. We show that it is feasible to deterministically create bright sub-Poissonian states of light using a compact (cm-sized) integrated CDP device and, furthermore, reach other useful nonclassical output states. The generation of nonclassical light occurs in the regime when a conventional single-photon loss is still not affecting significantly the process of the state generation.

The outline of the paper is as follows. Firstly, in Sec. II, we present the CDP scheme and describe the hierarchy of models which can be used to account for the state dynamics. We reveal how the description of the whole network can be reduced to the description of just a single superposition mode (“single-mode model”). In Sec. III, we show that the regime of interest for us is the initial stage of the dynamics, when strong non-classicality can be obtained. We develop an analytic model for evolution of the Mandel parameter and also show that at this stage the dynamics can be well captured by the linearization approach for the quantum correction to the field amplitude. The single-mode model allows to present a simple and illustrative analysis of the PhoG device feasibility providing necessary parameters of the set-up, and ways to optimize it.

In Secs. IV, V we analyze dynamics of the system beyond the simplest single collective mode model. We pay special attention to the two-mode model (Sec. IV) which exhibits a number of nonclassical features not captured by the single-mode model, such as entanglement generation. After that in Sec. VI we proceed to the complete network showing that for the initial stage of dynamics the predictions made within the single-mode model are indeed true. Strong photon-number squeezing is shown to emerge for moderate lengths of the CDP circuit with realistic losses and nonlinearity. Also, to assess the impact of specific effects of short-pulse propagation in our system of coupled waveguides, we dynamically evolved the spectral and temporal properties of the pulse via the Nonlinear Schrodinger Equation, and analysed their influence on the effectiveness of the NCL mechanism taking into account the combined effect of chromatic dispersion, self-phase modulation, and self-steepening of the pulse.

II Hierarchy of models

In the case of a large number of interacting quantum modes, a large number of photons and in the presence of nonlinearity, an exact description of the quantum state dynamics of our system is a formidable task. To tackle the problem, we have developed the hierarchical approach to model the photon gun and optimize its structure. The complete CDP scheme is reduced to much simpler systems of a few modes as shown in Fig. 2. These systems are much more tractable and can be analyzed, optimization can be carried out, approximations can be verified and then the obtained results are validated by the approximate numerical analysis of the complete CDP circuit.

The original scheme underlying all the subsequent discussion is depicted in Figs. 1(b,c): two waveguides $a$ and $b$ are coupled to the third waveguide $c_{0}$ (but not to each other), and this waveguide $c_{0}$ is coupled to the “tail” of next-neighbour coupled single-mode waveguides $c_{j}$, $j=1\ldots N$. In our consideration, we assume identical single-mode waveguides described by the corresponding bosonic creation and annihilation operators. We assume that some initial state is created in the network, and then this states evolves. Thus, the dynamics of this system is described by the following master equation for the density matrix $\rho$:

[TABLE]

where the Lindbladians $\mathcal{L}$ describing common single-photon loss with the rate $\gamma_{1}$ are given by $\mathcal{L}(x)y=xyx^{\dagger}-\frac{1}{2}x^{\dagger}xy-\frac{1}{2}yx^{\dagger}x$ for operators $x$ and $y$; the density matrix for all modes is $\rho$ and the Hamiltonian $H=H^{\rm int}+H^{\rm kerr}$ consists of two parts. The first part describes linear interactions between the modes

[TABLE]

where $g_{a,b}$ are couplings of the modes $a,b$, respectively, to the third mode $c_{0}$, and $g_{j}$ is the coupling between mode $c_{j-1}$ and $c_{j}$. The part $H^{\rm kerr}$ describes the self-Kerr interaction of each mode,

[TABLE]

$U$ is the Kerr nonlinear interaction constant.

For large number of waveguides $N$ and large initial number of photons, the complete model Eqs. (1,2) is practically intractable without approximations. However, we are interested in the regime when the “tail” formed by the waveguides $c_{j}$, $j=1\ldots N$, can be considered a dissipative reservoir rapidly guiding away field from the mode $c_{0}$ natcom ; vicencio . Assuming that the tail functions as a Markovian reservoir with the decay rate $\gamma_{c}$, we arrive at the three-mode model depicted in Fig. 2(b) and described by the following master equation for the three-mode density matrix:

[TABLE]

where $\gamma=\gamma_{1}+\gamma_{c}$, the Hamiltonian is $H_{3}=H^{\rm int}_{3}+H^{\rm kerr}_{3}$ with $H^{\rm kerr}_{3}$ as given in Eq. (3) and the interaction part is $H^{\rm int}_{3}=g_{a}a^{\dagger}c_{0}+g_{b}b^{\dagger}c_{0}+h.c$.

The three mode model of Fig. 2(b) is still too complicated for exact analysis in case of large number of photons. Further simplification can be reached assuming that the decay rate into the“tail” $\gamma_{c}$ is large enough to allow for the adiabatic elimination of the mode $c_{0}$. Let us introduce symmetric and anti-symmetric collective modes:

[TABLE]

with $G=\sqrt{g_{a}^{2}+g_{b}^{2}}$ and coupling constants $g_{a,b}$ assumed real. Then adiabatic elimination of mode $c_{0}$ leads to the following two-mode master equation mog2010 ; mog2011 :

[TABLE]

with $\Gamma=4G^{2}/\gamma$. The Hamiltonian $H_{2}=H^{\rm int}_{2}+H^{\rm self}_{2}$ is given by

[TABLE]

where $n_{\pm}=s_{\pm}^{\dagger}s_{\pm}$ are the photon number operators, and the coefficients read

[TABLE]

The two-mode model is depicted in Fig. 2(c).

Finally, if the state of the superposition mode $s_{+}$ decays to the vacuum much quicker than the typical time-scale of the dynamics of the superposition mode $s_{-}$, one arrives to the single-mode model of Fig. 2(d) described by the following master equation mog2010 ; mog2011 :

[TABLE]

where the decay rates are given by

[TABLE]

The master equation (8) takes account of the three decay channels for the anti-symmetric collective mode: conventional single photon loss ($\gamma_{1}$), two-photon loss ($\gamma_{2}$) and nonlinear coherent loss NCL ($\gamma_{3}$). We refer to the decay described by the Lindblad operator $n_{-}s_{-}$ as “nonlinear coherent loss” as the eigenstates of the operator $n_{-}s_{-}$ were named “nonlinear coherent states” vogel . The Lindblad operator $n_{-}s_{-}$ responsible for NCL can be considered as the annihilation operator of so-called $f$-deformed harmonic oscillator manko .

Note that the unitary part of the master equation, the last term in (8) connected to the self-Kerr interaction, does not manifest itself in the dynamics of the diagonal elements of the density matrix $\rho_{1}$.

The model (8) is much easier to tackle for large number of photons than the two-mode model (6). Moreover, it allows us to obtain some general conclusions about non-classicality which are valid even for numbers of photons so large that an exact solution becomes unfeasible even for the single-mode model.

III Single-mode model

In this section we analyze the simplest single-mode scheme of Fig. 2(d) to describe the generation of the bright sub-Poissonian state from the semi-classical coherent input, and outline the way to optimize the scheme by choosing the amplitude of the initial coherent state and parameters of the PhoG for realistic CDP circuits with significant linear loss and moderate Kerr nonlinearity. We establish an existence of universal parameters $X,Y$, which are determined by effective nonlinearity and thus allow us to characterize the scheme performance even for realistic small optical nonlinearities and, consequently, very large numbers of photons, when solving the master equation becomes hardly possible. Finally, we demonstrate practical feasibility of the scheme.

III.1 Nonlinear loss dynamics

There are two kinds of nonlinear loss present in the single-mode master equation (8). The first one is the two-photon loss described by the term $\mathcal{L}(s_{-}^{2})$ and the second kind is the NCL described by the term $\mathcal{L}(n_{-}s_{-})$. Both these kinds of losses are able to produce photon-number squeezing. However, their time-scales and asymptotic states are different. The two-photon loss asymptotically leads to a mixture of single-photon and the vacuum states ezaki ; alex . NCL asymptotically leads to the single-photon state mog2010 . We show here that for the initial coherent state with large amplitude, $|\alpha|\gg 1$ and the assumed equal decay rates, $\gamma_{2}=\gamma_{3}$, NCL leads to much faster decay of the average photon number in the anti-symmetric collective mode $s_{-}$ than the two-photon loss and thus to more rapid narrowing of the photon number distribution in $s_{-}$.

The degree of photon-number distribution squeezing can be conveniently described by the Mandel parameter Davidovich1996 , which for the superposition mode $s_{-}$ reads:

[TABLE]

Mandel parameter $Q={-}1$ corresponds to the perfect squeezing (Fock states), zero is for coherent states. This parameter can be inferred from dynamics of only the diagonal elements of the density matrix (8) given by the following equation

[TABLE]

where $\langle n|\rho_{1}(t)|n\rangle=\rho_{n}$ are the density matrix elements in the Fock-state basis.

Note that for the symmetric case, $g_{a}=g_{b}$, the NCL decay rate $\gamma_{3}$ is zero. For asymmetric coupling, $g_{a}\neq g_{b}$, we fix the parameter $g_{a}$ and start varying the parameter $x=g_{b}/g_{a}$ for $x\in[0,1]$. The maximal value of the rate $\gamma_{3}$ is then given by the condition $d\gamma_{3}(x)/dx=0$ and leads to the following result for the ratio of the interaction constants to maximize NCL:

[TABLE]

For this optimal coupling ratio, we have

[TABLE]

In further considerations, we adopt this ratio for estimations.

In Fig. 3 examples of the Mandel parameters and average photon number dynamics are given for different numbers of photons (several hundred) of the initial coherent state. Fig. 3 reveals a number of important features in the dynamics of the antisymmetric collective mode:

The Mandel parameter rapidly decreases at the initial stage of the dynamics. Only a small percentage of photons is lost when quite considerable negative values of the Mandel parameter are reached (less than $15\%$ to have significant photon-number squeezing with $Q<-0.25$, see Figs. 3(c,d)). 2. 2.

The practically reachable values of the Mandel parameter are approximately limited to $Q=-0.8$ (see also Ref. mog2013 ). This minimal value is shifted to smaller interaction length by increasing the average number of photons in the initial state. 3. 3.

Dynamics due to nonlinear coherent loss (NCL) occurs much faster than dynamics due to the two-photon loss. As seen in Fig. 3(a), even hundredfold increase in the two-photon decay rate, $\gamma_{2}$, does not influence much the initial stage of the $Q$-parameter dynamics. In general, the ratio between these two time scales is determined by the ratio between the coupling constants, $g_{b}/g_{a}$: $\gamma_{3}=\gamma_{2}\left(1-\left(g_{b}/g_{a}\right)^{2}\right)/\left(g_{b}/g_{a}\right)^{2}$ (see Eq. (9) and text around Eq. (12)). 4. 4.

NCL can outrun linear loss. Fig. 3(b) shows that even large linear loss weakly affects the initial stage of the dynamics. 5. 5.

Consider now the case of no linear loss. For different large initial numbers of photons, after some time the average numbers of photons tends to largely the same stationary value (Fig. 3(c)). It is a signature behaviour of NCL mog2010 .

III.2 Universal parameters

As we have shown, the two-photon loss only weakly influences the dynamics of photon number squeezing. Also, large photon-number squeezing is generated at the initial stage of dynamics when only relatively small number of photons is lost. These two facts allow us to introduce two dimensionless parameters allowing an estimation of the feasibility of the scheme. These two parameters are

[TABLE]

where $t_{fix}$ is the fixed interaction time (defined by the length of the PhoG). The parameter $X$ defines the reachable value of the Mandel parameter, $Q$, in absence of the linear and two-photon loss, and the parameter $Y$ defines the tolerable level of the linear loss when the value of $Q$ is still defined by $X$.

In terms of parameter $Y$, one has to have $Y<1$ for optimal $Q$ for a given $X$, Fig. 4 (a). Of course, even larger values of $Y$ still allow for considerable $Q$. However for that, larger values of $X$ are required.

Such a dependence of photon number squeezing on the universal parameters can be explained by simple analytic considerations. For large number of photons of the initial coherent state, the initial photon number distribution can be with high precision described by the Gaussian. Let us assume that for initial stage of the dynamics the distribution can still be described by the Gaussian

[TABLE]

where $\sigma(t)=(Q(t)+1)\langle n_{-}(t)\rangle$. In the Appendix A it is shown that Gaussianity of the photon-number distribution holds very well for $\langle n_{-}(t)\rangle\gg 1$. Assumption of a smooth photon number distribution slowly changing with $n$ (like the one given by Eq.(15) for large $\langle n_{-}(0)\rangle$) allows to introduce a continuous variable $n$, to make an approximation $\rho_{n+m}\approx\rho_{n}+m\frac{d}{dn}\rho_{n}+m^{2}\frac{1}{2}\frac{d^{2}}{dn^{2}}\rho_{n}$, and to calculate averages as

[TABLE]

Thus, one can obtain the following equation for the Mandel parameter

[TABLE]

Neglecting linear and two-photon losses, in the continuous approximation for Eq. (11) one obtains that $\langle n_{-}(X)\rangle\approx\langle n_{-}(0)\rangle/\sqrt{1+2X}$ (see the Appendix A). For the Mandel parameter we get:

[TABLE]

This situation is illustrated in Fig.4. The panel 4(a) shows that in the absence of linear loss neither the initial number of photons nor the interaction time by themselves are the parameters defining possible $Q$. Practically, it depends only on $X$. The black solid line in the panel 4(a) shows the result provided by the approximation (17) for several initial states with hundreds of photons. It is obvious that Eq.(17) very closely reproduce the exact solution of Eq. (11). Eq. (11) also shows why the influence of two-photon losses is so small for $\gamma_{2}$ close to $\gamma_{3}$: the time scales are determined by $\propto\gamma_{2}n^{2}$ for two-photon loss and $\propto\gamma_{3}n^{3}$ for NCL. For large photon numbers and $\gamma_{2}\approx\gamma_{3}$, the NCL obviously dominates.

The panel 4(b) shows an influence of linear loss (described by the parameter $Y$) on the reachable value of the Mandel parameter. For a given linear loss of $\gamma_{1}=10g$, for $\langle n(0)\rangle=900$ one has $Y<1$, and the curve describing the achievable values of the Mandel parameter (dash-dotted line) is practically coinciding with the approximation given by Eq. (17). For small $Y$ any influence of the linear loss is indeed small. For smaller initial number of photons, $\langle n(0)\rangle=100$, when $Y>1$, one has considerable difference with analytic predictions (dashed line in Fig.4(b)).

III.3 Feasibility

Universal parameters $X$ and $Y$ are decisive to ensure the optimal performance of PhoG as a deterministic source of sub-Poissonian light. For PhoG design we have a “rule-of-thumb”: we should aim to make $X$ as large as possible, and $Y$ as small as possible. In this section we demonstrate that it is feasible to reach $X>0.3$ and $Y<1$ for realistic waveguide structures.

First of all, we connect the nonlinear interaction constant $U$ in the expressions (9) with quantities commonly used for description of the Kerr-nonlinearity in waveguide structures. We assume that the pulse of the finite length propagates through the waveguide, and the pulse length $L_{\rm eff}$ is much smaller than the waveguide length; we introduce the propagation length as $z=ct_{\rm fix}/n_{\rm eff}$, where $n_{\rm eff}$ is the effective refractive index of the waveguide mode. Thus, the nonlinear interaction constant in our model can be written in the following way (see, for example, Refs. drummond ):

[TABLE]

where $n_{2}$ is the nonlinear refractive index and $V_{\rm eff}$ is the mode volume. For our model, we assume $V_{\rm eff}\approx A_{\rm eff}L_{\rm eff}$, where $A_{\rm eff}$ is the effective area. The effective area is defined as

[TABLE]

where $E(\vec{s})$ is the normalized transversal field profile,

[TABLE]

Here we do not consider influence of losses on the nonlinear interaction constant coso .

In a fiber-like waveguides it is customary to introduce the nonlinear fiber coefficient:

[TABLE]

where $T_{\rm eff}$ is the pulse time-duration in vacuum.

Now let us consider several waveguide arrangement with realistic parameters.

III.3.1 Large waveguides in bulk glass

Let us take coupling parameters as used in the CDP circuits in Ref. natcom , $g_{a}\approx 200m^{-1}$. The decay rate of the $c_{0}$ mode Fig. 2 can be approximated by $\gamma_{c}\approx 4\sqrt{g_{a}^{2}+g_{b}^{2}}$, which for optimal coupling ratio results in $\Gamma\approx 216.5$m*-1* (decay rate for the symmetric collective mode $s_{+}$, Fig. 2). For intended linear loss of about $0.5$ dB/cm, we have $\gamma_{1}\approx 11.5$m*-1*, and $\gamma_{3}\approx 0.0011$m${}^{-1}\times U^{2}$. The dimensions of the waveguides as used in CDP circuits in Ref. natcom are approximately $4\mu$m$\times 4\mu$m. For a conservative estimate, let us assume an effective modal area an order of magnitude larger than typical for single-mode fibers with the similar core dimensions around $1000$ nm (see, for example, Corning® HI 1060 fiber with the mode-field diameter of about 6.2 $\mu$m at $1060$  nm), $A_{\rm eff}=300\mu m^{2}$ (which would be more than three times larger even for typical modal areas at 1550 nm). Thus, taking $n_{2}=3\times 10^{-18}$W*-1m2 and $n_{\rm eff}=2.59$ typical for IG2 glass ig2 , for a $100$ fs pulse at $1060$ nm we have $\gamma_{3}\approx 8\times 10^{-18}$m-1*. Condition $Y=1$ of Eq. (14) requires energy levels of $1.2\times 10^{9}$ photons per pulse or $224$pJ per pulse. The feasibility of these energy levels in the context of our set-up is confirmed by the recent work robert2018 , where femtosecond pulses with energies more than $10$ nJ at $1030$ nm has been used for writing waveguides. Using all these parameters, for a $3$ cm waveguide we have $X\approx 0.33$, which is sufficient to reach high degree of photon-number squeezing $Q<-0.5$. As follows from the condition (14), for lower effective modal area one would need proportionally lower energies per pulse. For $A_{\rm eff}=30\mu$m2 only $25$pJ per pulse is sufficient.

III.3.2 Fiber-like systems

Now let us use for the PhoG set-up the parameters from the recent fiber-based interferometric scheme Ref. aruto for sub-Poissonian light generation. For the used polarization maintaining fiber (Nufern, HP-780) one has ${\bar{\gamma}}^{\rm NL}=8.51\times 10^{-3}$W*-1m-1* for $808$ nm. Using $n_{\rm eff}=1.45$, we arrive at $\gamma_{3}\approx 9.2\times 10^{-19}$m*-1*. It is nearly an order of magnitude lower than for the bulk waveguides. However, much lower loss of $3.5$ dB/km implies that the condition (14) is easier to satisfy. For that, $3\times 10^{7}$ photons per pulse (or less than $10$pJ per pulse) are needed. If $2\times 10^{9}$ photons per pulse are used, $10$cm of the fiber are required to reach $X=0.37$ and squeezing of about $Q<-0.5$ (which is more or less in agreement with results obtained in Ref. aruto ).

III.3.3 Silicon nanowires

Finally, let us consider silicon nanowires arrangement as a possible platform for PhoG generators. Si-nanowire waveguides have typically very high nonlinearity and loss, and small lengths. We assume loss of $5$ dB/cm and moderate nonlinearity, ${\bar{\gamma}}^{\rm NL}=300$W*-1m-1* with other parameters being the same as above sato . Notice that even for such loss level, $\Gamma$ is larger than $\gamma_{1}$, where $\gamma_{1}\approx 115$ m*-1*. For $n_{\rm eff}=3.5$ and wavelength of $1064$ nm, we obtain $\gamma_{3}\approx 7.75\times 10^{-11}$m*-1*, which is much larger than for all the previously considered cases. To satisfy $Y=1$, we require only $1.3\times 10^{6}$ photons per pulse or less than $0.3$pJ pulse energy. For this number of photons a $3$ mm waveguide would be sufficient to achieve $X=0.345$ and $Q<-0.5$.

III.4 Linearization

Here we demonstrate that, over the timescales when large photon number squeezing occurs, at the initial stage of dynamics the evolution of the Mandel parameter can be adequately described using a linearization approach. This approximation is typical for consideration of nonlinear waveguide systems. Essentially, it is a linearization on a quantum correction to the creation/annihilation operators, $s_{-}=S_{-}+\delta s_{-}$, where the quantity $S_{-}=\langle s_{-}\rangle$ is the classical amplitude, $\delta s_{-}$ are the quantum fluctuations which obey the same commutator as $s_{-}$, and $\langle\delta s_{-}\rangle=0$. The resulting system of equations allows us to find quantities $\langle s_{-}\rangle$, $\langle s_{-}^{\dagger}\rangle$, $\langle s_{-}^{2}\rangle$, $\langle s_{-}^{\dagger}s_{-}\rangle$ and $\langle s_{-}^{\dagger 2}\rangle$ and circumvents the computational complexity required in the case of exact solution for large photon numbers in a large number of modes.

Our starting point is Eq. (8), from which we may derive the following system of coupled equations:

[TABLE]

with coefficients $c_{1}=\left(-\gamma_{1}/2+i\varsigma_{1}+i\varsigma_{3}\right)$, $c_{2}=\left(\gamma_{2}-\gamma_{3}+2i\varsigma_{1}\right)$, $c_{3}=\left(-\gamma_{1}-\gamma_{2}-\gamma_{3}+4i\varsigma_{1}+2i\varsigma_{3}\right)$, $c_{4}=\left(-2\gamma_{2}-5\gamma_{3}+4i\varsigma_{1}\right)$ and $c_{5}=\left(-2\gamma_{2}-\gamma_{3}\right)$.

This system of equations is not closed, and so cannot yet be solved. To proceed, we must perform the linearization procedure in order to reduce operator products in Eq. (III.4) to at most second-order, thereby closing the system of equations and allowing us to obtain its solution. The linearization can be done in a standard manner using the cumulant expansion (see Appendix B). The linearized forms of Eqs. (III.4) are excplitly shown in Appendix C.

The solution of the linearized forms of Eq. (III.4) are used to approximate $Q$, and in Fig. (5) this result is compared to the exact method from the previous section. The linearization approximation (dashed lines) accurately predicts the evolution of $Q$ over the initial stages of evolution, and actually underestimates $|Q|$ in the later stages, although including realistic $\gamma_{1}$ allows the approximation to remain accurate as the nonclassical output state is pushed towards the vacuum.

Since when $\gamma_{1}=0$ the linearization approximation remains accurate over the timescales of interest (see the panel Fig. 5(a)), and since the presence of realistic loss makes the approximation increasingly accurate (see the panel Fig. 5(b)), we may confidently apply the linearization approach over the parameter regimes of physical interest – short times and realistic loss – even in the case where a fully quantum treatment would be intractable, i.e. large $\langle n\left(0\right)\rangle$ or large number of modes.

IV The symmetric two-mode model

In the previous section, we have shown that two-photon loss is far less efficient in producing photon-number squeezing than NCL. However, it is worth pointing out that even symmetric system exhibiting only two-photon loss can be of use. To show it, let us recourse to the two-mode model of Fig. 2(c). First of all, even a symmetric system is able to produce entangled states of waveguide modes asymptotically, when the state of the symmetric superposition mode $s_{+}$ is the vacuum (Sec. IV.1). The second and rather non-trivial feature of the symmetric system is the possibility of producing few-photon entangled states by driving the rapidly decaying collective mode $s_{+}$ (Sec. IV.2). This scheme can be also used for probabilistic generating of single photons.

IV.1 Modal entanglement

Generally, the state produced by the two-photon loss is quite far from Gaussian (asymptotically, the initial bright coherent state is driven by the two-photon loss toward the superposition of the vacuum and single-photon state ez1 ; ez2 ). However, one can still obtain a lower bound on the entanglement between modes from the covariance matrix for modes $a,b$: $\sigma_{kl}=\frac{1}{2}\langle d_{k}d_{l}+d_{l}d_{k}\rangle-\langle d_{k}\rangle\langle d_{l}\rangle$, with the vector $\vec{d}=\left(1/\sqrt{2}\right)\left[a+a^{\dagger},i\left(a^{\dagger}-a\right),b+b^{\dagger},i\left(b^{\dagger}-b\right)\right]$ Adesso2007 ; Weedbrook2012 . This bound is provided by the (Gaussian) logarithmic negativity $\mathcal{N}$ Weedbrook2012 ; Adesso2007 , which quantifies violation of the PPT criterion Simon2000 . The logarithmic negativity $\mathcal{N}$ is defined as $\mathcal{N}=\max{\left\{0,-\log\lambda\right\}}$, where $\lambda$ is the smallest symplectic eigenvalue of the partially transposed matrix $\sigma_{kl}$. Elements of this matrix are provided by the linearized approach considered in the previous section.

For example, for $\langle a^{\dagger}a\left(0\right)\rangle=\langle b^{\dagger}b\left(0\right)\rangle=2500$, $U=2g$, $\gamma_{c}=15g$, $\gamma_{1}=11.5g$ and symmetric coupling $g_{a}=g_{b}=60g$, the system evolves to logarithmic negativities $\mathcal{N}\approx 1.25$ between modes $a$ and $b$ within $0.01$ $gt$, while modes $a$ and $b$ each contain $18$ photons. Thus symmetric PhoG generates entanglement between modes $a$ and $b$ (though is unable to produce large photon-number squeezing in the initial stages of the dynamics).

IV.2 Photon pairs generation

Despite the quick decay of the symmetric mode $s_{+}$, it still can create a non-classical state in the mode $s_{-}$. Indeed, for the symmetric case of Fig. 6, $g_{a}=g_{b}\equiv g$, we have $\varsigma_{5}=0$ and the Hamiltonian part describing interaction between superposition modes is

[TABLE]

The Hamiltonian $H^{int}_{2}$ for the case describes the four-wave mixing process of transferring two photons of initially excited decaying mode $s_{+}$ to the initially empty mode $s_{-}$. Thus the symmetric PhoG of Fig. 6 can be designed to reproduce the four-wave mixing process which is known to allow for the generation of entangled photon pairs kumar ; tak . But notice that in our case there is no problem of separating the generated state from the pump. For sufficiently long PhoG the driving excitation of the mode $s_{+}$ completely decays leaving photons only in the mode $s_{-}$.

The process of the pair generation dynamics with weak coherent initial states is illustrated in Fig. 7 with the photon-number distribution of the mode $a$. Initial coherent states with equal amplitudes are excited in the waveguides $a$ and $b$ (corresponding to $s_{+},s_{-}$ as shown in Fig. 6). For small propagation lengths, the remnant of the initial excitation dominates (Fig. 7(a)). However, it quickly decays and for larger times predominantly two photons in the anti-symmetric mode are observed (Fig. 7(b)). Just two photons in this mode corresponds to the following entangled superposition of photons in the waveguides $a$ and $b$ at the output of the PhoG:

[TABLE]

One can see the photon-number distribution corresponding to single- and two-photon part of the generated state (23) in panel Fig. 7(b). The ratio of the single-photon component to the two-photon component is approximately 2.05 for $gt=0.75$ . Notice that by increasing amplitudes of the initial states, a squeezed vacuum state is created in the mode $s_{-}$.

Another important difference with the conventional four-wave mixing pair-generating schemes kumar ; tak is that here we have the two entangled photons in two spatially separated modes. As can be seen in Ref. natcom , at the PhoG output, the spatial distance between the modes exceeds considerably the mode cross-section.

V Three-mode model

As has been seen in the previous Sections, analysis of the single-mode model allows us to give simple and straightforward prescriptions for the practical realizations of the PhoG. The question is how close predictions given by the single-mode model approach the results for more complicated modal arrangements.

In this Section we analyze dynamics of a three-mode PhoG schematically shown in Fig. 2(c). Considering this model, we have in mind two purposes. Firstly, the three mode scheme can be realized in practice as it is, i.e., by arranging a strong loss in the waveguide $c_{0}$. It is feasible, for example, in semiconductor planar waveguide structures balk . Secondly, the three-mode model still allows for a comparison between the exact solution and the single-mode approximation for a moderately large numbers of photons (few tens). One could also use it to verify the linearization approximation, and check the optimal ratio of the interaction constants $g_{b}/g_{a}$ for minimizing the Mandel parameter with a fixed interaction time/length.

To verify the transition from the three modes model (4) to the single mode model (8), the evolution of the mean photon number $\langle n_{{-}}(t)\rangle$ and the Mandel parameter $Q_{{-}}(t)$ of the antisymmetric collective mode $s_{{-}}$ was simulated using three modes model and quantum trajectories approach traj . Obtained results were then compared with the same quantities computed using the single mode model with parameters taken from equations (9). Figure 8 illustrates how close these models are for the optimal ratio $g_{b}/g_{a}$ (12) and large enough losses in the third, lossy, mode. Solid lines in the figure correspond to the three modes model and dots correspond to the single mode model. The figure was obtained for $g_{b}/g_{a}=\sqrt{2}-1$, $U=0.012g_{b}$, $\gamma_{c}=6.04g_{b}$, $\gamma_{1}=0$.

Thus, the three modes model confirms the main evolution characteristic stemming from the single mode analysis, namely, rapid photon number squeezing for relatively small decrease in photon number in the initial stage of evolution.

VI Multi-mode model

Finally, we consider a complete multi-mode model depicted in Fig. 2(a) for the realistic parameters described in Sec. III, and show that predictions made in that section still hold true for a multi-mode circuit of dissipatively coupled bosonic modes. We consider the multi-mode PhoG device with the two signal modes coupled to a long tail, which is described by Eq. (1).

We take the physical parameters from Sec. III.3.1 for coupled waveguides in bulk glass, and the number of tail modes to be $N=28$. Since the Kerr nonlinearity parameter is small, $U\sim 10^{-8}$, in order to reach the regime of strongly sub-Poissonian light we have to take a bright input coherent state containing $\mathcal{O}\left(10^{9}\right)$ photons. To model such large numbers of photons and long tail lengths we recourse to the linearization method described in Sec. III.4 along the lines demonstrated in the Appendix B.

This linearization technique is applied to the full master equation, Eq. (1). Since in the multi-mode model the Lindblad operators are only first-order in $a,b$ and $c_{j}$, only the third- and fourth- order linearization approximations are required. Therefore, we may reasonably expect the approximation to be more accurate in the multi-mode case than in the single-mode case, where linearization on sixth-order expectation values was required. We thus derive a closed system of equations for expectations up to second-order. An example of the evolution of both $Q$ and $\langle n_{-}\rangle$ is shown in Fig. 9.

When a coherent state containing $1.2\times 10^{9}$ photons is initialized in mode $s_{-}$ (solid red lines in Fig. 9), the system quickly evolves to a strongly sub-Poissonian state, even in the presence of a realistic loss rate $\gamma_{1}$. A brighter initial state results in larger $|Q|$ over shorter timescales, and in greater robustness of the sub-Poissonian output to linear loss $\gamma_{1}$.

For no linear loss, $\gamma_{1}=0$ (dashed line in Fig. 9, bottom), $Q\approx-0.8$ remains steady until a later time which is proportional to $N$. Thus, the tail effectively acts as a Markovian reservoir, which justifies its adiabatic elimination in Eq. (4) natcom . Crucially, therefore, the long tail of modes combined with the self-Kerr interaction corresponds to the effects of nonlinear coherent loss into a Markovian reservoir (compare Fig. 3 and Fig. 9), and we note that the timescales over which this correspondence holds may be increased by increasing the tail length $N$.

Since $s_{-}$ is a linear combination of modes $a$ and $b$, Eq. (5), the input coherent state with amplitude $\alpha$ into $s_{-}$ may be generated by initializing modes $a$ and $b$ with coherent states of appropriate amplitude. Alternatively, we may consider the initial condition $\langle a\rangle=\alpha$ and all other modes as vacuum. In this case, the behaviour of $s_{-}$ is qualitatively equivalent to Fig. 9, but a slightly different optimal ratio $g_{b}/g_{a}$ should be considered.

Finally, we note that in the case of $\gamma_{1}=0$, the signal modes are highly sensitive to whether $N$ is odd or even, and for even $N$ the decay of mode $s_{-}$ is inhibited. Indeed, the decay rate $\gamma^{\prime}$ into an effective reservoir changes depending on $N$, and may be understood as an interference effect resulting from the unitary coupling between tail modes affecting whether the excitation returns to signal mode $a$ or mode $b$ from the tail. This effect of the reservoir, which quickly vanishes at even small levels of linear loss $\gamma_{1}$, will not manifest itself over the timescales of interest provided that $N$ is large.

We have demonstrated that the multi-mode model – which is the model most closely aligned to a physical implementation via waveguides inscribed in bulk glass, Fig. 2(a) – can generate a strongly sub-Poissonian output from a classical input state into mode $s_{-}$ over the short timescales of initial evolution. This corroborates results from Sec. III and the prediction of optimal coupling parameters $g_{a,b}$ (12) when mode $s_{-}$ is initially excited.

To assess the impact of specific effects of short-pulse propagation in the system of coupled waveguides, we dynamically evolved the spectral and temporal properties of the pulse and analysed their influence on the effectiveness of the NCL mechanism. In particular, we considered the combined effect of chromatic dispersion, self-phase modulation (Kerr effect), and self-steepening of the pulse via the Nonlinear Schrödinger Equation (NLSE) for the system of waveguides in a $\chi^{\left(3\right)}$ nonlinear medium. The resulting system of coupled propagation equations for all waveguides in the PhoG structure was numerically solved using the split-step Fourier method Agrawal . For a specific comparison we used the signature of the NCL mechanism, the nonlinear decay behaviour of photon number $\langle n_{-}\rangle$ in the anti-symmetric mode $s_{-}$ displayed in the inset of Fig. 9 (a), which may additionally be interpreted as an intensity-dependent loss. In Fig. 10 such decay can be clearly seen when evolving an initially Gaussian pulse according to the NLSE. Linear loss was neglected in this simulation to ensure that the decay was caused by the NCL mechanism, though we note that even when linear loss is included we observe strong agreement between the nonlinear decay behaviour of the full pulse propagation (NLSE) numerics and the quantum multi-mode model considered in Fig. 9.

The pulse propagation was modelled for a $100$ fs initially Gaussian pulse in IG$2$ glass with a central wavelength of $1550$ nm (see Ref.ig2 , and also Sec. III.3.1). Initial pulse energies of $50$, $100$, $150$ and $200$ pJ were chosen, in order to compare with the inset of Fig. 9 (a). The nonlinear parameter $\bar{\gamma}^{\rm NL}$ was set at $0.6~{}\textrm{W}^{-1}\textrm{m}^{-1}$, considering that the pulse energy was $200$ pJ. The resulting evolution is displayed in Fig. 10 and we observe the characteristic nonlinear decay, which is a good evidence that NCL occurs for these realistic physical parameters even when the full pulse dynamics is considered. Linear loss was neglected in this simulation to ensure that the decay was caused by NCL mechanism.

It was found that $g_{a}=210\textrm{m}^{-1}$, $g_{b}=360\textrm{m}^{-1}$ and $g_{c}=430\textrm{m}^{-1}$ best reproduces the signature behaviour of NCL (different decaying rates for the different initial energies) as shown in Fig. 10. The optimal ratio between interaction constants is only weakly influenced by the effects of the pulse propagation. The spectral broadening of the pulse is estimated to be about 10nm for 100 fs input pulse and 2-3 cm device at 1550 nm central wavelength. Provided that the device length is sufficient for the shortest wavelength of the pulse to undergo NCL decay (and for the example shown in Fig. 10 it is sufficient), the effects of such pulse spectral broadening do not affect the manifestations of the NCL mechanism.

VII Conclusions

We have discussed a design of a family of devices PhoG based of coupled single-mode waveguides in bulk nonlinear glass. We have shown that an asymmetric structure can indeed function as a deterministic generator of bright sub-Poissonian states of light. The sub-Poissonian light generation with considerably large photon number squeezing (up to the values of the Mandel parameter of about -0.8) can occur at the initial stages of the dynamics. In this regime conventional linear loss (which is quite high in nonlinear glass) still can be overcome by sufficiently intensive coherent input, and the nonclassical states can be produced. A hierarchy of models has been derived, from the original multi-mode waveguide network down to the simplest single-mode model exhibiting both two-photon- and nonlinear-coherent- loss. We have developed an analytic approach to the single-mode model and discovered two dimensionless parameters allowing one to design the waveguide system with aim to reach a maximally nonclassical state for a given length of the structure. Feasibility of the suggested device has been analyzed for three different eligible practical structures: a system of waveguides in bulk glass, multi-core fiber, a set of coupled nano-wires. We have estimated the number of photons required to reach large nonclassicality with the femtosecond input pulses and demonstrated that this number remains well within the operational possibilities. Further, we have developed an approximation linearized on the quantum correction to the classical solution and demonstrated its closeness with the exact solution of the single-mode problem. Both the three-mode approximation and the complete multi-mode model lead to very close results in the initial stage of the dynamics.

In summary, suggested CDP circuits indeed allow for implemention of the engineered nonlinear loss to produce strongly non-classical states in realistic photonic structures. Our design of the PhoG with a system of waveguides laser-written in the bulk glass can serve as a practical recipe for realizing such structures. Being comparatively inexpensive and easy to produce, these generators can find applications for different detection, measurement and metrology tasks bio ; imaging and as quasi-single photon sources in various quantum technology applications, such as, e. g., quantum key distribution wang1 .

Acknowledgements

The authors acknowledge support from the EU Flagship on Quantum Technologies, project PhoG (820365). D.M., A.S. and A. M. also acknowledge support from the National Academy of Sciences of Belarus program “Convergence”, and the BRRFI project F18U-006. We thank Benjamin Brecht for the help in developing the code for pulse propagation in waveguide network.

Appendix A On Gaussianity of the photon-number distribution of the PhoG generated state

Here we derive approximate analytical expressions for dynamics of mean photon number, Mandel parameter and higher-order moments of the photon number distribution for the single-mode model. We start from the master equation (8) of the main text. For the photon number distribution, $p(n)=\langle n|\rho|n\rangle$, Eq. (8) yields

[TABLE]

where $\xi_{i}$ is the number of the photons lost per each dissipation event ($\xi_{1}=\xi_{3}=1$, $\xi_{2}=2$); $f_{1}(n)=n$, $f_{2}(n)=n(n-1)$, and $f_{3}(n)=n(n-1)^{2}$.

For a function $g(n)$, one can define its average as

[TABLE]

Eq. (24) implies that the average satisfies the following equation:

[TABLE]

For example, the mean photon number $\mu=\langle n\rangle$ corresponds to $g(n)=n$ and satisfies

[TABLE]

For further calculation of the photon number distribution moments, it is also useful to calculate the derivative for a function $g(\delta n)$ of the photon number deviation $\delta n=n-\mu$. Here, one needs to take into account that $\mu$ also varies with $z$:

[TABLE]

In further calculations, $n$ in Eq. (28) can also be expressed in terms of $\delta n$: $n=\mu+\delta n$.

Let us assume that the initial state is a coherent one with the mean photon number $\mu(z=0)=n_{0}$:

[TABLE]

Central moments of the input state are equal to

[TABLE]

Taking into account that Poisson distribution does not have heavy tails and its central moments are mainly determined by the central part of the distribution, one may expect that, at least for the beginning of the state evolution, $p(n)$ will have non-zero values for $|\delta n|\lesssim const\sqrt{\mu}$, and, therefore,

[TABLE]

Let us consider such $z$ that the conditions (31), together with

[TABLE]

are still satisfied. Substituting definition of $\zeta_{2}$ into Eq. (28) and using Eq. (31) for $k\geq 3$, one can derive the following equation:

[TABLE]

According to the derived equation, the quantity $\zeta_{2}$ tends to reach the value

[TABLE]

Taking into account that $\zeta_{2}(z=0)=1$, one can expect that the condition (31) will be satisfied during most of the system evolution, until the mean photon number $\mu$ becomes small.

Assuming that a single type of dissipation is prevailing, one can derive the values

[TABLE]

for linear, two-photon, and nonlinear coherent loss respectively. The corresponding values of Mandel parameter $Q=\zeta_{2}-1$ are

[TABLE]

Now let us assume that for the considered $z$ the conditions

[TABLE]

and

[TABLE]

are still satisfied (for $z=0$ one has $\zeta_{3}=1$ and $\zeta_{4}=1/\sqrt{\mu}$) together with Eq. (31) for $k\geq 5$.

Therefore, one can derive the following equation for $\zeta_{3}$:

[TABLE]

The quantity $\zeta_{3}$ tends to the value

[TABLE]

In a similar way one can show that $\zeta_{4}$ satisfies

[TABLE]

The derived equation has two important implications. First, it shows that the condition $\zeta_{4}=O(1)$ remains valid while $\langle\delta n^{5}\rangle/\mu^{5/2}=O(1)$, which in its turn is a consequence of Eqs. (33) and (34). Second, in contrast to lower-order moments, the equation for $\zeta_{4}$ is not disentangled from higher moments and cannot be solved independently of them. Still, the condition $\zeta_{4}=O(1)$ is sufficient for validity of the previously derived Eqs. (39) and (34).

For a single type of dissipation, one can derive the values

[TABLE]

for linear, two-photon, and nonlinear coherent loss respectively, which correspond to the following values of skewness $\zeta_{3}/(\sqrt{\mu}\zeta_{2}^{3/2})$:

[TABLE]

Similarly, for excess kurtosis one can derive the estimate

[TABLE]

The results show that the photon number distribution, despite becoming strongly sub-Poissonian, remain Gaussian with high accuracy until the mean photon number $\mu$ becomes small.

Appendix B Derivation of linearization approximations

In this section we will derive the linearization approximations which are used in Sec. III.4 and Sec. VI. These approximations become increasingly necessary for modelling the system when large photon numbers or large numbers of modes must be considered.

Consider arbitrary quantum operators $A,B,C$. We will explicitly show the derivation for linearization of $\langle ABC\rangle$, but the approach may be readily generalized to higher-order products. We will derive replacements which will allow us to approximate expectations of products of three or more operators, using only first- and second- order terms e.g. $\langle A\rangle$, $\langle AB\rangle$.

We may expand each operator $A,B,C$ into a classical “mean-field” term and a quantum fluctuation term: $A=\langle A\rangle+\delta A$. We take $\langle\delta A\rangle=0$ which ensures that the mean-field term $\langle A\rangle$ is meaningful and consistent. Then substituting these expansions into $\langle ABC\rangle$,

[TABLE]

A key tool we require is the cumulant expansion for generic operators $\left\{O_{1},\dots,O_{n}\right\}$

[TABLE]

where $\mathbb{P}$ denotes all disjoint partitions of the set of operators, $|\mathcal{P}|$ denotes the number of blocks in partition $\mathcal{P}$, and $p$ iterates over each block in the partition. For example

[TABLE]

We perform our linearization assumption on Eq. B by setting $\mathcal{C}\left(\delta A,\delta B,\delta C\right)=0$, which implies therefore that $\langle\delta A\delta B\delta C\rangle=0$ since $\langle\delta A\rangle=\langle\delta B\rangle=\langle\delta C\rangle=0$. Finally, using $\delta A=A-\langle A\rangle$ we arrive at our final expression,

[TABLE]

Higher-order expectations of products of operators may be calculated in the same way, with the only requirements assumed about fluctuations being the zero-mean condition $\langle\delta A\rangle=\dots=\langle\delta Z\rangle=0$ and the linearization approximation $\mathcal{C}\left(\delta A,\dots,\delta Z\right)=0$.

Appendix C Linearized single-mode model

By applying the linearization approximations derived in Appendix B to the system of coupled ODEs (III.4) we arrive at the following closed system of ODEs:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $n_{-}=s_{-}^{\dagger}s_{-}$.

This system is solved numerically for $\langle s_{-}\rangle$, $\langle s_{-}^{\dagger}\rangle$, $\langle s_{-}^{2}\rangle$, $\langle s_{-}^{\dagger 2}\rangle$, $\langle n_{-}\rangle$ and the results are shown as dashed lines in Fig. 5.