Waves in spatio-spectral and -temporal coherence of evolving ultra-intense twin beams

Waves in the spatio-spectral and -temporal coherence of evolving ultra-intense twin beams are predicted: Twin beams with low intensities attain maximal coherence in the beam center until certain threshold intensity is reached. Then the area of maximal coherence moves with increasing intensity from the beam center towards its edges leaving the beam center with low coherence (the first coherence wave). For even larger intensities, a new coherence maximum is gradually built in the beam center with the increasing intensity and, later, it again moves towards the beam edges forming the second coherence wave. Rotationally-symmetric twin beams are analyzed within a three-dimensional model that couples spectral and spatial degrees of freedom. Relation between the twin-beam coherence and its local density of modes during the nonlinear evolution is discussed.

Parametric down-conversion, in parallel with its applications in classical nonlinear optics 1 , nowadays serves as a common source of three types of highly nonclassical states 2 : states of individual entangled photon pairs, intense twin beams composed of macroscopic numbers of photon pairs and states with squeezed phase and photon-number fluctuations (squeezed, sub-Poissonian light). The states of individual photon pairs have been intensively studied for already more than 40 years 3 . Several milestones in modern physics were reached with their help including the experimental demonstration of the Bell-inequalities violation 4 and experimental teleportation of the polarization state of a photon 5 . Also the squeezed states that are provided by the degenerate variant of parametric down-conversion in resonators have been investigated extensively 6,7 . This has resulted in the construction of sources of squeezed light that exhibit the reduction of vacuum fluctuations better than 15 dB 8 . Both these types of states are at present, apart from their applications in ultra-precise metrology, considered as promising for quantum-information processing, especially in connection with optical communications. On the other hand, much lower attention has been paid to intense twin beams (TWB) that emerge in parametric down-conversion in the running-wave configuration and contain macroscopic, i.e. large, numbers of photon pairs [9][10][11][12][13][14][15][16][17][18][19] . Their sub-shot-noise intensity correlations between the signal and idler beams that constitute a TWB have been demonstrated [20][21][22][23][24][25] and even exploited for sub-shot-noise imaging 26 . Also tight spatial correlations between the signal and idler beams that originate in phase matching of the nonlinear interaction have been analyzed in different configurations [27][28][29][30][31] . Applications of TWBs aiming at improving the experimental precision have been suggested in spectroscopy 32 and quantum interferometry 33 . Also TWBs with squeezed spectrally broad-band phase fluctuations have been suggested 34 .
The generated TWBs are naturally highly multi-mode, both in their spectra and spatial transverse modes 9,30,[35][36][37] . The need to generate TWBs composed of ideally only one spatio-spectral mode triggered investigations of both spectral and spatial coherence of intense TWBs 29,30,36 . An increase of both spectral and spatial coherence with an increasing TWB intensity has been experimentally confirmed for the beams with lower intensities 30 . On the other hand, it has been observed that, for the sufficiently intense TWBs, degradation of coherence occurs 36 . This degradation originates in considerable depletion of the pump beam and we speak about ultra-intense twin beams in this case. This behavior has been confirmed by direct numerical solution of the nonlinear Maxwell equations 36 . A model of intense parametric down-conversion based upon the decomposition of optical beams into the Schmidt dual modes shed more light into the dynamics of formation of such ultra-intense TWBs by considering the evolution of individual modes inside the interacting pump, signal and idler beams 38,39 . It revealed that an increase (decrease) of coherence of a TWB results from the reduction (increase) of the number of effectively populated modes inside the TWB, which occurs as a consequence of the dynamics of individual modes' triplets each being Figure 1. A signal (or idler) beam in its transverse wave-vector plane at four subsequent stages during its evolution in a sufficiently-long nonlinear crystal: (a) broad weak beam with only weakly developed coherence in the central part, (b) narrow stronger beam with well developed coherence in the central part, (c) broad strong beam with well developed coherence at the edges, and (d) narrow even stronger beam with well developed coherence in the central part. The degree of coherence is quantitatively indicated by the widths of red ellipses. The first coherence wave in the TWB evolution occurs between stages (b) and (c). The second coherence wave begins to propagate from stage (d). Qualitatively the same evolution occurs also in the plane in which the radial transverse wave-vector k is substituted by frequency ω.
we not only provide this generalization we also incorporate into the model the mutual coupling between the spectral and radial wave-vector degrees of freedom in a rotationally-symmetric TWB. This allows us to analyze also spectrally non-degenerate TWBs.
At the end of Introduction, we would like to note that the theory developed here for parametric down-conversion 1 can be modified such that nonlinear resonant interactions involving four-wave mixing in cold atomic ensembles [42][43][44] are described. Compared to parametric down-conversion, much higher effective nonlinear coupling constants occur for resonant interactions and also parameters of the generated optical beams are more appropriate when subsequent absorptions by atoms are considered. High values of sub-shot-noise intensity correlations between the beams were observed in such systems [45][46][47] . Combinations of several nonlinear resonant processes were also used to generate more complex multi-mode states endowed with quantum correlations [48][49][50][51] . Whereas the properties of intense TWBs arising in parametric interactions are more appropriate for metrology applications, intense TWBs and their generalizations reached in (multiple) resonant interactions in atomic ensembles are suitable for implementations of various quantum-information protocols. Metrology applications whose measurement uncertainties are suppressed below the classical limit 20,34 due to the TWB sub-shot-noise intensity correlations may further profit from additional stability due to the large number of photon pairs inside intense TWBs 32,33 .
The paper is organized as follows. The suggested model is described in section Three-dimensional model of an intense twin beam. Quantities that characterize an intense TWB including its coherence are determined in section Physical quantities characterizing twin beams. TWB intensities, numbers of modes and intensity profiles are discussed in section Twin-beam intensities, intensity profiles and numbers of internal modes. TWB coherence as it evolves with the increasing pump power and typical profiles of intensity auto-and cross-correlation functions are analyzed in section Spatial, spectral and temporal coherence of twin beams, coherence waves. Section Twin-beam coherence and density of modes is devoted to the relation between the TWB coherence and its density of modes. Evolution of the coherence of a TWB propagating along a long crystal is described in section Coherence of a twin beam propagating along the crystal. Section Conclusions summarizes the main findings.

