Continuous variable multimode quantum states via symmetric group velocity matching

Configurable and scalable continuous variable quantum networks for measurement-based quantum information protocols or multipartite quantum communication schemes can be obtained via parametric down conversion (PDC) in non-linear waveguides. In this work, we exploit symmetric group velocity matching (SGVM) to engineer the properties of the squeezed modes of the PDC. We identify type II PDC in a single waveguide as the best suited process, since multiple modes with non-negligible amount of squeezing can be obtained. We explore, for the first time, the waveguide dimensions, usually only set to ensure single-mode guiding, as an additional design parameter ensuring indistinguishability of the signal and idler fields. We investigate here potassium titanyl phosphate (KTP), which offers SGVM at telecommunications wavelengths, but our approach can be applied to any non-linear material and pump wavelength. This work paves the way towards the engineering of future large-scale quantum networks in the continuous variable regime.


I. INTRODUCTION
Photonic quantum technologies are grounded on the remarkable ability of engineering quantum states of light in different non-linear optical processes [1,2]. Besides the versatility in producing different quantum states, optics has the advantage of offering scalability in the number of systems by exploiting multimode electro-magnetic fields. This is even more relevant when we consider a continuous variable (CV) encoding of quantum information since it is possible to deterministically generate large multimode entangled states, where entanglement is established between amplitude and phase quadratures of different light modes [1,[3][4][5][6]. CV multimode entangled states are resources for measurement based quantum computing [7,8], quantum simulations [9], multi-party quantum communication [5,10], and quantum metrology [11,12].
For all the applications that rely on the transfer of information over appreciable distances, telecommunication wavelengths offer the most reliable solution. In this regime single-mode quantum states with large amount of squeezing [13,14], and on-chip few-modes squeezed states have been reported [15,16]. Second-order non-linear waveguides are promising candidates for the controlled generation of both single-and multimode quantum states at telecom wavelengths [17][18][19][20], but a direct characterization of CV multimode quantum resources, like cluster states, in these systems is still missing. Waveguides offer in fact large nonlinearities and long interaction lengths, along with the possibility to tailor the modal decomposition of the quantum state [21], but losses and modal engineering have to be carefully controlled in order to measure the desired quantum effects via homodyne detection used in CV protocols.
In this work we study the generation of CV multimode entangled states at telecom wavelengths via parametric down conversion (PDC) in non-linear waveguides. We aim at obtaining a controlled number of modes with sufficient squeezing, which can coherently interfere in order to generate the wanted entangled states via schemes exploiting both the spectral [1] and the temporal [4] degrees of freedom of a femtosecond laser source [22]. We identify type II PDC in periodically poled non-linear waveguides in the so-called symmetric group velocity matching (SGVM) configuration [23] as the most promising scheme.
The paper is structured as follows. In Section II, we set the theoretical grounds, define the physical variables and discuss the optimal configuration of the non-linear waveguide. In Section III, we revisit the concepts underlying SGVM and show the properties of type 0, I and II PDC processes, concluding that type II is best suited for the target application. Section IV applies the results to potassium titanyl phosphate (KTP) and yields realistic fabrication parameters. Section V finishes with the conclusions.

A. Squeezed modes
The parametric interaction in a χ (2) non-linear crystal involving a pump field and two output fields, usually called signal and idler, can be expressed in appropriate conditions (see for example [24]) via an effective hamilto-nianĤ that governs the evolution of the quantum operators in the interaction picture. This effective hamiltonian can be written (assuming a classical undepleted pump), where J i,j are the elements of the so-called joint twophoton amplitude andâ † i andb † j are creation operators of, respectively, the signal and idler photons in spatiotemporal modes f i (r, t) and l j (r, t), each forming an orthonormal mode basis.
The discrete form of the hamiltonian in Eq. (1) is due to the implicit discrete mode basis {f i (r, t), l j (r, t)} in which it is expressed. In this work, we will restrict ourselves to temporal modes, assuming a single spatial mode common for all the fields involved, a configuration that is typically realised in waveguides. Under these conditions we write f i (r, t) → f i (t) where ω is the frequency and F.T. stands for Fourier Transform.
In the continuous frequency basis, the effective hamiltonian of Eq. (1) is written as: where now the joint two-photon amplitude is a twodimensional function of the signal and idler frequencies ω s and ω i , known as the joint spectral amplitude [25], and the creation operators create field excitations at different frequencies.
The joint spectral amplitude, J(ω s , ω i ), contains information about the pump field that enters the waveguide and the linear and non-linear optical properties of the waveguide itself. One can perform a decomposition of this function in a set of modal functions {h k (ω s )} and {g k (ω i )} [26]: This decomposition is also known as the Schmidt decomposition, and the set of complex values {λ k } the Schmidt coefficients, fulfilling k |λ k | 2 = 1. By inserting Eq. (3) into Eq. (2) we obtain: The Hamiltonian of Eq. (4) generates two sets of eigenmodes that are EPR entangled with each other (twin modes), where each couple is in the so called two-mode squeezed state [24,27,28], and which are defined spectrally by the same sets {h k (ω s )} and {g k (ω i )} of the Schmidt decomposition in Eq. (3). These eigenmodes are also called temporal modes [19]. Eq. (5) is then simply a basis change from the continuous frequency basis to the discrete Schmidt basis.
From the continuous variable picture of the process, we define the normalized quadrature operators associated to the temporal modes as: The EPR couples are then characterized by the following relations of quadrature variances are the variances of the quadratures in the case of vacuum states, that here are taken equal to one according to the quadrature definition in Eq. (6). This means that in the case of high gain, i.e., high value for the Schmidt coefficients λ k , the quantum state of each EPR couple is an approximation of the so called quadrature entangled EPR state [29], characterized byq A,k =q B,k and p A,k = −p B,k . It can be shown [30] that the effective number of modes, K, i.e., the effective number of pairs of functions from {h k (ω s ), g k (ω i )}, can be calculated as This quantity is often called the Schmidt number, and it will be a figure of merit in our analysis. The temporal modes form eigenbases of the signal and idler parts of the total field. One can show [24] that the eigenvalues of the total joint spectral amplitude are doubly degenerate, with eigenvectors {s k (ω)} such that: which are mutually uncorrelated modes that present squeezing. These modes are called squeezed modes (or "supermodes" [31]), and their creation operators are de- . The effective hamiltonian in terms of the supermode creation operators is then: The corresponding variance of the quadrature variables, with operators {q S,k ,p S,k }, defined in the same way as in Eq. (6), are ∆ 2 (q S,k ) = e −2λ k and ∆ 2 (p S,k ) = e 2λ k . They are squeezed, and their level of squeezing, 2λ k , is proportional to the corresponding Schmidt eigenvalue λ k . Hence, in the temporal mode basis, the signal and idler parts of the total field are individually diagonalized and there exists entanglement between them, while in the squeezed mode basis, the modes are uncorrelated and carry squeezing as a quantum resource.
It is important to note that the transformation of Eq. (8) physically corresponds to the interference of the signal and idler modes on a balanced beamsplitter.
In the following, we focus on the calculation and characterization of the squeezed eigenmodes s k (ω) and their squeezing eigenvalues λ k , since they constitute a natural multimode quantum state where continuous variable quantum information can be encoded and processed. In particular, the Bloch-Messiah reduction [32] ensures that all multimode Gaussian states, i.e., the states whose quadratures are characterized by Gaussian statistics and that can be generated by hamiltonians that are at most quadratic in the bosonic operators like the ones in Eqs.
(1), (2) (4) and (9), can be written as a set of single mode squeezed states and two multiport interferometers. In our context, this means that we could generate deterministically any highly entangled multimode Gaussian state if we input the squeezed eigenmodes into an appropriate multiport interferometer, generating a configurable source of cluster states that could then be used as the resource for measurement based quantum computation, multiparty quantum communication or quantum metrology.

B. Use of waveguides and mode detection
Until now we have considered the general case of nondegenerate parametric down conversion and the emergence of the squeezed modes. We will now specify and justify the configuration that will be treated in the rest of the paper.
We will consider the use of non-linear waveguides, instead of bulk crystals, for their ability to confine light over a long length (in particular over the coherence length of the pump laser) and their discretization in the number of available spatial modes. Both effects enhance the nonlinear interaction by several orders of magnitude. This enables the generation of single-pass squeezing [21]; no optical cavities are required. The use of waveguides also means that we will work in the so-called collinear configuration, where signal and idler fields exit the crystal in the same beam.
Another important reason for the use of waveguides is that, as we will see later, the engineering of the physical parameters of the waveguide opens up new possibilities for engineering the PDC, and hence the generated modes, which will benefit our target applications.
Furthermore, the measurement of the modes would be performed via coherent detection, in particular homodyne detection, where the quadratures of the targeted field can be accessed by interfering the field with a local oscillator. Hence, in order to measure the quadratures of the different squeezed modes, the local oscillator should be shaped spectrally to match the corresponding modes {s k (ω)}. However, this becomes a difficult task if the squeezed modes take arbitrary shapes. In Eq. (8) we mentioned that the squeezed modes are obtained by the interference of signal and idler temporal modes, which can themselves be approximated as Hermite Gauss modes [20]. In order for this interference to be efficient and for the squeezed modes to inherit the approximate Hermite Gauss functional form, both signal and idler temporal modes should be as similar as possible, which translates in having indistinguishable fields after the non-linear interaction in the waveguide.
The spatial indistinguishability is given due to the collinear configuration imposed by the use of waveguides, while in the bulk crystal scenario a precise alignment of the crystal would be needed. The spectral indistinguishability constrains to work in the so-called degenerate case, where both central frequencies for signal and idler are equal (and equal to half the pump central frequency by conservation of energy): ω 0 s = ω 0 The final degree of freedom of the fields is their polarization. Depending on the output polarization in the PDC process, three types of parametric interactions can be found. In type 0 and type I, the polarization of signal and idler fields is the same, while in type II the fields have orthogonal polarizations, making them distinguishable. However, we will show in the next section that, due to their dispersion properties, both type 0 and type I prove unpractical. We will consequently restrict ourselves to type II processes. We note that the polarization distinguishability can be erased by rotating the polarizations of the signal and idler by 45 degrees and interfering them on a polarizing beamsplitter. At the output of the beamsplitter, one then obtains the squeezed modes [21]. We also note here that, due to the different polarizations of the signal and idler modes, there may be a temporal walk-off after the waveguide, which can be corrected for with appropriate compensation crystals. Fig. 1 schematically shows the configuration discussed above: degenerate collinear type II parametric down conversion in a non-linear waveguide.

III. SYMMETRIC GROUP-VELOCITY MATCHING
A. Joint spectral amplitude structure The mathematical structure of the joint spectral amplitude, or JSA, of Eq. (2), J(ω s , ω i ), can be expressed as: where the subscripts s/i/p stand for signal, idler, and pump, and ω p = ω s + ω i . The function α p (ω s + ω i ) is the pump envelope function in the spectral domain, which we take to be a Gaussian with a certain width, denoted by w p (see Fig. 1). Finally, φ(ω s , ω i ) is the phasematching function. It is governed by the waveguide's optical properties and the geometry of the problem. This function does not generally have properties of symmetry. The phasematching function can be written as: where L is the waveguide length and the function ∆k(ω s , ω i ) is the wavevector mismatch where the k are the corresponding wavevectors. In general, Eq. (12) is a vectorial equation, but we recover its scalar version because of the collinear condition. It is possible to match two arbitrary signal and idler frequencies so that ∆k = 0 with a well-known technique known as quasi-phasematching [33], which introduces a new term to Eq. (12) by modulating the χ (2) coefficient along the waveguide with a poling period Λ. Eq. (12) including the quasi-phasematching modulation can be written as: such that the selection of the poling period makes the mismatch zero for a certain pair of desired frequencies.
In this work, the poling period will correspond to the case of matched signal/idler central frequencies.
B. Type 0 and I PDC Indistinguishability of the signal and idler fields is required for their proper interference after the non-linear interaction. This implies that the joint spectral amplitude has to be symmetric under the exchange of the signal and idler labels, i.e., Since the pump function α p (ω i + ω s ) is naturally symmetric under this exchange, the fields are indistinguishable if the phasematching function shares this symmetry, i.e., if φ(ω s , ω i ) = φ(ω i , ω s ). From Eq. (11) we therefore require the condition on the wavevector mismatch: ∆k(ω s , ω i ) = ∆k(ω i , ω s ).
In the case of type 0 and type I PDC, it is easy to see that the symmetry condition is fulfilled, since k s (ω s ) = k i (ω i ) in these cases.
As in [23], we write the Taylor series of the mismatch function around the central frequencies ω 0 s = ω 0 i = ω 0 p /2 to second order: The coefficients are related to the derivatives of the wavevector with respect to the frequency as: where we wrote ω p instead of ω 0 p in the evaluation of the derivatives for clarity in the notation. The constant term ∆k 0 = ∆k(ω 0 s , ω 0 i ) can be made zero via quasiphasematching, as described previously.
The coefficient γ s (γ i ) associated with the linear term in the expansion is related to the pump and signal (idler) group velocities, while the coefficients δ s (δ i ) associated with the quadratic term is related to the pump and signal (idler) group velocity dispersion.
In type 0 and type I PDC, γ i = γ s ≡ γ, δ i = δ s ≡ δ, etc. In addition, usually w p ≤ 10 nm, such that the linear term dominates. Hence, we can approximately write: The phasematching can be approximated by a function of ω s + ω i , and therefore so can the joint spectral amplitude Since both the pump envelope and the phasematching are functions of the sum of the frequencies, they are aligned, making a −45 • angle, in the 2-dimensional coordinate system with axes (ω s ,ω i ). The JSA is then non zero in a large spectral region, yielding spectrally broad temporal modes that are not well-suited for homodyne detection, since one needs good spectral overlap with a local oscillator. Fig. 2 a) shows this situation. Overcoming this problem necessitates either very broad local oscillators and/or higher order term contribution in the mismatch function, both of which is unpractical.
However, looking at Eq. (15), we obtain identical modes if we manage to have the condition γ s = −γ i , which is called the symmetric group velocity matching (SGVM) condition [23]. This translates into having the following condition on the group velocities (recalling that the group velocity is v g (ω) = ∂ω/∂k): In that case, the Taylor expansion of Eq. (15) gives: with γ ≡ γ s = −γ i . Hence, the mismatch function is antisymmetric under the exchange of signal and idler labels: ∆k(ω s , ω i ) = −∆k(ω i , ω s ). By Eq. (11) and the even parity of the sinc function we obtain a symmetric phasematching function, and hence a symmetric joint spectral amplitude (up to first order in the Taylor series of Eq. (15)). Under this condition, signal and idler fields are indistinguishable (except for their polarization). From Eq. (19), the joint spectral amplitude is approximately of the form In Fig. 2 b) we see that the phasematching function now makes an angle of +45 • in the (ω s ,ω i ) coordinate system, hence making a 90 • -angle with the pump envelope function. Therefore, the JSA is non zero only in a small spectral region around the central frequencies.
In this situation, the coherent detection of the squeezed modes is experimentally realizable and the interference of signal and idler fields is approximately perfect after the change of polarization of one of the two.
In conclusion, all the reasoning exposed until now leads us to consider the configuration of degenerate, collinear, type II PDC in non-linear waveguides. We will also consider being around the symmetric group velocity matching condition, which means: • The poling period Λ provides quasi-phasematching for the central frequencies ω 0 p and ω 0 s = ω 0 i = ω 0 p /2. • The inverse group velocities at the central frequencies satisfy the condition γ s = −γ i , where these coefficients are defined in Eq. (16).
• The higher order terms in the Taylor expansion of Eq. (15) are sufficiently small compared to the first order term in the entire frequency range considered (see B for insights into this point).
We will also call the wavelength fulfilling these conditions as symmetric group velocity matching wavelength λ SGVM .

