Properties of bright squeezed vacuum at increasing brightness

A bright squeezed vacuum (BSV) is a nonclassical macroscopic state of light, which is generated through high-gain parametric down-conversion or four-wave mixing. Although the BSV is an important tool in quantum optics and has a lot of applications, its theoretical description is still not complete. In particular, the existing description in terms of Schmidt modes with gain-independent shapes fails to explain the spectral broadening observed in the experiment as the mean number of photons increases. Meanwhile, the semiclassical description accounting for the broadening does not allow us to decouple the intermodal photon-number correlations. In this work, we present a new generalized theoretical approach to describe the spatial properties of a multimode BSV. In the multimode case, one has to take into account the complicated interplay between all involved modes: each plane-wave mode interacts with all other modes, which complicates the problem signiﬁcantly. The developed approach is based on exchanging the ( k , t ) and ( ω, z ) representations and solving a system of integrodifferential equations. Our approach predicts correctly the dynamics of the Schmidt modes and the broadening of the angular distribution with the increase in the BSV mean photon number due to a stronger pumping. Moreover, the model correctly describes various properties of a widely used experimental conﬁguration with two crystals and an air gap between them, namely, an SU(1,1) interferometer. In particular, it predicts the narrowing of the intensity distribution, the reduction and shift of the side lobes, and the decline in the interference visibility as the mean photon number increases due to stronger pumping. The presented experimental results conﬁrm the validity of the new approach. The model can be easily extended to the case of the frequency spectrum, frequency Schmidt modes, and other experimental conﬁgurations.