three-Dimensional Model of an Intense twin Beam
Parametric down-conversion as the second-order nonlinear process is characterized by tensor χ (2) of second-order nonlinear susceptibility that mediates an effective nonlinear interaction among the signal (index s), idler (i) and pump (p) beams. Assuming propagation of these beams along the z axis an interaction momentum operator G int appropriate for this process is written as follows 52,53 :ˆ∫ and r = (x, y, z). The positive-frequency part of a classical pump electric-field amplitude is denoted as + E p ( ) similarly as the negative-frequency part of a signal [idler] electric-field operator amplitude is referred as ]. Symbol H.c. replaces the Hermitian conjugated term, ε 0 means the permittivity of vacuum, and symbol: shorthands the χ (2) tensor with respect to its three indices.
We decompose the electric-field amplitudes ] of the interacting beams into harmonic plane waves with wave vectors k a (ω a ) and frequencies ω a : c is the speed of light in vacuum. We invoke the paraxial approximation for the propagating beams and so we identify the harmonic plane waves by the transverse parts ⊥ k a of their wave vectors and frequencies ω a . An intense pump beam with amplitude ξ p , central frequency p 0 ω , propagating along the z axis and having Gaussian profiles both in the transverse plane (beam radius w p ) and temporal spectrum (pulse duration τ p ) is assumed to impinge on a nonlinear crystal of length L: for the incident pump-beam amplitude. In Eq. (3), we assume that the pump-beam transverse profile does not change along the z axis if the nonlinear interaction is neglected. This 'cylindrical approximation' means that the crystal length has to be smaller than the Rayleigh range of the used pump Gaussian beam.
The complex nonlinear dynamics of three interacting beams can successfully be treated when we introduce suitable spatio-spectral modes in each beam that approximately ' diagonalize' the interaction among three beams and that replace the usual harmonic plane waves 39,40 . Such spatio-spectral modes' triplets can be found by first considering the Schmidt dual modes in the signal and idler beam revealed by the Schmidt decomposition 54 of a two-photon spatio-spectral amplitude Φ si of a state containing one photon pair (for details, see below). A pump-beam mode that corresponds to a given pair of signal and idler modes and belongs to a common modes' triplet is obtained by considering one nonlinear quantum event with simultaneous annihilation of a signal and an idler photon and creation of a pump photon in a classical fashion 40 : A signal photon with its classical electric-field www.nature.com/scientificreports www.nature.com/scientificreports/ amplitude in the form of a given mode function together with its idler twin with the classical electric-field amplitude given by its mode function undergo classical nonlinear interaction that provides, according to the rules of classical nonlinear optics, the profile of classical electric-field amplitude of a pump photon. This physically motivated approach has two drawbacks. First, the obtained pump-mode profiles are not normalized to unity, as needed in any quantum model. Second, the pump-mode profiles determined for different Schmidt dual signal and idler modes are not exactly orthogonal. The first drawback is easily removed by re-normalizing the pump modes to unity 40 and simultaneously changing the corresponding nonlinear effective coupling constant for the corresponding modes' triplet. As for the second drawback, numerical analysis of non-orthogonality of pump modes reveals that the largest overlap is between the nearest neighbors and drops down fast as the distance between modes (expressed by the difference of modes index numbers) grows. According to our experience with down-conversion pumped by pulses with duration in units of ps, the greatest overlaps of pump mode functions of the nearest neighbors are lower than 15%, the overlaps for the first couple of more distant modes are only in units of %, and the remaining overlaps are practically zero. This behavior arises from the typical oscillatory behavior of pump-mode profiles, a qth 1D mode has 2q + 1 peaks in its intensity profile for q = 0, 1, …. We note that a qth 1D signal (or idler) mode has only q + 1 peaks in its intensity profile. Omission of these overlaps is the essence of the used approximation. It leads to the picture of many independent modes' triplets with their own internal dynamics. As there exists no parallel quantum model of intense TWBs, validity and applicability of the used model can only be judged by comparing its predictions with the experiment and classical numerical simulations that mimic quantum behavior. The model was successfully used to interpret the properties of experimental TWBs in ref. 38 . On the other hand, it was shown in ref. 36 that the properties of experimental TWBs accord with the results of classical numerical simulations based on the solution of the nonlinear Maxwell equations. In our opinion, highly multi-mode character of intense TWBs importantly contributes to the applicability of this approximation.
Contrary to the case of spectrally degenerate TWBs analyzed in refs. 39,40 , the consideration of spectrally non-degenerate TWBs requires a general procedure for obtaining the dual Schmidt modes. The procedure begins with the determination of the state |ψ si 〉 of one photon pair leaving the nonlinear crystal. For the crystal that extends from z = −L to z = 0 the state |ψ si 〉 is obtained as a first-order perturbation solution of the Schrödinger equation assuming the incident vacuum state |vac si 〉 in the signal and idler beams: a a a creates a photon in mode ω ⊥ k ( , ) a a and  stands for the reduced Planck constant. The explicit formula for the two-photon amplitude Φ si introduced in Eq. (4) can be found in ref. 55 .
In the model, we assume rotational symmetry (with respect to the z axis along which the pump beam propagates) of the analyzed system. This means that the model is suitable for type-I parametric down-conversion in collinear or not-to-far-from collinear geometry for sufficiently wide pump-beam transverse profiles, as discussed in detail in ref. 12 . This symmetry implies that the general two-photon amplitude Φ si depends only on the difference ϕ s − ϕ i of azimuthal angles of the signal and idler modes [ The azimuthal functions Φ m occurring in the azimuthal decomposition (5) are determined along the formula: To arrive at a Schmidt-mode decomposition 56,57 of the azimuthal two-photon amplitudes Φ m that admits factorization in its variables as a suitable approximation, we rotate the orthogonal coordinate systems ω  Fig. 2(b). The corresponding transformations are expressed as follows: The azimuthal two-photon amplitudes Φ m written in the transformed variables can then be approximately replaced by their factorized forms Φ ∼ m : The subsequent Schmidt decompositions of functions , then allow us to arrive at the complete decomposition of the general two-photon amplitude Φ ∼ expressed in the transformed variables in the following 'factorized' form: In Eq. (11), the global Schmidt numbers λ mlq are given as www.nature.com/scientificreports www.nature.com/scientificreports/ mlq m ml mq the azimuthal Schmidt number λ m is determined along the formula Using the Schmidt functions  u a ml , and f a mq ,  occurring in Eqs. (9) and (10), respectively, and those defined in the azimuthal angles ϕ a , m m s, s s i, i i the overall Schmidt functions g a mlq ,  introduced in Eq. (11) are expressed by the relation: a a a a defined for the harmonic plane-wave modes by the fields' operators a a mlq , of the overall Schmidt modes g a,mlq (expressed in the original variables ⊥ k a and ω a ), , a = s, i, given in Eq. (2) 55 : Considering an impinging classical strong pump beam that depletes during its propagation in the nonlinear medium and the signal-and idler-beam electric-field operator-amplitude decompositions written in Eq. (17), we can approximately replace the original momentum operator G int given in Eq. (1) by the following momentum operator Ĝ int av 55 : The overall coupling constant ∼ K is given as where the overall pump-beam amplitude is determined for photon numbers. A pump-beam amplitude A p,mlq occurring in Eq. (18) characterizes the pump mode associated with an (mlq)th signal and idler dual Schmidt mode. We assume that the overall number of pump photons present in the incident pump beam is distributed into these modes linearly proportionally to the squared overall Schmidt coefficients mlq 2 λ , i.e. the incident pump-beam amplitude A mlq p,  of mode mlq related to normal ordering of field's operators is given as This approach follows from the perturbation solution of the Schrödinger equation and the position of the Schmidt coefficients in the decomposition of the two-photon amplitude Φ si . Approximative ortho-normality of the pump mode functions plays an important role in this approach. As discussed above the obtained pump modes are not normalized to unity and their re-normalization refines the model 40 . The corresponding re-normalization of the coupling constant ∼ K is obtained as follows: A polychromatic signal photon with averaged energy ω s 0 and mode function ω ϕ . Such pump-photon mode function is easily determined in the crystal output plane (polar coordinates r a and ϕ a , a = p, s, i) and time domain: In Eq. (20), the signal-and idler-beam mode functions g r t ( , , ) a mlq rt a a a , ϕ are given for the rotationally symmetric beams according to the formula we express the suggested re-normalization by the following substitutions: guarantees the correct number of pump photons, as it can be verified by the perturbation solution.
If the rotation angle α [or more precisely the angle α′ that transforms the variables δ ⊥ k a and δω a /c in agreement with the transformation in Eq. (7)], that expresses the strength of coupling between the spatial and spectral variables, is small and we neglect the part of re-normalization associated with the beams transverse plane (for sufficiently wide beams), we arrive at the simplified formula  (18) for parametric down-conversion with the incident vacuum signal and idler beams was found in ref. 39 in the generalized parametric approximation. In this approximation, independent evolution of the interacting modes' triplets is first treated classically (using symmetric ordering of fields' operators) which provides the pump-mode amplitude A p,mlq of an (mlq)th triplet along the crystal in the form:  Then, in the second step, the linear Heisenberg equations for the signal and idler fields' operators that incorporate the classical pump-mode solution, , are solved. The solution is obtained in the following general form: The overall number K n of modes describing the TWB can be determined, e.g., by the formula derived for a multi-mode thermal field 55,58 : As an interesting alternative, the number K of modes can also be obtained directly from the statistics of photon pairs, as suggested in ref. 55 :    symbol : :  normally orders fields' operators and x x x Δ = − 〈 〉ˆˆ denotes the fluctuation of operator x. We note that, in parallel to the correlation functions in Eqs (38) and (39), the spectrally and spatially averaged intensity correlation functions can be defined and more easily determined (for details, see ref. 39 ). Substituting the relations (16) and (29) into Eqs (38) and (39) we arrive at the formulas for intensity auto-and cross-correlation functions:  Temporal properties of a TWB observed in the far field belong to important experimental characteristics. They are quantified by the appropriate intensities and intensity auto-and cross-correlation functions that are determined by the formulas analogous to those written above. Instead of the mode functions g k ( , , ) a mlq kt a a a a ml a a a mq a a a a , , ,