D. Index of refraction
The dispersion relation where c is the speed of light in vacuum and n the index of refraction, encapsulates the optical properties of the crystal and could depend on several physical variables apart from frequency, such as temperature, T , and waveguide characteristics. The latter include the waveguide section size, with height h and width w, and spatial mode order inside the waveguide, defined by two integer numbers n 1 and n 2 . If one is interested in modeling the functional form of the index of refraction with these waveguide parameters, approximations must be made. A wide-spread approximation is the so-called metallic waveguide approximation [34], in which one assumes the waveguide to be surrounded by perfectly conducting edges. This is the approximation that will be made in this work. More refined modeling of the index of refraction involves more advanced computational methods solving Maxwell equations inside and outside the waveguide, like Mercatili's method [35] or finite element methods [36].
Using the metallic waveguide approximation, the index of refraction is written as n(λ, T, n 1 , n 2 , w, h) = n(λ, T )+ λ n 1 + 1 2h where the function n(λ, T ) depends on the non-linear material under consideration, and is given by the empirical Sellmeier equations (for e.g. KTP see [37]).
In this work, we will restrict ourselves to the case in which only the fundamental spatial mode travels through the waveguide, which translates into n 1 = n 2 = 0. This is exactly the case in a perfect single-mode waveguide and, for shorter wavelengths, this can be achieved by careful mode-matching in the input coupling.
Furthermore, we will have a different index of refraction of the form of Eq. (22) for different light polarizations. For instance, in the case of uniaxial crystals, we would have two indices of refraction (for the ordinary and extraordinary axis), and in biaxial crystals, like KTP, we would have one index of refraction for every spatial direction.
These indices define the different wavevectors for pump, signal and idler fields in each case. We can therefore study the symmetric group velocity matching condition via the dependence of the indices of refraction on wavelength, temperature and waveguide dimensions.