I. INTRODUCTION
At a high parametric gain, parametric down-conversion (PDC) and four-wave mixing (FWM) generate a bright squeezed vacuum (BSV).The BSV is a nonclassical state of light without a coherent component (displacement) but with a large (macroscopic) number of photons per mode.This state has strong photon-number correlations (twin-beam squeezing) [1][2][3][4], quadrature squeezing [5], multimode structure [6][7][8], and polarization entanglement if the generated photons have orthogonal polarizations [9].The BSV is a promising tool for a lot of applications in quantum optics and metrology: imaging [10][11][12][13][14], quantum state engineering [15,16], nonlinear interferometry [17][18][19], super-resolution, and phase sensitivity beyond the shot-noise limit [20,21].Due to the high Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.mean photon number and the large number of modes involved, a theoretical description of a BSV is complicated.
Several works on a description of a multimode BSV are based on the coupled differential equations for the signal and idler plane-wave operators under the plane-wave pump approximation [22][23][24][25].This approximation implies that each plane-wave signal mode interacts with only one plane-wave idler mode.This simplifies the problem significantly and provides analytical expressions for the output operators.In other works, equations similar to the classical propagation equations were derived and their solutions based on the Green function method were suggested [26].Integrodifferential equations for PDC with a fixed spectral profile of the pump were written in Refs.[27][28][29][30].Multimode PDC in the frequency domain was studied in Refs.[31,32] using the Magnus expansion.
The broadband-mode approach for the temporal domain based on the independent Schmidt modes was introduced in Ref. [27] and developed in Ref. [28]; the Schmidt mode approach for the spatial domain was developed and applied to experiment in Ref. [33].This approach describes several effects observed for a BSV and is very convenient for the analytical treatment of the problem.However, the existing Schmidt mode theory neglects the energy mismatch between the pump, signal, and idler photons and therefore leads to gain-independent shapes of the Schmidt modes.For this reason, it cannot describe the broadening of the spectrum, which is observed in experiment [34] as the BSV gets brighter due to the increase in the parametric gain (stronger pumping).
Moreover, the behavior of BSV properties with the increase of the parametric gain can be completely different depending on the geometry of experiment.For PDC in a single crystal, the spatial intensity distribution has a typical "sinc-squared" shape at low parametric gain, and it broadens as the parametric gain increases (see the results below).Similar behavior was observed for the frequency spectrum in Ref. [34].In contrast, a two-crystal configuration with an air gap in between leads to a complicated interference pattern of intensity with side lobes, which gets narrower with increasing parametric gain [17,33,35].This configuration is known as the SU(1,1) interferometer [18], and it has recently attracted a lot of attention due to its metrological applications [21,[36][37][38][39].
In this work, we present a new theoretical approach to the description of the spatial properties of the spatially multimode BSV taking into account the mode correlations.Our approach is based on exchanging the (k, t ) and (ω, z) representations and solving the high-dimensional system of integrodifferential equations for plane-wave operators.This approach describes various features of the BSV, such as the intensity distribution and the shapes of the Schmidt modes, as well as their evolution with increasing parametric gain, both in the case of a single crystal and in the case of a two-crystal configurations.In full agreement with the experiment, the theory predicts the broadening of both the intensity distribution and the Schmidt mode shapes with increasing gain in the case of a single crystal and the reduction of the side lobes in the two-crystal configuration.The suggested approach does not include any limitations on the pump waist width and the number of modes, as it was the case in the previous considerations, and it does not assume that the Hamiltonian commutes with itself at different moments of time.
The paper is organized as follows.Section II describes the theoretical approach and applies it to the case of highgain PDC in a single crystal.The Schmidt modes at variable parametric gain are considered in Sec.III.Section IV deals with the two-crystal configuration.In all sections, experimental results are also presented and compared with the theory.Finally, Sec.V is the conclusion.

II. HIGH-GAIN PDC IN A SINGLE CRYSTAL
The Hamiltonian of PDC in a crystal with a quadratic susceptibility χ (2) (r) is given by [23] where Ês,i are electromagnetic field operators for signal/idler photons, the pump is assumed to be a classical beam with a Gaussian envelope, propagating along the z axis, 2σ 2 e i(k p z−ω p t ) , with the full width at half maximum (FWHM) of the intensity distribution being 2 √ ln 2σ .By using the quantization of the electromagnetic field, ,where a † k s,i are the creation plane-wave operators, C k s,i are the coefficients of the decomposition, the Hamiltonian becomes Here we neglect the dependence of the coefficients C k s,i on k s,i and suppose that the interaction strength ˜ , involving χ (2) (r), the pump field amplitude, and other parameters is a constant.In the existing Schmidt mode approach, unlike the current approach, the frequency mismatch ω s + ω i − ω p is neglected, which makes the Hamiltonian independent of time and ultimately leads to gain-independent Schmidt modes.
For simplicity we consider a 2D model, using only one transverse coordinate x.In systems with the radial symmetry, without loss of generality, the 2D case can be easily extended to a 3D case by taking into account the additional integral over the y coordinate.After integration over x and substituting dk s,i = dq s,i dk sz,iz , where q s,i are the transverse wave vectors and k sz,iz are the longitudinal wave vectors of signal (idler) radiation, the Hamiltonian can be represented in the form ) This Hamiltonian is written in the momentum-time (k, t ) representation.In this picture, the Heisenberg equation of motion for the signal plane-wave operators takes the form × e i(k p −k sz −k iz )z e i(ω s +ω i −ω p )t a † i (q i , k iz , t ), (4) and similarly for the idler operators.The operators a † s,i (q s,i , k sz,iz , t ) and a s,i (q s,i , k sz,iz , t ) defined before are the slowly varying parts of creation and annihilation operators.In other words, these operators are solutions to the Heisenberg equation in the rotating frame of reference.Taking into account the free-propagation Hamiltonian for the signal and idler fields, one can obtain solutions in the fixed frame of reference.We called these solutions fast varying components.The fast varying components are connected with the slowly varying parts as a † s,i (q s,i , k sz,iz , t ) = e iω s,i t a † s,i (q s,i , k sz,iz , t ), a † s,i (q s,i , ξ, ωs,i ) = a † s,i (q s,i , ξ, ωs,i )e −ik sz,iz ( ωs,i ,q)ξ .(5) Here we assume long pump pulses and short crystals, so that the pulse covers the entire crystal at the same time.The Fourier transformation allows one to pass from (q, k z , t ) to (q, z, ω) representation.For example, for the fast varying part of the idler creation operator, the Fourier transformation is After substituting ( 6) and ( 5) into (4) and integrating its left and right parts from τ 0 = −∞ to τ = +∞ over the interaction time, Eq. ( 4) takes the form In real experiments, the final time τ corresponds to the end of a long interaction.Please note that k iz ( ωi , q i ) is a function of ωi and q i , while k iz is the integration variable.The operators at the final a s (q s , k sz , τ ) and initial a s (q s , k sz , τ 0 ) moments of time describe the boundary conditions.Due to these conditions, we have an equality between the operators in the end (beginning) of the interaction and operators corresponding to further (previous) free propagation in the linear medium.Outside the nonlinear medium, k and ω are connected by the dispersion law ω = ω(k) such that a s,i (q s,i , k sz,iz , t ) does not depend on t and a s,i (q s,i , z, ω s,i ) does not depend on z.In this case, Eq. ( 6) with the substitution k sz = ω s u ksz leads to the link between boundary conditions in different representations: a s (q s , k sz , τ ) = u ksz a s (q s , L, ω s ), a s (q s , k sz , τ 0 ) = u ksz a s (q s , 0, ω s ) [23], where u ksz is the projection of the group velocity vector of the signal photon on the z axis and L is the length of the nonlinear medium.
In the right-hand side of Eq. ( 7), integration over time and longitudinal momentum leads to the δ functions, the first of them defines the idler frequency through the signal and pump frequencies.
The δ functions allow one to take integrals over d ω and d z and simplify Eq. (7): where = ˜ /u kz and we assume u ksz = u kiz = u kz .In Eq (9), the idler frequency is defined through the signal and pump frequencies due to Eq. ( 8); to emphasize it, we wrote the frequency dependence explicitly for k iz in the brackets.Differentiation of the left-and right-sides of Eq. ( 9) with respect to the L leads to the coupled integrodifferential equations for the signal/idler annihilation/creation operators, where A similar equation can be written for the idler frequency, The solution to the system of coupled integrodifferential equations ( 10) and ( 11) can be found in the form a s (q s , L, ω s ) = a s (q s ) + dq s η(q s , q s , L, )a s (q s ) + dq i β(q s , q i , L, )a † i (q i ), + dq s β * (q i , q s , L, )a s (q s ), (12) where η(q s , q s , L, ), β(q s , q i , L, ) are functions of the transverse wave vectors, crystal length and the interaction strength, and a s (q s ) = a s (q s , L = 0, ω s ), a † i (q i ) = a † i (q i , L = 0, ω i ) are the initial plane-wave operators.In the case of an arbitrary pump and a high gain, the functions η(q s , q s , L, ) and β(q s , q i , L, ) can be found only numerically.However, in the low-gain regime, Eqs. ( 12) coincide with the solutions given by the Schmidt mode approach [33].In addition, with the spatially broad pump, Eqs. ( 12) coincide with the well-known analytical solution for the plane-wave pump approximation [23].
Finally, by solving the system of integrodifferential equations (10) and (11) in the form of Eqs.(12), various characteristics of the BSV can be found.For example, the mean photon-number distribution over transverse wave vectors is In what follows, we assume small angles of emission θ s,i , so that q s ≈ k s θ s , and consider angular intensity distributions.Please note that, to generalize our approach, we have considered signal and idler photons distinguishable in at least one degree of freedom (for example, noncollinear emission) and marked their operators with different indices s and i.Nevertheless, this approach can be simplified to the case of degenerate type-I PDC under the substitution a s = a i and with taking into account the corresponding commutation relation.
To find the connection between the theoretical parameter and the measured experimental gain, we have modelled the single plane-wave mode by calculating the total photon number in the collinear direction (q s = q i = 0) from Eq. ( 13) as a function of and fitted this dependence by the well-known dependence for the single-mode regime, y = B sinh 2 (A ), where A and B are the fitting parameters.Then the parametric gain is defined as G = A .A similar procedure was performed in the experiment.Using a pinhole, the dependence of the total intensity in the collinear direction on the parametric gain was measured and fitted by the function y = B e sinh 2 (A e √ P), where A e and B e are the fitting parameters, P is the pump power.In this case, the experimental gain is defined as G e = A e √ P, and the theoretical and experimental gains have to coincide.
To compare the predictions of the described theory with the experiment, we considered a 2-mm-thick BBO crystal and a pump laser with the wavelength 354.7 nm and with a beam waist of FWHM 170 μm.The angular intensity distributions of type-I PDC were calculated using Eq. ( 13) for different values of the parametric gain G, as shown in Fig. 1.For ideal experimental conditions the phase mismatch δk = k p (θ axis , ω p ) − 2k s,i (ω s,i ) = 0, with θ axis being the angle between the pump wave vector and the optic axis of the crystal.Due to imperfect crystal alignment, the phase mismatch can be slightly nonzero.Such a deviation is hard to fix in the experiment but has a considerable effect on the shape of the intensity distribution, which is shown in panels (a) and (b) of Fig. 1.
It is clearly seen that with increasing parametric gain the angular intensity distribution broadens.The broadening is directly connected with the fact that the angular variables and the parametric gain both enter the functions β and η.In contrast, in the Schmidt mode approach, the angular variables and the parametric gain are separate: the eigenmodes depend only on the angular variables, the eigenvalues depend only on the parametric gain, the broadening is not expected.The origin of the broadening can be clearly seen under the plane-wave pump approximation, where strong correlations between the signal and the idler photons take place: q s = −q i .FIG. 2. The theoretical FWHM of the BSV intensity distribution vs parametric gain for a 170 μm FWHM pump (cyan solid line) and plane-wave pump (blue dashed line) calculated according to the approach of Refs.[23,24], and the experimental data (red dots).
In experiment, a BSV was obtained through PDC pumped by the third harmonic radiation (wavelength 354.7 nm) of a pulsed Nd:YAG laser.The pulsed radiation (pulse width 18 ps, repetition rate 1 kHz) is essential to reach the highgain regime.The intensity distributions were recorded with a charge-coupled device (CCD) camera in the Fourier plane of a lens with the focal length of 100 mm.The spectral filtering was performed using a band-pass filter with the transmission centered around the wavelength 710 nm and with a bandwidth of 10 nm.
The dependence of the FWHM of the spatial intensity distribution of the BSV on the parametric gain for δk = 4530 m −1 (obtained by fitting the experimental distributions and fixed for further calculations) was calculated and compared with the experimental data, see Fig. 2. In the low-gain regime, the FWHM of the BSV intensity distribution coincides with the value calculated using the first-order perturbation theory [40].As the parametric gain increases, the FWHM monotonically grows.The same tendency is observed in the experiment (red points) and is in good agreement with the theoretical dependence (cyan line).The blue dashed line, calculated for the case of a plane-wave pump, predicts a slower broadening of the angular distribution than the one with a focused pump.Note that in the gain-independent Schmidt mode approach, the FWHM decreases with increasing gain, contrary to experiment.