Twin-beam intensities, intensity profiles and numbers of internal modes
To reveal the behavior as well as coherence of TWBs with the increasing pump power we consider a typical experimental configuration based on non-collinear spectrally non-degenerate type-I down-conversion in a BBO crystal of length L = 4 mm. The crystal is pumped by a pump beam of radius w p = 500 μm, central wavelength i 0 θ = . deg. We note that similar experimental configurations tuned for spectrally degenerate non-collinear down-conversion were used in refs 36,38 where the experimental generation of ultra-intense TWBs in the regime with pump depletion was investigated by measuring the spatially and spectrally averaged intensity correlation functions. According to these experimental investigations, ultra-intense TWBs are expected for pump powers in tens or hundreds of mW in similar configurations when usual nonlinear crystals are used. For nonlinear structures with higher effective χ (2) nonlinearities, ultra-intense TWBs are reached for lower pump powers inversely proportional to the squared χ (2) coefficient. We also note that the previous theoretical investigations of intense TWBs 39,40,55 were performed for spectrally degenerate TWBs as a consequence of limitations of the used model. Contrary to this, we consider here spectrally non-degenerate TWBs to demonstrate the range of applicability of the developed model. We also provide the comparison with the previously analyzed degenerate cases. Below, we analyze different TWBs occurring at the output plane of a crystal with ä fixed length that are generated by pump beams with the increasing intensities. This corresponds to the way how the experimental investigation of TWBs is performed. As shown below, a TWB intensity increases with the increasing pump-beam intensity. The properties of TWBs consecutively obtained this way are similar to those reached when we allow a pump beam with given initial intensity to propagate through a sufficiently long crystal. We note that the solution of the nonlinear interaction for one modes' triplet depends on the crystal length L via the quantity P L involving the pump power P 39 . Properties of such propagating TWB including its coherence have been discussed in Introduction and schematically depicted in Fig. 1.