E. Mode overlap
To wrap up this section, we define the quantity that will quantify the similarity between signal and idler temporal modes as the overlap integral o n . The n-th order overlap between two modal functions, h n (ω) and g n (ω), is defined as: where N is a normalization factor, such that the overlap is adimensional and its maximum value is 1, if and only if h n (ω) = g n (ω).

IV. RESULTS
In our analysis, we will focus on KTP as the nonlinear material for the waveguide. Lithium niobate (LN) is another important non-linear material, since commercial LN waveguides already exist. However, calculations with LN in our configuration lead to the conclusion that the waveguide width and height should be in the order of 1 µm for the symmetric group velocity condition to hold. This is smaller than the telecom wavelengths considered here (around 1.55 µm). The metallic waveguide approximation breaks down at this scale and more elaborate techniques must be used in order to model the index of refraction of these small structures. Therefore, we will only consider KTP in this work, which naturally provides SGVM at telecom wavelenghts.

A. SGVM wavelength
First, we performed calculations to compute the wavelengths fulfilling the SGVM condition. This first section shows the existence of this wavelength for KTP for reasonable physical parameters.
For this, we fixed the temperature and the waveguide dimensions to typical values and computed the dispersion properties of the crystal, in particular the different group velocities and group dispersions in KTP, from the indices of refraction of Eq. (22), as a function of wavelength. Thanks to the metallic waveguide approximation, the derivatives can be computed analytically in this case. In this way the coefficients of Eq. (16), as well as the corresponding poling period, are calculated as a function of wavelength. Fig. 3 shows the Taylor coefficients associated with the group velocities (γ s and −γ i ). They coincide at the SGVM phasematching wavelength, λ SGVM , where γ s = −γ i .
Once λ SGVM is obtained, it is important to check the contribution of the higher order terms in Eq. (15). This is discussed in B.