III. SCHMIDT MODES
The BSV radiation is strongly multimode.This multimode structure is important for a lot of applications [11] but, at the same time, is difficult to analyze.The most useful way to describe the multimode BSV radiation is by introducing a system of normalized orthogonal Schmidt modes.Within the Schmidt mode basis, each signal mode is only correlated with a single matching idler mode, which greatly facilitates the analysis.
In a simplified Schmidt mode approach [33], the shapes of the BSV Schmidt modes do not depend on the parametric gain.In addition, the natural mode competition mechanism, i.e., low-order modes acquiring larger weights at larger gain, leads to the reduction in the number of Schmidt modes with increasing gain.These two statements lead to the angular (and frequency) narrowing of the intensity distribution, which is in contradiction with the observed broadening.The broadening can be only understood in the framework of the new approach considered here, and, as shown below, is connected with the mode widths changing with the gain.
The complex function β in Eq. ( 12) can be written using the Schmidt decomposition with respect to the transverse wave vectors, β(q s , q i , L, ) = n n e iφ n u n (q s , )ψ n (q i , ), (14) where gain-dependent eigenfunctions u n (q s , ) and ψ n (q i , ) for the signal and idler beams, respectively, are labeled with the index n, while n represent gain-dependent weights of the decomposition, and φ n are constant phases.Similarly, the function η can be decomposed with eigenvalues ˜ n , constant phases ϕ n , and eigenfunctions v n (q s , ) and ξ n (q s , ): From direct numerical decomposition of β(q s , q i , L, ) and η(q s , q s , L, ) the following relations between eigenfunctions are ensued: v n (q, ) = u n (q, ), ξ * n (q, ) = ψ n (q, ), and the absolute values of all functions are equal: 2 ), where λ n are eigenvalues of the Schmidt decomposition of the two-photon amplitude [33]) but are getting closer with increasing the gain and become equal in the high-gain regime ).This behavior of the renormalized weights is in full agreement with the gainindependent Schmidt mode approach [33].Using the Schmidt decompositions in Eqs. ( 14) and ( 15), we can introduce new photon creation/annihilation operators for the collective spatial Schmidt modes (the Schmidt mode operators) of the radiation as a result of the nonlinear interaction, The operators A † n and B † n have the same form but they are related with the signal and idler plane-wave creation operators, respectively.Equations ( 12) can be written in terms of the Schmidt operators, Equations ( 17) clearly show that the output signal and idler plane-wave operators are connected with the same functions u n (q, ) = v n (q, ), which is because we assumed frequency degeneracy.
The input/output relations (17) are similar to the ones of Ref. [33].However, in Eqs. ( 17) not only the weights n , ˜ n but also the functions u n (q s , ), v n (q i , ) depend on the parametric gain .Thereby, the output operators are now defined by the Schmidt modes whose shapes are gain-dependent.Using Eqs.(17), the intensity distribution Eq. ( 13) can be written in a simple form as a sum of the squared absolute values of the Schmidt modes with the corresponding weights: The Schmidt eigenmodes and eigenvalues of the BSV can be reconstructed from the covariance of its intensity distribution [41].Indeed, consider the sum of the contributions of the signal and idler radiation for a fixed gain, i.e., I (q) = I s (q) + I i (q), the covariance of intensities measured at positions q and q is defined as Cov(q, q ) = I (q)I (q ) − I (q) I (q ) . ( Calculation of the covariance distribution in terms of the Schmidt modes, using the input/output relations of Eqs. ( 17), leads to where we suppose that n = ˜ n for high gain.The first two terms of Eq. ( 20) are related to the autocorrelation of intensity fluctuations, respectively, of the signal and idler beams, while the third one represents the cross-correlation between the signal and idler radiation.
For simplifying the reconstruction of the Schmidt modes, in experiment we eliminated the cross-correlations between the signal and idler radiation by filtering a wavelength slightly shifted from the degeneracy point, so that the detected signal photons did not have idler matches.In this case, in the covariance (20) only the first term should remain, containing the signal Schmidt modes u n (q).Simultaneously, if the filtered wavelengths are still rather close to degeneracy, one can assume u n (q) = v n (q).The Schmidt modes and weights were found by performing the singular value decomposition (SVD) of the square root of the covariance distribution [42].
In experiment, we filtered the signal radiation using a band-pass filter with the central wavelength 700 nm and a bandwidth of 10 nm.Since the length of the nonlinear crystal was small enough (2 mm), we could neglect the effect of spatial walk-off.Accordingly, the BSV radiation was axially symmetric and one could assume the factorability of the x and y degrees of freedom.Therefore the analysis was restricted to the intensity profiles along the x direction.For a better signal-to-noise ratio, we integrated the intensity distributions in the y direction within the range from −2 to 2 mrad of the angle q y /|q|, where q = (q x , q y ).Around 2000 single-shot intensity profiles I (θ ) with θ = q x /|q| were measured, and the 2D covariance distribution Cov(θ, θ ) was calculated.
Figure 3 shows the experimental distribution of the covariance for G e = 6 in normalized units.The autocorrelation part is distributed along the main diagonal, where θ = θ .Conversely, the cross-correlation, which would lead to nonzero covariance values along the complementary diagonal θ = −θ , is absent due to spectral filtering.In order to get rid of the background noise, a 2D fit was performed on the covariance distribution.Given the high-gain version of the intensity profile of the PDC radiation [24] and the fact that the covariance distribution along the main diagonal behaves as the squared intensity profile, the function used for the fit was with A, B, C, D being fitting parameters and E being an experimentally determined quantity dependent on the phase mismatch.The inset of Fig. 3 demonstrates the fitted distribution for G e = 6, which is indeed equivalent to the experimental one.The white dashed line in the main figure represents the half-maximum level of the fitted covariance distribution and shows a good agreement with the half-maximum level (in cyan) for the experiment.The broadening of the covariance distribution with the increase of the parametric gain is shown in Fig. 4. The experimental FWHM follows the predicted trend from the theory for both the main and the complementary diagonals.For low gain, the theoretical dependence coincides with the covariance obtained through the first-order perturbation theory.The dependence of the main diagonal on the gain (blue solid curve (Please note the different axis scales.)The solid lines represent theoretical calculations, while the points stand for the experimental data.The dashed blue line corresponds to the main diagonal of covariance calculated under the plane-wave pump approximation [23,24]. in Fig. 4) has a minimum.This minimum is also observed in the case of a plane-wave pump (blue dashed line in Fig. 4) and qualitatively separates the low-and high-gain regimes.
The shapes and the weights of the Schmidt modes can be obtained through the SVD of the function √ Cov(θ, θ ).This procedure has been performed for both the theoretical, Eq. (20), and the fitted experimental covariance distributions.The results of the reconstruction for different gain values are shown in Fig. 5 for the first and the second Schmidt modes.The shapes of the modes are close to the Hermite functions and their widths depend on the gain.The general tendency is the broadening of the Schmidt modes with increasing gain, this broadening being more pronounced for the higher-order modes.The theoretical results show a good agreement with the experimental data.

IV. TWO-CRYSTAL CONFIGURATION
The method described above can be extended to the twocrystal configuration, with the two crystals separated by an air gap of length d, known in the literature as the SU(1,1) interferometer [7,18,[36][37][38].The earlier models, for example, the gain-independent Schmidt mode approach, describe sufficiently well the BSV at the output of an SU(1,1) interferometer, including the redistribution of the mode weights and the narrowing of the spatial distribution.Qualitatively an agreement between theory and experiment is observed.Quantitatively, there is a small disagreement in the width of the spatial distribution, especially as the gain increases [33].This disagreement is eliminated in the new model.
In the two-crystal case, one should take into account that during the free propagation in the air gap the signal, idler and pump photons acquire an additional phase.This phase creates a factor of e i k d standing in front of the integral over the second crystal, where k = k air p − k air s − k air i is the wavevector mismatch in the air [43].Considering the problem step by step, we modify Eq. ( 7) by taking into account the δ function conditions in Eq. (8): The boundary conditions connect the beginning and the end of the interaction in different representations: a s (q s , k sz , t = 0) ∼ a s (q s , L = 0, ω s ), a s (q s , k sz , t = τ ) ∼ a s (q s , 2L, ω s ).Differentiation of Eq. ( 22) in L gives da s (q s , 2L, 23) and for the idler creation operator, where a † i (q i , L, ω p − ω s ) and a s (q s , L, ω s ) can be found by solving Eqs.(10) and (11) for the single crystal.
The intensity distribution in the presence of the air gap is completely different from the single-crystal case and depends on the length of the air gap.Due to the different refractive indices of the pump and BSV photons in the air, an additional phase is acquired in the gap and the intensity of light emitted in the collinear direction oscillates from minimum to maximum as d increases.Also, an increase in d leads to more and more frequent interference fringes in the intensity distribution [Fig.6(a)].
In the experiment, we pumped two nonlinear crystals with the third harmonic of a pulsed Nd:YAG laser (repetition rate 50 Hz, wavelength 354.7 nm, and pulse duration 29.4 ps) with a FWHM diameter of approximately 0.3 mm.The two BBO crystals (3-mm thick, cut for degenerate type-I PDC) were aligned, in turn, for degenerate phase matching.A dichroic mirror and a color-glass filter suppressed the pump after the nonlinear interaction.A band-pass filter selected a 10-nm bandwidth (FWHM) of the PDC around the wavelength 710 nm.A lens brought the PDC to the momentum space, where a CMOS camera was introduced.On the camera, the background was subtracted and the data was acquired for 200 ms.The pump energy per pulse was measured before the crystals with a calibrated energy meter.The distance between the crystals was varied by changing the position of the first crystal using a translation stage.Neutral density filters were used to avoid the saturation of the camera.The ring patterns measured with the camera were then transformed into polar coordinates.To increase the accuracy, the radial profiles were obtained by averaging out the polar plots over the azimuthal angle.Due to the radial symmetry, this procedure is equivalent to fixing one of the Cartesian coordinates in theoretical calculations and calculating the intensity distribution over the other coordinate.
As the pump power increases, the parametric gain in each crystal grows (the gain after two crystals with a gap grows nonmonotonically).Figures 6(a)-6(c) shows the resulting spectra both calculated (blue) and measured (red), with the parametric gain and the pump pulse energy shown in each panel.(In the low-gain case, the measurement was not possible because of the small intensity.)Apart of a small shift in the fringes, increasing the pump power leads to the reduction of the side peaks.Therefore the envelope of the angular distribution gets narrower as the parametric gain grows, in contrast to the single-crystal case.Moreover, from Figs. 6(a)-6(c), it is clearly seen that the visibility of the interference fringes drops down with increasing pump power.Note that with stronger pumping the Kerr effect, which leads to an additional phase (mostly manifested in the collinear direction [41], becomes more pronounced and was taken into account here. In Figs.6(d)-6(f), one can observe that with increasing the distance between the crystals, the total width of the angular distribution for the same pump power (or for the same gain in each crystal) is reduced, as reported in Ref. [33].This happens due to diffraction, which leads to the reduction of the angular width of the BSV that overlaps with the pump and is amplified in the second crystal; diffraction is more pronounced for larger distances.As it was mentioned above, the second mechanism leading to the narrowing of the spatial intensity distribution in the two-crystal configuration is the increase of the pump power (the parametric gain in each crystal).Simultaneously, these two mechanisms diminish the number of the Schmidt modes in the system, allowing one to create different shapes of intensity distributions and to control the number of modes in the BSV.

V. CONCLUSION
We have presented a new theoretical approach to describe the spatial properties of a BSV generated through high-gain PDC.In this approach, we derived and solved the integrodifferential equations for plane-wave operators without limitations on the pump waist width, number of modes and commutation of the Hamiltonian at different moments of time.The developed approach correctly captures a lot of features of the BSV.On the one hand, it is compatible with the Schmidt mode representation.On the other hand, it properly describes the broadening of the angular distribution with increasing parametric gain.As a result, the new treatment correctly predicts the dependence of the Schmidt mode widths on the parametric gain.The model describes different experimental configurations: the singlecrystal case and the configuration of two crystals with an air gap between them.For the verification of our theoretical model, we have performed several experiments, both with a single-crystal and with a two-crystal PDC sources [SU(1,1) interferometer].The presented experimental results are in good agreement with performed theoretical calculations.Our model gives a deep insight into the properties of highgain PDC, its mode structure and the origin of nonclassical correlations.
The gain-dependent weights n and ˜ n are different in the low-gain regime ( n = sinh 2 ( √ λ n ) and ˜ n = 4 sinh 4 ( √ λ n

FIG. 3 .
FIG. 3. Experimental distribution of the covariance in normal-ized units for the parametric gain G e = 6.Cross-correlations are removed by using a filter selecting only the signal radiation.For comparison, the white dashed line denotes the 0.5 level of the fitted distribution, shown in the inset.

FIG. 4 .
FIG. 4. Dependence of the FWHM of the covariance main (blue) and complementary (orange) diagonals on the parametric gain.(Pleasenote the different axis scales.)The solid lines represent theoretical calculations, while the points stand for the experimental data.The dashed blue line corresponds to the main diagonal of covariance calculated under the plane-wave pump approximation[23,24].

FIG. 5 .
FIG. 5. (a) The first and (b) the second Schmidt modes for different parametric gain values: G = 5.0 (red), 5.7 (blue), and 6.6 (black).Solid lines represent theoretical calculations.Dashed lines stand for modes retrieved from the experiment.

FIG. 6 .
FIG. 6.The normalized intensity distributions of the BSV in the two-crystal configuration with an air gap.[(a)-(c)] The distance between the crystals is 10.66 mm, the gain in each crystal is (a) G = 0.01, (b) 1.1, and (c) 2.45.[(d)-(f)] The gain is fixed, G = 3.15, the distance between the crystals is (d) 5.58, (e) 10.66, and (f) 23.36 mm.Blue curves are theoretical calculations, red curves present the experimental data.The legends also show the values of the pump energy per pulse.