twin-beam intensities and numbers of internal modes. We first analyze the behavior of 'integral'
parameters of TWBs. This gives us a general framework for detailed discussion of the behavior of correlation functions on one side. On the other side, the comparison of properties of the analyzed spectrally non-degenerate TWBs with those being spectrally degenerate 39 allows us to elucidate the role of the difference of the (central) signal and idler frequencies in forming a TWB during the nonlinear dynamics.
For our spectrally non-degenerate TWBs, the overall number N of photon pairs increases when the pump power P increases. Whereas this increase is exponential for lower pump powers (due to the properties of parametric down-conversion with an un-depleted pump beam), it gradually becomes linear for sufficiently high pump powers. As evident from the curve in Fig. 3(a) that gives the number N s of signal photons as a function of pump power P, both areas are well separated. Explanation of this behavior was given in ref. 39 , in relation to the number K of effectively populated TWB modes: During the initial TWB amplification from the signal-and idler-beam vacuum state, the spatio-spectral modes with greater nonlinear coupling constants ≡ ∼ K KA mlq mlq ps are amplified faster than the modes with lower coupling constants K mlq and this behavior gradually suppresses the modes with lower coupling constants in the TWB. As a consequence, the number K of TWB modes decreases for lower pump powers P, until certain threshold pump power P th is reached [see Fig. 3(b)]. For the pump powers P greater than the threshold power P th , the number K of modes again increases. The reason is that the TWB modes with the greatest coupling constants K mlq are already at the stage of evolution in which they return back their energy to the www.nature.com/scientificreports www.nature.com/scientificreports/ corresponding pump mode. This diminishes the role of such modes in the structure of the TWB and it allows the modes with lower coupling constants 'to be seen' again in the TWB. For the analyzed TWB, this threshold occurs around P th = 250 mW and, as the curves for the numbers K and K n of modes plotted in Fig. 3(b) show, there occurs also the second threshold around the power P th,2 ≈ 8 W 39 . This second threshold is observed at the pump powers for which the TWB modes with the greatest coupling constants, that loose their energy after the first threshold power P th and later take back this energy, start again to loose their energy. We can see in Fig. 3(b) that the reduction of the number K of modes is smaller for the powers around P th,2 compared to that around P th and, moreover, there occur two groups of populated modes with different values of their coupling constants in this case (for details, see ref. 39 ). The presence of two groups of modes around P th,2 results in weakening of the effects of narrowing the TWB profiles and increasing the TWB coherence compared to those observed in the area around P th , as discussed below.
The reduction of the number K of TWB modes naturally leads to narrowing of TWB intensity profiles both in radial wave-vector and frequency domains [see Fig. 3(c) for the widths ΔI s,k of signal-beam radial wave-vector intensity profiles]. This narrowing has its origin in the properties of the Schmidt modes defined in the (k ⊥ , ω) plane: The smaller the coupling constant of a modes' triplet the wider the mode functions in the (k ⊥ , ω) plane (for typical intensity mode profiles, see ref. 59 ). It holds both in the spectral and temporal domains, as well as in the transverse wave-vector domain, that the intensity profiles of a qth mode exhibit q + 1 local peaks and their overall widths increase with q. The reduction of the number of TWB modes thus causes simultaneous narrowing of both spectral and temporal intensity profiles for the pump powers around P th , as documented in Fig. 3(d). We note that this behavior, that stems from the nonlinearity of the analyzed interaction, qualitatively contrasts with the property of the (linear) Fourier transform that assigns longer pulses to narrower spectra.
The described behavior is qualitatively similar to that found for spectrally degenerate TWBs 39 . Thus, the spectral non-degeneracy, that qualitatively influences the properties of weak fields composed of individual photon pairs, plays only a minor role in forming the basic properties of ultra-intense TWBs. This is related to the fact that the dynamics of a TWB effectively evolves in greater number of orthogonal Schmidt spectral modes, each being spread over a large part of the spectrum.
Twin-beam intensity profiles. Now, we pay attention to the evolution of spatial, spectral and temporal TWB intensity profiles towards ultra-intense beams. To qualitatively discuss this evolution, we compare in Fig. 4 the shapes of spectrally resolved intensity profiles ω ⊥ I k ( ) k s, s in the radial wave-vector direction of the signal beam for three specific cases. The shapes of a weak signal beam, as more-less formed by spontaneous parametric down-conversion, are  www.nature.com/scientificreports www.nature.com/scientificreports/ plotted in Fig. 4(a). In the frequency domain and our configuration, the signal-and idler-beam spectra match each other: They together form a common spectrum of down-conversion observed in a given propagation direction. The shapes of intensity profiles ω ⊥ I k ( ) k s, s drawn in Fig. 4(b) for a TWB around the threshold power P th exhibit considerable narrowing in the whole (k ⊥ , ω) plane. These profiles broaden in the (k ⊥ , ω) plane for the pump powers P exceeding P th , as documented in Fig. 4(c) for a TWB generated for P = 1.5 W. Moreover and importantly, the spectral profiles attain a two-peak structure, which is a consequence of the inverted flow of energy in the modes with the greatest coupling constants. We note that the spectral modes with even mode numbers (q = 0, 2, …) have central peaks in their intensity profiles and greater coupling constants compared to their closest lower neighbors with odd mode numbers and central dips in their profiles. On the other hand, the intensity profiles ω ⊥ I k ( ) k s, s in the radial wave-vector direction remain single peak.
The comparison of shapes of temporal intensity profiles I t ( ) t k s, s resolved in the radial wave vector and plotted for the signal beam in Fig. 5 reveals qualitative similarity with the already discussed spectral shapes of Fig. 4. Also here, there occurs narrowing of the profiles around the threshold power P th and two-peak structure of the temporal profiles for powers P exceeding P th . This two-peak structure is a consequence of the properties of temporal mode functions that are qualitatively similar to those of the spectral mode functions discussed above.