B. Telecom wavelength
In the calculations above, a wavelength fulfilling the SGVM condition was proven to exist for KTP, Fig. 3. We will now be interested in varying the other physical variables affecting the index of refraction in Eq. (22), namely temperature and waveguide dimensions, to get the condition at a certain desired wavelength.
In particular, we are interested in obtaining SGVM at 1550 nm, given that telecom wavelengths present low losses in optical fibers, and there exist optimized of-theshelf components at this wavelength, opening the way to practical applications.
It is important to remark the geometrical choice of the problem, because results would be different depending on the orientation of the crystal's optical axes with respect to the pump field polarization. In particular, the Kleinman symmetry [38] provides the supported processes for all combinations in propagation and polarization directions. The only configuration in which we have results at wavelengths around 1550 nm is the one that corresponds to a crystal cut (i.e., a vertical polarization direction) in the z direction and propagation direction along the x direction, where the indices of refraction for every spatial direction are ordered like in [37]. Therefore this is the case considered in this work. More details can be found in A.
As announced above, we will now consider the dependence of the index of refraction on temperature and waveguide dimensions.
a. Temperature: Fig. 4 shows the temperature dependence of the SGVM wavelength for KTP. For this calculation we fixed the waveguide dimensions such that λ SGVM 1550 nm.
We obtain a change in the SGVM wavelength of ∆λ SGV M 4 nm in a range of temperatures of ∆T 230 o C. This result indicates that temperature has a very limited effect on the change of λ SGVM .
According to these results, it seems reasonable to fix the temperature and focus on the tuning of the waveguide dimensions. We therefore fixed the temperature to be room temperature for the rest of the work, which is natural when considering future technology applications, where involved temperature control may be undesirable.
b. Waveguide Dimensions: Fig. 5 shows the SGVM wavelength as a function of the waveguide width and height for KTP. As the waveguide size increases, the surface flattens, making apparent the diminishing impact of the waveguide dispersion. We also find this behaviour in Eq. (22), where the waveguide dimensions contribution to the index of refraction disappears as h, w → ∞. The green area in Fig. 5 marks values of λ SGVM that lie in the interval between 1549 and 1551 nm. Therefore, we conclude that to work at 1550 nm in KTP, the waveguide should have a size of about 9 by 9 µm in width and height. From Fig. 4 we know that this is almost independent of temperature in the temperature interval considered.

C. Engineering of the quantum states
Next, we performed simulations of the PDC process for KTP to compute the signal and idler modes of our system with the associated squeezing eigenvalues under the SGVM condition. We remark that this exact procedure could be done for any non-linear optical material if its dispersion properties are known.
Our pump function is modeled as a (squarenormalized) Gaussian function in the following way: where the parameter w p is the pump width. The pump central wavelength is λ 0 p = 775 nm, corresponding to a central frequency ω 0 p = 2πc/λ 0 p . Furthermore, in our results, we will express the pump width in length units (wavelength bandwidth instead of frequency bandwidth), by redefining w p as 2πc/w p .
The phasematching function is that of Eq. (11), where we fix the waveguide length. The wavevector mismatch is calculated using Eqs. (21) and (22), where we also set values for the temperature and waveguide dimensions in With these two functions we compute the joint spectral amplitude, or JSA, according to Eq. (10).
We then perform the Schmidt decomposition of Eq. (3) by computing the singular value decomposition of the JSA, obtaining two matrices containing the signal and idler temporal modes and a diagonal matrix with the Schmidt coefficients {λ k }. We recall that the distribution of the squeezing eigenvalues is equivalent (except for a global factor of 2) to the distribution of the Schmidt coefficients. From these coefficients we obtain the Schmidt number according to Eq. (7) and we also compute the temporal modes' full width at half maximum (FWHM). This yields the spectral width of the temporal modes, which is important for designing their coherent detection. Finally, we use Eq. (23) to obtain the overlap between the signal and idler modes in order to quantify their similarity.
The main goal of this analysis is to obtain the dependence of the quantum fields on input parameters that will be specified in the following. For the purpose of this work, we can see the computation just described as a set of input variables that are experimentally controllable, and output variables after the parametric interaction. The main input variables are the waveguide length, the waveguide dimensions (height and width), and the pump width, and the main output variables are the Schmidt number, the modes FWHM and their overlap. We then perform distinct simulations to obtain the corresponding dependencies. Note that the temperature could be taken as another input variable. However, following the discussion earlier, we consider it to be a fixed variable (at room temperature), since we have seen that it would have a minor effect on the output modes.
For the sake of clarity, we construct Table I, where each column represents an input variable and each row a corresponding output variable. The table elements show general dependencies that will be discussed in detail in the following.
The range of values adopted for the input variables is the following: we consider pump widths from 1 to 12 nm, waveguide dimensions between 4 and 10 µm, and waveguide lengths between 5 and 30 mm. The pump width can be modified by tuning the pump field that enters the waveguide, while the waveguide dimensions and length are fixed during the waveguide fabrication process.
The rest of this section will explore Table I in detail. As a general comment, the dependencies are based on the following observations: • The waveguide length, L, changes the bandwidth of the phasematching function as sinc L −1 .
• The width of the pump, w p , although selfexplicative, changes the width of the pump function.
• The waveguide dimensions mainly change the angle between the pump function and the phasematching function, because they change the optical properties (index of refraction) of the material. This is why they are the critical parameter in order to obtain the SGVM condition and why it is interesting considering them as an engineering parameter for the quantum states. We have also seen in this work that even though temperature also takes part in the index of refraction, its contribution is negligible in the usual experimental range of values. We remark that given fixed waveguide dimensions and temperature, the angle between the phasematching function and the pump function is totally defined by the choice of the non-linear material.
The results of this section can be understood keeping in mind these three observations, which are schematically shown in Fig. 6. For the sake of clarity, we will refer to the elements of Table I like matrix elements. For instance, the element 2.1 refers to the second row and first column, which corresponds to the dependence of the mode FWHM on the pump width.
First row of Table I, the Schmidt number: Fig. 7 shows the Schmidt number as a function of pump width for different waveguide dimensions (elements 1.1 and 1.2 of Table I). The Schmidt number, and thus the number of modes, grows linearly with the pump width, while the waveguide dimensions have almost no influence. The waveguide length was set to L = 10 mm for this simulation. This result can be interpreted in the following way: as the pump width increases, the pump envelope does too, and hence its intersection with the phasematching envelope. This creates frequency correlations between signal and idler fields and hence increases K. This is true because in this set-up the pump envelope is bigger at any time than the phasematching envelope. In general, there would be a minimum value of K when both envelopes have equal widths.  Table I. Fig. 8 shows the Schmidt number as a function of the waveguide length for different pump widths. In this case we are changing the phasematching envelope, and we observe the minimum value of K when both the pump envelope's and the phasematching envelope's widths are equal. This is why we see that the waveguide length in which the minimum value is reached depends on the pump width. Also, the minimum value itself would approximately equal 1 when the intersection of α p and φ was at 90 • that is, at the SGVM condition. This is interesting for applications in single photon sources, where K = 1 ensures maximum purity of heralded single pho-  Table I. tons.
After the minimum point, increasing the waveguide length decreases the width of the phasematching envelope, and so we are effectively in the same situation as in Fig. 7, and we observe the same type of behaviour, which is a linear increase in K. It is important to note that in Fig. 8 we see that one can control the slope of the increasing linear regime by controlling the pump width, which opens the possibility of modifying the number of modes just by adjusting the pump field. In a realistic scenario, the maximum waveguide length would be limited by the optical losses within the waveguide and the size of the available substrate material. State of the art KTP waveguides have losses of 0.25 dB/cm, which would point towards waveguide lengths in between 10 and 20 mm. This is in accordance with the maximum size of KTP substrates of roughly 30 × 30 mm 2 .
Equivalently, the Schmidt number increases rapidly by decreasing the waveguide length before the minimum, because in that case the 1/L dependence on the phasematching width is more important, creating strong wiggles in the sinc function in Eq. (11), and in turn producing strong correlations between the output fields, hence increasing K.
Second row of Table I, the modes bandwidth: We remark that, given the nature of the Hermite Gauss modes, the mode FWHM of interest is only the first mode (that is approximately Gaussian), as the rest would grow in width as √ n, where n is the mode order. Therefore knowing the first FWHM is sufficient to approximately characterize all of them. Fig. 9 shows the width of this mode as a function of pump width, which corresponds to the element 2.1 in Table I. The waveguide dimensions were set to w, h = 6 µm, intentionally away from the SGVM condition. We  Table I.
can see that indeed signal and idler modes have different widths for all pump widths and hence their interference would not be perfect with these waveguide dimensions, as expected from the results of the SGVM analysis. Fig.10 shows the modes' FWHM as a function of the waveguide dimensions for pump widths of 2 and 12 nm. Both signal and idler modes coincide at waveguide dimensions near 9 µm (we considered square waveguides), as expected from the SGVM analysis before. We show the plots for two different pump widths to remark again that the similarity of signal and idler modes is independent of the pump width, as it is related to the waveguide's phasematching function. These relations correspond to the element 2.2 in Table I.  Table I.
Finally, Fig. 11 shows the modes FWHM as a function of waveguide length for w p = 3 nm and w, h = 6 µm, again intentionally away from the SGVM condition. Both widths decrease with length, which would point towards the use of the longest possible waveguide if we are to measure the fields by coherent detection. This relation corresponds to element 2.3 in Table I.
Third row of Table I, the modes overlap: Fig. 12 shows the overlap of the first six temporal modes as a function of pump width, leaving fixed the rest of the variables. We can see that the overlap decreases with increasing mode order; the reason is discussed in Appendix C, in which we showed that this effect is due to two factors: the increasing width of the modes with the mode order, and the mathematical properties of the Hermite Gauss modes themselves. The pump width does not change the overlap, even though we have just seen that it changes the relative modes' FWHM by up to 25%. This counter-intuitive behaviour arises due to the stability of the relative width of two Gaussian functions when their relative widths change. The mathematical result is that for a relative width difference of 25%, i.e., if the width of one Gaussian is 1.25 times the width of the other, then the change in their overlap does not even reach 1%. This explains why we obtain a constant overlap for different pump widths, although the modes' widths change significantly.
In a similar way, Fig. 13 shows that the waveguide length has no effect on the mode overlap. The Schmidt numbers, and hence effective mode numbers, for different lengths are marked in the plot. We remark that the real number of non zero modes (that we compute with the singular value decomposition of the JSA), is always higher than the Schmidt number K. For example, in Fig. 12, at a pump width of 4 nm, we have K = 2.30 and 4 Schmidt modes, while when K = 3.54, at a pump width of 6 nm, we already have 6 Schmidt modes. It is therefore worth noting that the Schmidt number K is an indicative of the number of Schmidt modes present, but not their exact quantity [30]. In particular, if m is equal to the real number of Schmidt modes, then K ≤ m. These two figures correspond to the elements 3.1 and 3.3 in Table I.  Table I.  Table I.
Completing Table I, Fig. 14 shows the mode overlap as a function of the waveguide dimensions. We observe that the individual temporal mode overlaps have a common maximum at waveguide dimensions around 9 µm, as we computed in the SGVM study. Also, for every waveguide dimension, the similarity between modes decreases as the mode order increases (again the reason for this can be  Table I. found in Appendix C). This dependence corresponds to element 3.2 in Table I.