spatial, spectral and temporal coherence of twin beams, coherence waves
Now we come to the central results of our article discussing the behavior of spatial, spectral and temporal intensity auto-and cross-correlation functions in different areas of the TWB profiles as they evolve with the increasing pump power.
Local coherence of twin beams as it depends on pump power, coherence waves. The suppression of modes with larger mode numbers m, l and q (and thus smaller coupling constants) with the increasing pump power P, that narrows the intensity profiles, broadens the intensity auto-and cross-correlation functions. This can be understood by the fact that the complex mode functions belonging to the modes with small mode numbers m, l and q change slowly their phases along their intensity profiles. As follows from the comparison of curves in Fig. 6(a-d) showing the widths Δ ω A k s, , Δ ω A k s, , Δ ϕ A k s, , and A t k s, Δ of the signal-beam intensity auto-correlation functions in in turn frequency, radial wave-vector, azimuthal wave-vector, and time domain, this occurs more-less synchronously in all three directions of the signal-beam 'phase-space' spanned by transverse wave vector and frequency, or time instead of frequency. According to the curves in Fig. 6(a-d), in the centers of intensity profiles the maxima of the widths of intensity auto-correlation functions are found around the threshold power P th . This accords with the results of ref. 39 where the spatially and spectrally averaged intensity correlation functions were analyzed. However, the curves drawn in Fig. 6(a-d) reveal that these maxima gradually move towards the greater pump powers P when we move away from the centers of intensity profiles towards their edges. Even at the edges of the intensity profiles, the reached maximal widths of the intensity correlation functions are comparable to those observed at the centers of intensity profiles, but at greater pump powers. Moreover, these maxima occur for TWBs effectively composed of much larger numbers K of modes compared to that determined at the threshold power P th . This means that the improved coherence at the TWB edges is caused by some form of local phase synchronization of the Schmidt modes. This behavior forms coherence waves in the TWB 'phase-space' , as discussed in Introduction.
To understand systematically the behavior of local threshold powers P th in the whole signal-beam 'phase-space' , we determine their values from the behavior of intensity auto-correlation functions A computed along two orthogonal directions in the plane k ( , ) s s ω ⊥ that coincide with the axes δ ⊥ k s r and δω s r of the rotated coordinate system [see Fig. 2(c)]. In each point in these directions, the local threshold powers P th can be inferred from the correlation functions defined in three orthogonal directions, i.e. when we consider them as functions of radial wave vector (k s ⊥ ), frequency (ω s ) and azimuthal wave-vector angle (ϕ s ). The curves plotted in Fig. 7(a,b) for the cuts along the k s r δ ⊥ and δω s r axes show that the local threshold powers P th belonging to the correlation functions www.nature.com/scientificreports www.nature.com/scientificreports/ in three orthogonal directions are close to each other. This is true also for the local threshold power P th derived from the temporal intensity auto-correlation functions A t k s, at the δ ⊥ k s r axis, as evidenced in the curves of Fig. 7(b). We can roughly say that the smaller the local intensity I s in a given point in the signal-beam 'phase-space' the greater the local threshold power P th (in all orthogonal directions). This suggests an idea that certain local Schmidt modes can be defined around an analyzed point in the TWB 'phase-space' and the evolution of local coherence can be explained from the point of view of the evolution of populations of these local modes in the same vein as we did in the case of the whole TWB. This idea generalizes the approach that estimates the number of modes as the Fedorov ratio 60 , i.e. from the spread of local intensity correlation functions.
Typical profiles of intensity auto-and cross-correlation functions. Before we conclude this section, we discuss typical profiles of the intensity auto-and cross-correlation functions that have been analyzed above. It holds that the averaged intensity auto-correlation functions A in the transverse wave-vector and frequency domains are broader than their cross-correlation counterparts C for low pump powers P and they tend to coincide for greater powers P 39 . This behavior is observed also for the investigated not-averaged intensity auto-and cross-correlation functions in the areas of TWB 'phase-space' where the local phase-matching conditions are close to the ideal ones found in the central part of the TWB. We note that the ideal phase-matching conditions occur in the 'phase-space' at positions ω ϕ ϕ = and parameter ϕ varies from 0 to 2π. Under these conditions, the intensity correlation functions form either a single smooth peak or small oscillations may be observed at the shoulders of this peak, as documented in Fig. 8(a,b) for the normalized radial wave-vector intensity auto-( ω A k s, ,n ) and cross-( ω C k s, ,n ) correlation functions. However, if there occurs certain local phase mismatch in the investigated 'phase-space' point, the relation between the auto-and cross-correlation functions qualitatively changes. The cross-correlation functions, that are strongly sensitive to the local phase mismatch, shift their maxima into a different than expected point in the 'phase-space' where better phase-matching conditions are found. In this case, the cross-correlation functions C are considerably broader than the corresponding auto-correlation functions A. The normalized azimuthal intensity cross-correlation  www.nature.com/scientificreports www.nature.com/scientificreports/ ,n together with their much narrower azimuthal auto-correlation counterparts ϕ A k s, ,n plotted in Fig. 8(c,d) for the points outside the TWB central part may serve as an example. In this case, some of the cross-correlation intensity profiles C k s, ,n ϕ have a two-peak structure.

twin-beam coherence and density of modes
In this section, we analyze the coherence of three typical TWBs with different intensities and elucidate the relation between their coherence and densities of modes.
twin-beam coherence. We consider a weak TWB (P = 1 × 10 −8 W), a TWB generated around the threshold power (P = 0.5 W) and a TWB for which P = 1.5 W and mutually compare their properties. Similarly as above, we analyze the TWB coherence along the axes δ ⊥ k s r and s r δω of the rotated coordinate system. As examples, the widths A k s, Δ ω and Δ ω A k s, of signal-beam intensity auto-correlation functions given in radial wave-vector and frequency domains, respectively, are drawn in Fig. 9(a,b) for three considered TWBs.
Whereas the signal-beam coherence is more-less uniform in the signal-beam 'phase-space' for a weak TWB, it considerably varies in the 'phase-space' of more intense TWBs. The analyzed more-intense TWBs generated with the pump powers 0.5 W and 1.5 W represent in certain sense two boundary cases. Whereas the maximal coherence in the less-intense TWB is observed in the center of its signal-and idler-beam profiles, the maximal coherence is reached at the edges of the profiles of the more-intense TWB.
When we consider the TWBs with powers P increasing from the power P = 0.5 W to the power P = 1.5 W their area of maximal coherence moves from the signal-and idler-beam centers towards their edges. This is documented in Fig. 10 where the widths Δ ω A k s, , Δ ϕ A k s, and A k s, Δ ω of signal-beam intensity auto-correlation functions determined in turn in frequency, azimuthal wave-vector and radial wave-vector domains are plotted along the axes k s r δ ⊥ and s r δω as functions of pump power P. For the TWBs generated by even greater pump powers P > 1.5 W, the area of maximal coherence being at the edges of the TWB with P = 1.5 W moves further away from the beam center and disappears in the tails of beam profiles for sufficiently large powers P (see Fig. 10). As follows from the graphs in Fig. 10, qualitatively the same behavior of TWB coherence as observed for pump powers around the  www.nature.com/scientificreports www.nature.com/scientificreports/ first threshold power P th is found for the pump powers around the second threshold power P th,2 . Also here coherence maxima occur in the beam center for the pump powers P below the threshold power P th,2 and they gradually move away from the beam center as the pump power P increases for P > P th,2 . We can speak about coherence waves describing the propagation of coherence maxima across the beam profiles. The first two coherence waves are observed for the pump powers P > P th and P > P th,2 . The coherence waves propagate in the plane spanned by the signal (as well as idler) transverse wave-vector and frequency axes. The propagation of these waves reveals a detailed structure (origin) of coherence oscillations found along the power axis P for the averaged intensity correlation functions in ref. 39 .
Densities of twin-beam modes. The change of coherence across a TWB profile can be interpreted as the change of the local densities K a of modes, a = s, i, defined as the ratio of the volume in 'phase-space' occupied by beam a and the local coherence volume: ω ⊥ and the normalization constant K should be set such that the density K a integrated over the whole beam profile gives the number  of modes determined by Eq. (35). According to the formula (44) the better the local beam coherence the smaller the local beam density K a of modes. Thus, in agreement with the observation of the previous subsection, close-to-uniform densities of modes characterize weak TWBs. Intense nonlinear interaction then introduces changes into the densities of modes first reducing them in the central part of a TWB. For pump powers obeying P > P th , increase of the TWB intensity causes gradual movement of the area with the reduced densities of modes towards the edges of the TWB profile and then this area disappears in the TWB profile tails. Similar scenario is observed for greater pump powers P comparable or larger than the second threshold power P th,2 . This means that there exist waves in the densities of modes that are synchronized with those in the TWB coherence. To illustrate this behavior we plot in Fig. 11(a,b) the normalized signal-beam densities K s n of modes for the cuts along the axes δ ⊥ k s r and s r δω of the rotated coordinate system. We remark that the normalized densities K k s ,n of modes considered as a function of radial transverse wave vector ⊥ k s in Fig. 11(b) show systematic increase of their values with the increasing k s ⊥ . This reflects the rotational symmetry of the TWB that requires greater numbers of azimuthal modes for larger values of ⊥ k s .