V. CONCLUSIONS
In this work, we have studied the generation of multimode squeezed states in single-pass PDC in non-linear optical waveguides. Such states are appealing candidates for the realization of numerous quantum information applications, including measurement based quantum computation and multipartite quantum communication. For long distances and integrated application, it is desirable to generate these states at telecom wavelengths.
We have first set the theoretical grounds for the PDC process in a non-linear waveguide, leading to the generation of uncorrelated squeezed states, the "supermodes", and we identify them as a natural platform for continuous variable quantum information science. Then, we explored the different PDC configurations.
In type 0 and type I PDC the generated signal and idler fields are indistinguishable. Nevertheless we have shown that these schemes are not practical for CV protocols. Although highly multimode states are generated, the generated modes have indeed a too large spectrum to be accessible via homodyne detection.
Consequently, we targeted type II PDC as the suitable scheme. In this case, the process generates a number of two-mode squeezed states, that, upon interference on a beamsplitter, can be turned into the single-mode squeezed resource states. For this to work efficiently, signal and idler fields have to be spectrally indistinguishable. We showed that this requirement enforces the use of SGVM in PDC.
We then investigated the use of KTP, a nonlinear material which naturally provides SGVM for wavelengths around the telecommunications regime. We studied the number of generated squeezed states and their similarity as function of the PDC pump spectral width and waveguide length. In addition, we explored the influence of the waveguide dimensions on the generated PDC state. Waveguide dimensions are usually only set to ensure spatially single-mode operations. However, it turns out that they play a crucial role for achieving spectrally indistinguishable signal and idler fields at a desired wavelength.
Our calculations reveal that a waveguide width and height of 9 µm facilitates SGVM at a signal and idler central wavelength of 1550 nm. In this configuration, both the waveguide length and pump spectral bandwidth can be used to tune the number of generated states, while having no detrimental impact on signal-idler similarity whatsoever.
In conclusion we have explored symmetric group velocity matching in non-linear waveguides for the tailored generation of scalable quantum resource for continuous variable based technologies. Moreover we have identified the suited processes and materials along with realistic parameters for their generation at telecommunication wavelengths.
However, we remark that with other polarization combinations or other non-linear materials, different wavelengths fulfilling the symmetric group velocity matching condition could be found even for the interesting case of waveguide dimensions compatible with spatial singlemode propagation.