Coherence of a twin beam propagating along the crystal
Before we conclude, let us rephrase the above findings obtained for TWBs generated from a crystal of fixed length by varying the pump power in terms of the development of the coherence of a TWB during its propagation in a sufficiently-long crystal. Considering, e.g., the signal beam in its 'phase-space' spanned by the transverse wave vector and frequency (plotted along the propagation z axis), its shape resembles a 'blurred' torus (rotationally symmetric around the z axis) with certain radius at low pump powers. This torus is gradually deflating as the signal beam propagates through the crystal, attains its minimal volume, then it is inflating until certain position in the crystal is reached. At this position it starts to deflate again. Starting from the input crystal face the TWB coherence increases inside the whole torus, the area of maximal coherence lies in the center of the torus. At the position where the volume of the torus is minimal (corresponding to the first threshold power P th ), the area of maximal coherence starts to move towards the surface of the torus. It crosses the surface of torus in the area inside the crystal where the volume of torus is already large, continues its propagation outside the torus and gradually disappears. This represents the first coherence wave. Later in the crystal, similar behavior occurs. The area with coherence maxima is again formed in the center of the torus, but at certain position in the crystal (corresponding to the second threshold power P th,2 ) it starts to move towards the surface of torus forming the second coherence www.nature.com/scientificreports www.nature.com/scientificreports/ wave. Specific stages in the signal-beam evolution are depicted in Fig. 1 considering two 2D cuts of the torus by hyper-planes k ( , ) s s ϕ ⊥ and (ω s , ϕ s ).