Appendix B: Higher order contribution in the mismatch
We have seen that one of the conditions for being under SGVM is that the second and higher order terms of Eq. (15) should be negligible compared to the first order contribution. As we are performing a two-variable Taylor expansion, the higher order contributions are expected to get bigger as we get further away from the central frequencies ω 0 s = ω 0 i . Inserting Eq. (22) in Eq. (21) we can numerically compute the wavevectors for pump, signal and idler, hence allowing the calculation (up to numerical error), of the wavector mismatch ∆k, by Eq. (13), where the poling period Λ would be computed beforehand for phasematching the central frequencies.
The first order term, that we will denote here as F (ω s , ω i ), can be calculated as: where the coefficients have been defined in Eq. (16).
Therefore, the higher order terms, denoted as O(ω s , ω i ) can be computed by: In general, the higher order contributions depend on waveguide dimensions and temperature. In our case, we are interested in seeing the contribution of these terms under the phasematching conditions, i.e., when γ s = −γ i , so that our waveguide dimensions are fixed in order to satisfy the condition.
As we have seen in the paper, temperature does not play an important role in the wavelength where the first order term provides the phasematching condition. We have checked that temperature does not change appreciably the relative values of the higher order terms with respect to the first order term either.
Also, we will show here only the case of KTP, although the same procedure could be performed for any nonlinear material.
FIG. 15. Left to right: numerical mismatch, ∆k(λs, λi), first order term, F (λs, λi) and higher order terms, O(λs, λi) as a function of signal and idler wavelengths, λs, λi. Fig. 15 shows the three functions ∆k, F and O = ∆k − F under phasematching conditions for KTP as a function of the signal and idler wavelengths (recall that λ = 2π/ω). We can see that the mismatch function is effectively well described by the first order term in the wavelength region considered (that is bigger than the Schmidt modes bandwidth calculated in the paper).
We therefore conclude that, for KTP, and under our configuration, the higher order terms in the mismatch can be discarded in practice and the SGVM condition holds.

Analitical mode overlap calculation
In Fig.12, Fig.13 and Fig.14, we can see that the overlap between the modes decreases as the mode order increases.
The reason for this is two-folded: first, as the mode order increases, their spectral width does too, reaching frequencies further away from the central frequency, where the Taylor expansion of Eq.(15) is less precise and hence the modes are expected to present more differences between them. Their overlap is expected to be smaller, as we confirm.
The second reason resides in the mathematical structure of the modes itself. If we approximate the phasematching function (sinc function), by a gaussian function, the Schmidt modes can be calculated analytically, giving two sets of Hermite-Gauss modes of the form: HG n (x) = 1 n! √ π2 n w H n (x) exp(−(x − x 0 ) 2 /2w 2 ) (B3) where H n (x) is the n-th order Hermite polynomial and w is the associated width of the zeroth-order function, which is a gaussian. Our Schmidt modes, depicted in Fig.1, are not Hermite-Gauss, but they approximate them fairly well.
If we try to compute the n-th order overlap defined in Eq.(23) from two sets of Hermite-Gauss modes with different widths, w and w , we find, using Eq.(B3): (B4) and therefore the problem is reduced to finding the integral: for the solution of this integral one can define both the whole family of two variable Hermite polynomials H n (x, y) and the family of two-indices Hermite polynomials H n,m (x, y, w, z|τ ) [40]. By the method of the generating function, one finds that the solution to the integral L in Eq.(B5) is an evaluation of the two-index (with the same index n and n) Hermite polynomial in specific coordinates depending only on a and b. For typical values of relative widths for signal and idler Schmidt modes computed in our simulations, the overlap between modes decreases slightly and linearly. If we also model the fact that our modes are expected to be less similar as the mode order increases, due to their larger spectral width, by making their relative widths to increase with the mode order, then we obtain Fig.16. In this Figure, the relative widths for the firsts modes (gaussians) is set to w/w = 0.95. From this value, we decrease their relative widths by 2% for every higher mode order. We observe a decreasing overlap function with the mode order, in very good accordance with the results from the manuscript. We conclude that both the dispersion and the mathematical nature of the modes are responsible for the observed behaviour of the overlap as the mode order increases.