Conclusions
We have predicted spatial, spectral and temporal waves in the coherence of ultra-intense twin beams during their evolution in the nonlinear interaction. Our prediction is based upon the analysis of a three-dimensional model of an intense rotationally-symmetric twin beam that was developed in the transverse wave-vector and frequency domains and then extended to the time domain. The model is based on the Schmidt decomposition of two-photon spatio-spectral amplitude. It treats the nonlinear evolution of parametric down-conversion in terms of suitable approximately independent modes' triplets in the regime with pump-beam depletion. Intensity profiles of twin beams were determined as they change with the increasing pump power in a crystal of fixed length together with the intensity auto-and cross-correlation functions that vary across the twin-beam profiles. It was revealed that the maximal coherence of a twin beam is observed in the twin-beam central part provided that the generating pump power is lower than the threshold power. For greater pump powers above the threshold power, the rotationally-symmetric area of maximal coherence moves from the center of the signal (or idler) beam profile towards its edges as the pump power grows. These two-dimensional coherence waves across the signal (as well as the idler) beam profile drawn in the radial wave-vector and frequency plane are accompanied by two-dimensional waves in the local density of twin-beam modes. We believe that the predicted waves in the spatial, spectral, and temporal coherence of ultra-intense twin beams and the provided explanation will stimulate further theoretical as well as experimental investigations of the coherence of intense optical beams and open the door for applications of such beams, e.g., in metrology.

Data Availability
All data generated and analyzed during this study are included in this published article.