{'fraction_non_alphanumeric': 0.03840884285813808, 'fraction_numerical': 0.00938160369690918, 'mean_word_length': 4.39041181624734, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mode structure and photon number correlations in squeezed quantum pulses

The question of efficient multimode description of optical pulses is studied. We show that a relatively very small number of nonmonochromatic modes can be sufficient for a complete quantum description of pulses with Gaussian quadrature statistics. For example, a three-mode description was enough to reproduce the experimental data of photon number correlations in optical solitons [S. Spalter et al., Phys. Rev. Lett. 81, 786 (1998)]. This approach is very useful for a detailed understanding of squeezing properties of soliton pulses with the main potential for quantum communication with continuous variables. We show how homodyne detection and/or measurements of photon number correlations can be used to determine the quantum state of the multi-mode field. We also discuss a possible way of physical separation of the nonmonochromatic modes.


I. INTRODUCTION
For a complete quantum description of an optical pulse, one has to use a multimode density matrix. The question is, how many modes are necessary for such a description. If one works with monochromatic modes, then an infinite number of modes would be needed, which is impractical. On the other hand, one can construct different sets of modes as linear combinations of the monochromatic modes [1]. These modes (which are however nonmonochromatic) may be more suitable for the quantum description of pulses. In particular, it would be useful to find such a modal structure so that the quantum state of the pulse can be described by the density operator̺ = ̺ exc ⊗̺ vac , where̺ exc is some nontrivial density operator of a few modes and̺ vac is the vacuum state density operator of all the remaining modes. Thus, we could work with a concise description of the pulse by means of̺ exc .
Our work is motivated by the question how to completely describe the quantum state of solitons in optical fibers. Such soliton pulses are used in various schemes of quantum information processing (for a review see, e.g., [2]). Can a single-mode description of a soliton pulse be sufficient? When explaining squeezing in the Kerr medium, one sometimes plots a picture of phase space where a quantum uncertainty circle is deformed into an ellipse-a description which clearly corresponds to a single mode situation. This approach is useful to qualitatively understand the generation of continuous variables entanglement as in [3], but do we use all the quantum features of a given pulse if we treat it as a single mode object? Should one, instead, work with many more modes of a pulse in the hope that each of them can become a useful resource for quantum communication? Multimode correlations of photon numbers in a soliton were observed [4]. To numerically calculate such correlations, quantum variables of a soliton were treated on a position grid of typically ∼ 10 2 − 10 3 points [5,6]. Are there hundreds of useful modes available in a pulse, or would a better choice of the mode functions show that only a few modes (perhaps just one?) are in some nontrivial quantum state? Or are just four quantum operators introduced by Haus and Lai [7] sufficient to describe all the relevant phenomena?
Knowing the answer to the question of what is the complete quantum description of a pulse would be very helpful in quantum information processing: one could fully utilize the squeezing and entanglement properties of our sources. After forming a soliton pulse in a fiber, one can determine its multimode quantum state in the most comprehensive way. One can then optimize the medium properties to achieve maximum squeezing, or maximum purity of an entangled state. Working with pairs of correlated pulses, one could apply a proper measurement scheme and use entanglement criteria for multimode bipartite Gaussian states [9] to check wheteher the pulse pair is entangled or separable. One can also better understand the influence of the medium on the propagating pulse: provided that the output pulse is deformed, what happens with the quantum information carried by the pulse? Is it washed out by decoherence processes (or perhaps by an eavesdropper), or is it just unitarily transformed into other modes of the same pulse? A simple and correct measurement and description of the multimode state is highly desirable.
This work is organized as follows. In Sec. II we introduce nonmonochromatic modes and deal with the transformations between different sets of quadrature operators. Sec. III studies basic properties of the photon statistics of multimode fields: mean photon numbers and covariances between different modes. In Sec. IV we discuss a homodyne scheme for a complete determination of Gaussian states of nonmonochromatic multi-mode fields. In Sec. V we compare the photon number squeezing available via applying a proper local oscillator modulation and via using spectral filtering. Sec. VI suggests a way for an optimal selection of the nonmonochromatic mode functions. In Sec. VII we discuss a possible way how to physically separate individual nonmonochromatic modes of the pulse. Discussion and conclusion are presented in Sec. VIII. Most of the mathematical details are discussed in Appendixes A-E, and in App. F we compare our approach to that of Haus and Lai [7].

II. NONMONOCHROMATIC MODES
Under the term "mode" we understand a single degree of freedom of the electromagnetic field; a mode can be described by a pair of bosonic operators. It can be monochromatic (evolution described by a single frequency) or nonmonochromatic. Any state of the field can be treated in different mode decompositions. In a given decomposition the state is called single mode if all modes except of one have vacuum statistics of their operators, in the opposite case the state is multimode. In this section we deal with the transition between different mode decompositions.

A. Bosonic operators
Let us assume polarized optical field propagating in a given direction. The monochromatic modes of the field are denoted by the corresponding frequency ω, and the bosonic operators of these modes satisfy the commutation relations Let us assume a complete orthonormal set of functions We define a new set of operatorsb k aŝ which satisfy the commutation relations Thus, the functions f k (ω) can be used to define a new set of nonmonochromatic modes. The inverse transformation of the nonmonochromatic modes to the monochromatic ones readsâ

B. Example
Let us consider a pulse with a normalized frequency envelope f (ω); |f (ω)| 2 dω = 1. Let us define an orthonormal system of mode functions as in Eqs. (2) and (3) with f 1 (ω) ≡ f (ω). Let the quantum state of the first modeb 1 be a coherent state |β 1 ,b 1 |β 1 = β|β 1 , and let all the other b-modes be in vacuum |0 k ,b k |0 k = 0, k > 1. By means of the transformation (6) we find that using the monochromatic modes, the quantum state is a multi-mode coherent state ω |f (ω)β ω . Even though the single-mode and multimode coherent states are related to the same object, i.e., a pulse in a coherent state, the single-mode description is clearly more convenient for handling and for understanding the field structure. Of course, for other than coherent states the transformations of the state in different mode-systems are not so straightforward, but still can provide us with a very convenient way of the quantum state description.

C. Quadrature operators
We define the hermitian quadrature operatorsx(ω), p(ω) asx and similarly the b-mode quadrature operatorsX k ,P k aŝ These operators obey the same commutation relations as the quantum mechanical position and momentum operators withh = 1, i.e., and The transformation between the two sets of quadrature operators can be written in the form with ξ = x, p, and µ = X, P , and the transformation matrix Z is As can be checked, matrix Z satisfies the orthogonality relations Sometimes it is useful to work with a discrete set of frequencies ω j corresponding to frequency bins of width ∆ω. Quadratureξ j of the jth bin is obtained aŝ where the transformation matrix elements Z µξ kj are For brevity, we can join the quadraturesx(ω) andp(ω) into a single vectorξ, and similarly the quadraturesX k andP k into a single vectorμ, so that the transformations (15) and (16) [or (21)] are written in the matrix multiplication form Here the superscript T means matrix transposition, and the rows of the matrix Z are represented by indices µ and k, and the columns by the indices ξ and ω (or j).

D. Mean quadratures and variances
Important properties of multimode quantum states are given by the vectors of the quadrature mean values, ξ j ≡ Tr[̺ξ j ] and µ k ≡ Tr[̺μ k ], and by the variance matrices V and V ′ , given as where {. . .} stands for the anticommutator, {Â,B} ≡ AB +BÂ, and . . . represents averaging Â ≡ Tr(̺Â). The relationship between these quantities can be written in the matrix multiplication form as In general, to transform between µ, V , and ξ, V ′ , one has to know the complete transformation matrix Z, i.e., the whole set of mode functions f k (ω). However, if we assume that only a few nonmonochromatic modes are excited while the rest is in vacuum state, then the explicit form of the empty mode functions does not play any role. (It was also illustrated in the Example II B where only one mode function f (ω) was enough to transform the state into the multi-mode case.) This can be used for a simplified calculation of the transformed quantities, as discussed in App. A.

III. PHOTON STATISTICS OF MULTIMODE GAUSSIAN STATES
Among all possible states, the class of Gaussian states is one of the most important-both from the theoretical and experimental point of view. The basic property of these states is that their Wigner function (as well as the Q-function or the characteristic function) have Gaussian form (see, e.g., [8,9,10]). Examples of Gaussian states are coherent states, squeezed coherent states, thermal states, and squeezed thermal states. Gaussian states are typically observed in most experiments with optical pulses. In this paper we confine ourselves to Gaussian states; the main advantage used here is that a Gaussian state remains Gaussian in any mode decomposition and a relatively small number of parameters is necessary for its description.

A. Mean photon numbers and variances
Let us first study the relation between quadrature moments and photon number moments in discrete modes. The photon number operator in kth moden k can be expressed by means of the quadrature operators aŝ The mean photon number in kth mode is thus and the mean product of photon numbers is which can be written in terms of the quadrature means and variances as cov (n k , n l ) =x kxl V ′ xk,xl +x kpl V ′ xk,pl +p kxl V ′ pk,xl +p kpl V ′ pk,pl with the special case ∆n 2 k ≡ cov(n k , n k ). A very important role is played by the so called normally ordered covariance C kl ≡ : ∆n k ∆n l : = cov(n k , n l ) − δ kl n k .
The normally ordered covariance is used to define the normalized correlation matrix as in [4], : ∆n k ∆n l : Values of the correlation matrix measured in [4] are shown in Fig. 1a. In the limit of continuous frequency modes, one can define photon number density n(ω) , variance density ∆n 2 (ω) = n(ω) , and the normally ordered covariance C(ω, ω ′ ) with its normalized version C (n) (ω, ω ′ ). Since a Gaussian state remains Gaussian in any mode decomposition, and since the parameters of a Gaussian stateξ and V can easily be transformed from one mode decomposition into another, one can also calculate the photon number mean values and correlations in arbitrary mode decomposition.

B. Photon correlations in multimode pulses
If we try to reproduce the experimental results of [4] as truly as possible with as few modes as possible, the first attempt would be to start with just a single nonmonochromatic mode. However, as shown in Appendix C, the covariance matrix cov(n k , n l ) of such a field would have the same sign for all values k and l. The sign would be positive if the single-mode state is super-Poissonian, or negative if the single-mode state is sub-Poissonian. On the other hand, the covariance matrix measured in [4] contains both positive and negative elements. Therefore, the quantum state of the soliton pulses measured in [4] cannot be described as a single nonmonochromatic mode.
Exact analytical expressions are rather lengthy if more than one mode are excited. However, relatively simple relations can be obtained if the following requirements are met: (i) The coherent amplitude of at least one of the modes is much greater than any of the variance matrix elements, (ii) only the X-quadratures have non-zero mean values, i.e.,P k = 0 for all k, and (iii) the mode functions are real. While the condition (i) is generally valid for all strong pulses, the conditions (ii) and (iii) are related to our choice of mode functions and their generalization is straightforward. In Eq. (34) the terms containingp k,l kl as a function of the wavelengths of the pulse spectrum. Results measured in [4] with omitted elements with very large fluctuations (> 0.8) (a), and the best fit using a reconstructed three-mode variance matrix V Xk,Xk ′ (see Sec. III D) (b).
vanish and the first term on the right is dominant; thus Eq. (34) becomes To calculate this quantity using the nonmonochromatic modes, we writē (17), and (22)), and for V ′ xk,xl we use Eq. (A7) from Appendix A. The mean photon number is using Eqs. (31), (B6), and the conditions (i)-(iii). The normally ordered covariance (35) thus becomes Since in very narrow frequency bins ∆ω, the mean photon number and photon number variance are approximately equal, ∆n 2 k ≈ n k ≪ 1, Eq. (39) can be used also to obtain ∆n 2 k , and the normalized correlation matrix (36) becomes Taking ∆ω → 0, we can work with the continuous quantity In Eq. (42) only the X-elements of the quadrature variances V occur. In other words, the measured photon number covariances are only influenced by the covariances of quadratures which are in phase with the quadrature of the strongly excited modes.

C. Reconstruction of the quadrature variances from the photon number covariances
In the preceding subsection we have seen how to calculate the normalized photon covariance from the quadrature variance matrix of the nonmonochromatic modes. We can also consider an inverse problem. Let us assume that the mode functions f k (ω) are known. Let us also assume that the normally ordered covariances C(ω, ω ′ ) are measured as in [4] so that the normalized covariance C (n) (ω, ω ′ ) can be determined. As can be seen from Eq. (42), this quantity is a linear combination of the products of the mode functions f k (ω). The coefficients in the linear combination are elements of the quadrature variance V . Thus, to obtain the elements of V , one can use the orthonormality of the mode functions f k (ω) to invert the linear dependence. We find Thus, from the spectral covariances of the photon numbers one can reconstruct the N (N + 1)/2 elements of the quadrature variance matrix, out of the total number of N (2N + 1) independent elements. Since Eq. (43) is linear, it allows for a direct estimation of the reconstruction error (see, e.g., [10]). If the normalized covariances were measured with precision ∆C (n) (ω, ω ′ ), then the error in reconstructing the elements of V can be estimated as Note that here, for the sake of brevity, we have assumed uncorrelated errors in the elements of C (n) (ω, ω ′ ); derivation of a more general formula taking into account correlated errors is straightforward. Let us stress that the formulas derived in here and in Sec. III B are valid for very narrow frequency bins where ∆n 2 k ≈ n k ≪ 1. However, in real experiments, it is often necessary to work with wider bins in which the photon statistics can differ from Poissonian so that using Eqs. (42) and (43) could cause significant errors. This situation can be easily taken into account by substituting the proper expression for ∆n 2 k in the formulas connecting the photon correlation matrix with the quadrature variance matrix. The equations (42) and (43) then become slightly more involved; since this generalization is straightforward, we do not include the formulas in this text.

D. Application to experimental data
To illustrate our method, we have used the experimental data obtained in [4] (reproduced in Fig. 1a). They were measured for soliton pulses propagating in a 2.7 m optical fiber. We have used a preliminary set of mode functions (starting with a sech function as the basic shape of the soliton) and applied the reconstruction formula (43). We have diagonalized the resulting variance matrix V Xk,Xk ′ and found that only three eigenvalues are substantially different from the vacuum value 1/2. The corresponding eigenvectors can be used to construct a new set of mode functions (see Fig. 2). In this set, the X-quadratures are independent of each other. Thus, we have found that with respect to the quadratures which are in phase with the coherent amplitude (X-quadratures in our convention), the pulse measured in [4] can effectively be described as a three mode field, the quadrature variances being V X1,X1 ≈ 0.29, V X2,X2 ≈ 1.39, and V X3,X3 ≈ 2.69 with the off-diagonal elements being zero. The smallest eigenvalue corresponds to squeezing −2.35 dB. The reconstructed variance matrix V Xk,Xk ′ can be used to calculate back the photon number correlation matrix C To estimate the precision of our results, we have Monte-Carlo generated 1000 "experimental" matrices C (n) kl with elements fluctuating with errors given by the error estimates of the original experiment. The squeezing value then fluctuated between −2.2 dB and −4.1 dB, with most results centered around −3.1 dB. Even though the error is too large to make a definite statement, these result suggest that larger squeezing is available than the measured ≈ −2.5 dB of this setup with spectral filtering. Let us note that the shape of the mode function corresponding to the squeezed quadrature (full line in Fig. 2) resembles the spectral filtering approach: the contribution from the middle of the spectrum is enhanced while the outer parts are suppressed. The question of photon number squeezing availability via spectral filtering and via local oscillator functions is studied in more detail in Sec. V.

E. Conclusion
We have seen that one mode description of a pulse is not sufficient to explain the experimental data of [4], whereas a three mode description gives a good agreement with the measurements. However, we cannot conclude from this that the pulse does not contain more than three modes. In other modes the P -quadratures can be excited which do not influence the observed photon statistics. Also, the measurement noise was too high so that weak excitation of some additional modes might be undistinguished from the data noise backgraund.
It is interesting to compare our approach with that of Haus and Lai [7] (see also [11,12,13]). In their case the quantum field is decomposed into a "soliton" part and a "continuum" part. The soliton field is described by four operators ∆n, ∆θ, ∆x, and ∆p, related to the soliton energy, phase, position, and velocity, respectively. As shown in App. F, these operators can be expressed by our quadrature operators of four modes, provided that the mode functions are properly selected. However, only two of these operators (∆n and ∆x) related to two our modes (f 1 and f 2 ) influence the photon statistics. Thus, confining ourselves to the four soliton operators of [7] would not be sufficient to describe the observed phenomena. To apply the formalism of [7], one would have to work with the full set of the soliton and continuum operators.

IV. COMPLETE DETERMINATION OF THE VARIANCE MATRIX
Even though one can obtain full information about the X-quadratures from the spectral correlations of photon numbers, one has no access to the variances V P k,P k ′ and V Xk,P k ′ . To get full information on the multimode Gaussian quantum state, one has to perform phase dependent measurements. Homodyne detection is one example of such a measurement; recently it was studied how to apply homodyne detection scheme to the quantum state reconstruction of a multimode optical field [14,15]. Because we confine ourselves to the Gaussian states, the task is easier than reconstruction of a general quantum state.
To find the variance matrix V , one can use the scheme of [15] with the local oscillator pulses shaped to the form of weighted combinations of the mode functions. It is very useful to subtract the coherent amplitude of the pulse by using a balanced Sagnac interferometer [16] as in Fig. 3. Here, two counterpropagating identical squeezed pulses are interfering at a 50%/50% beam splitter. In one of the outputs of the beam splitter the pulses interfere constructively and form a bright pulse, which is then used to form the local oscillator. In the other output where the pulses interfere destructively, squeezed vacuum (or a more general field with zero mean amplitude) is formed. The squeezed vacuum pulse has the same quadrature variance matrix V as each of the two counterpropagating pulses.
The bright pulse is shaped interferometrically so that its envelope has the form of different combinations of the mode functions f k (ω). Such a pulse is then used as a local oscillator in a balanced homodyne detector. Thus, if the local oscillator is f k (ω), one can obtain the value of V Xk,Xk , whereas with the local oscillator if k (ω) one can obtain V P k,P k . Knowing these values and using the local oscillator pulse of the form f k (ω) + f k ′ (ω) one can obtain the value of V Xk,Xk ′ , by using the local oscillator form i(f k (ω) + f k ′ (ω)) one can obtain the value of V P k,P k ′ , and by using the form f k (ω) + if k ′ (ω) one can obtain the value of V Xm,P m ′ . Altogether, N (2N + 1) different forms of the local oscillator are sufficient to obtain all the N (2N + 1) independent elements of the variance matrix V . Let us note that to increase the precision of the measurements, one can increase the number of different phases of the local oscillator; a Maximum-Likelihood method for parameter estimations of a single mode Gaussian state using many phases has been discussed in [17].
Having found the variance matrix, one has all the information about the Gaussian state. Then we can immediately see, e.g., what is the maximum available squeezing of the state: it is the value corresponding to the minimum eigenvalue of V . Typically, this value will be smaller than the minimum diagonal element of V , which means that the optimum squeezing is shared among different modes. For example, there is a quadrature squeezing in the fundamental mode f 1 due to the Kerr effect (assuming mode functions as in App. F), but the X 1 quadrature is also correlated to the X 3 quadrature, because of the correlation between the pulse width and energy. Therefore, better squeezing can be expected to occur in a combination of the modes 1 and 3 than in the isolated mode 1. Thus, if the scheme is used to produce squeezing, the measurement scheme can serve as a tool for a selection of the optimum local oscillator.
One can also see how much are individual modes entangled with each other. To take full advantage of this knowledge one has to be able to separate individual modes from each other and distribute them among Alice, Bob, and other entanglement consumers. This is, however, a rather nontrivial task; in Sec. VII we will briefly mention a possible approach to its solution. Let us assume that we want to prepare a pulse with the maximum photon number squeezing, i.e., the photon number fluctuates as little as possible. To quantify the photon number squeezing, one uses the Mandel Q-parameter [18] defined as where n refers to the photon number of the entire pulse. This quantity is negative for sub-Poissonian field and positive for super-Poissonian fields. It was suggested in several works [19] that spectral filtering of the pulse can lead to improvement of the photon number squeezing by blocking frequency bands with correlated photon fluctuations and letting through frequency bands with anticorrelated photon numbers. Here we show that a proper modulation of the coherent amplitude (playing the role of the local oscillator) leads to the optimum squeezing which cannot be overcome by the spectral filtering method.

A. Local oscillator modulation
Let us first assume that the quadrature variance matrix V is fixed while we can modulate the mean values of the quadraturesX n andP n ; this corresponds to the experimental situation as in Sec. IV and Fig. 3. The photon number variance is Assuming as in Sec. III B that the coherent amplitudes are much greater than the variance matrix elements, and that the mode functions are real (generalization to complex functions being straightforward), we can express the photon number variance as +2X mPn V Xm,P n +P mPn V P m,P n .
The mean photon number is so that the Q-parameter of Eq. (45) can be written, after some algebra, as Thus the Q-parameter can be expressed in the matrix multiplication form as whereX andP are column vectors with theX m andP m elements, and V XX , V XP , and V P P are the corresponding submatrices of the variance matrix V . The variance matrix is multiplied from the left and from the right with a unit vector, so that the Q-parameter is limited by where V (min) is the minimum eigenvalue of the variance matrix V and V (max) is its maximum eigenvalue. The minimum value of Q is reached if the vector of mean values of X and P is the eigenvector of V corresponding to its minimum eigenvalue. Setting the components of this vector is possible using the scheme of Fig. 3.

B. Spectral filtering
Let us assume a spectral filtering function 0 ≤ c(ω) ≤

1: the frequency component is completely transmitted if c(ω) = 1 and it is completely blocked if c(ω) = 0. This function transforms the mean quadratures and variances asx
The δ kl terms in Eqs. (56) and (57) follow from the quantum mechanical nature of the quadratures: partially blocking a frequency component means that the corresponding field is mixed with vacuum. Eqs. (54)-(58) can be used to determine ∆n 2 and n as in the preceding subsection, and thus to find the Q-parameter of the filtered field. After a straightforward algebra, one can express the new Q-parameter as with Again, the Q-parameter is calculated as a matrix product of the variance V multiplied from the left and from the right by a vector. This time, however, the vector is not of the unit length, so that Q f is limited by where A 2 is the square of the magnitude of the multiplying vector, To show that Q f can never be smaller than the minimum value achievable by local oscillator modulation, it is enough to show that A 2 ≤ 1, i.e., that the magnitude of the multiplying vector is not bigger than one. Proof of this inequality is shown in Appendix D.

C. Conclusion
We can see that no spectral filtering can improve the photon number squeezing below the minimum eigenvalue of the quadrature variance matrix, which is available via the local oscillator modulation approach. Let us stress that our method of finding the optimum local oscillator function is very simple and straightforward since it is based on linear algebra. It was suggested recently to use an adaptive algorithm to optimize the pulse shape for achieving optimum photon number squeezing [20]. In this method, the optimum was reached after 20 000 iterations. In our case, provided that the pulse is sufficiently well described by N = 5 modes, N (2N + 1) = 55 iterations is enough.

VI. SELECTION OF APPROPRIATE MODE FUNCTIONS
So far it was assumed that the set of mode functions was given and the question was about the statistics of the corresponding quadratures. But how should these mode functions be selected?
In principle, any orthogonal set of mode functions can be used for description of the pulse statistics. However, since the aim is to reduce the number of quantum variables necessary to describe the pulse, the functions should be carefully chosen. One possibility is to start from theory and assume some particular shape of the pulse -e.g., a hyperbolic secant soliton and to construct the orthogonal set from the typical perturbations. An example is given in App. F, Eqs. (F15)-(F18) where relationship to the soliton quantum fluctuations approach by Haus and Lai [7] is discussed.
Another approach does not assume any particular mode function form and is related directly to the experiment. It is useful to assume that only one mode has a nonzero coherent amplitude and a very small number of modes have quadrature fluctuations substantially different from the vacuum values. To find the optimum set of mode functions, one can use the following procedure.
1. Select f 1 (ω) as the measured classical envelope of the pulse (see App. E).
2. Determine the minimum value of variance to be still considered as different from vacuum. 6. Diagonalize V (2) to getṼ (2) = W V (2) W T with W an orthogonal matrix. Let us choose W such thatṼ (2) 11 is the largest element ofṼ (2) ;Ṽ

Construct a temporary set of orthogonal functions f
k,l . 7. The function f 2 is constructed as It can be checked that the variance of the X quadrature corresponding to this mode function is 11 .

8.
A new temporary set of N ′ − 2 mode functions is constructed from the old one as an orthogonal complement to f 1 and f 2 . By means of the transformation connecting the new temporary set of mode functions to the initial one, calculate the corresponding variance matrix. Construct reduced matrix V (3) by excluding rows and columns referring to quadratures of modes 1 and 2. Diagonalize V (3) and find mode function f 3 in the same way as in 6 and 7. Note that the variances of the quadratures of mode 3 are smaller than the variance X 2 2 .
9. Repeat this procedure of redefining the temporary mode set, transforming and diagonalizing the variance matrix, and defining a new mode function. After each repetition, the maximum variance of the (n + 1)st mode is smaller than the maximum variance of the nth mode. If for some n = N the variance of mode n+1 is sufficiently close to 1/2 (as defined in item 2), the quantum state of the (n + 1)st mode is indistinguishable from vacuum. The pulse can be described with the given precision as an Nmode object.
Let us note that the procedure of redefining the temporary mode set and rediagonalizing the reduced variance matrix each time when a new mode function is constructed is necessary. It is not possible, as one might be tempted, to find a set of mode functions simply by diagonalizing the measured variance matrix, since the SO(2N ) transformation corresponding to diagonalization is generally not a canonical transformation (see, e.g. [8] for more details). On the other hand, in the special case of pure Gaussian states which are specified by N (N + 1) real parameters one can find a set of uncorrelated modes. In this case our procedure would be terminated after the first diagonalization. The diagonalization of pure Gaussian multimode states into uncorrelated modes has recently been studied in [21].
The procedure of operational construction of mode functions as described above is, of course, not the only possible. It is, however, very useful for finding the minimum subspace of mode functions sufficient for the pulse description. Other sets of mode functions can be selected such that the variance matrix takes some special shape, e.g., some of the "canonical" forms studied in [8].

VII. SEPARATION OF MODES
It may be very useful to separate individual nonmonochromatic modes which form the pulse. The possibility of obtaining nonclassical correlations of optical pulses by partitioning the pulse in the spectral region was suggested in [22]. Similarly as with spectral filtering, this approach is not necessarily the optimum one for obtaining maximum entanglement from the source.
The basic idea for obtaining the optimum separation is to send different (nonmonochromatic) modes into different channels by using a special unitary transformation among them. It was shown in [23] that any discrete unitary operator can be constructed interferometrically. For this purpose we suggest to apply a scheme as in Fig. 4. In the first step one decomposes the pulse into quasimonochromatic components. The frequencies of these components are then shifted by means of acustooptical modulators (AOMs) so that each channel has the same central frequency. The channels then interfere on a 2N -port consisting of beam splitters and mirrors. By a proper choice of the 2N -port parameters one can man-age that most of a pulse in the f k (ω) mode leaves the 2N -port in the kth output channel. One can also optimize the 2N -port parameters to prepare a set of optimally entangled K modes (multipartite entanglement), sent to different channels. This would be a generalization of the proposal of [22] for partitioning soliton pulses to generate entangled states.
Let us note that the spectral filtering of solitons used to produce photon number squeezing [19] is a special kind of mode separation. However, as shown in Sec. V, in this case the separation is not optimized with respect to all relevant degrees of freedom. Therefore, better results can be expected in our general approach.

VIII. DISCUSSION AND CONCLUSION
To summarize our results we can state that: (i) For a complete description of the experimental results measured on optical solitons it appears that a very small number of nonmonochromatic modes is enough. The results cannot be described by means of a singlemode field, but already three nonmonochromatic modes were sufficient to reproduce the experimental data of [4]. With our approach it is easy to interpret the "butterfly" pattern of the measured photon covariances: three different non-monochromatic modes overlap and contribute with their fluctuating in-phase quadratures to the photon statistics (see Fig. 2).
(ii) Provided that the pulses are Gaussian (which seems to be a relevant assumption for most of the experimental situations), a complete quantum description of the pulse can be done by means of a multimode variance matrix. The elements of this matrix can be determined by means of homodyne detection with specially shaped local oscillator pulses.
(iii) Concrete form of the mode functions is a matter of choice. If the aim is to find the smallest number of modes, an operational method for constructing the mode functions is provided. In this case the first mode contains the coherent amplitude, whereas the rest of the modes have zero mean fields and their field variances decrease with increasing mode index. By selecting the minimum set of modes, one can substantially reduce the number of parameters necessary for a complete description of the physical situation.
(iv) Knowledge of the mode structure of the pulse and of the corresponding quantum state will be very useful for quantum information purposes: one can select the optimum shape of the local oscillator pulse to detect maximum squeezing or entanglement. This optimum finding is very straightforward and much faster than adaptive algorithms as in [20]. One can also relatively easily study the influence of the medium on the propagated pulses and on the quantum information they carry. By measurement of the multimode quantum state of the input and output pulses we can find the von Neumann entropy of the states. This would enable us to tell whether the observed pulse deformation corresponds most probably to a unitary evolution or rather to decay and dephasing, possibly caused by an eavesdropper.
(v) In principle, one can also separate individual modes to match the requirements of the pulse user, e.g., to prepare a single mode maximally squeezed field, or to extract a maximally entangled two-mode field, etc.
(vi) The pulse separation is generalization of spectral filtering which has also been used for observation of photon number squeezing [19]. We have shown that spectral filtering can never beat the coherent amplitude modulation approach in achieving better photon number squeezing. In our approach, we can understand why spectral filtering can help with squeezing, but we can also see its limitations.
(vii) Our approach of multimode description of solitons is different from the approach by Haus and Lai [7] (see App. F for more details). The formalism of [7] was developed to solve an idealized quantum soliton equation and to study the soliton dynamics. Its essential feature is the decomposition of the field into the "soliton" and "continuum" part, where the soliton field is completely described by four operators. Our approach ignores the dynamics (which is rather complicated in real life), focusing on the phenomenological description of the pulse. The four soliton operators of [7] can be expressed by means of our quadrature operators, but they alone appear to be insufficient for description of the observed phenomena. It is necessary to work with the complete set of "soliton" and "continuum" operators if one wishes to describe the observed pulses using the formalism of [7]. discussions. T.O. thanks to Prof. G. Leuchs for his kind hospitality at the Friedrich-Alexander University in Erlangen.

APPENDIX A: SIMPLIFIED CALCULATION OF TRANSFORMED MEAN QUADRATURES AND VARIANCES FOR FEW EXCITED MODES
Let us assume that only the first N modes of theb k system are occupied, while the rest is in vacuum, |0 k for k > N . The quadrature means and variances are thus and V µk,µk = 1 2 for k > N, The transformations (27) and (29) can then be written asξ and To get (A5) from (A4), the orthogonality and completeness relation (20) was used. As a special case, let us consider the variance matrix elements V ′ xj,xj ′ when the mode functions f k (ω) are real, as in Sec. III B. Eq. (A5) then becomes An N -mode Gaussian state with mean quadratures µ and variance matrix V has the Wigner function W (µ), where V −1 is the matrix inversion of V . Thus, a Gaussian state is fully determined by 2N 2 + 3N real parameters (2N mean quadratures plus N (2N + 1) independent elements of the symmetric matrix V ). Generally, a quadrature moment μ n k can be calculated as the integral of the Wigner function The symmetrical two-variable moments 1 2 {μ k ,μ ′ k ′ } and can be calculated as whereδ µk,µ ′ k ′ ≡ 0 if µ k and µ ′ k ′ commute with each other, andδ µk,µ ′ k ′ ≡ 1 if µ k and µ ′ k ′ are conjugate variables. Using Eq. (B1) one can analytically evaluate the integrals and express the quadrature moments by means of the parameters µ k and V ′ µk,µ ′ k ′ . We obtain and Note that although (B5), (B6), and (B9) are generally valid for all states, the equalities (B7), (B8), and (B10) only hold for Gaussian states. Corresponding expressions can be found also for the moments of the quadratures ξ(ω). Let us assume a nonmonochromatic mode defined by the discrete function f (ω k ), k = 1 . . . N , N k=1 |f (ω k )| 2 = 1. Let the state of this mode be Gaussian. We are interested in the photon number correlations between different frequency channels.
Theorem: All non-zero elements of a multi-mode photon correlation matrix of an effectively single-mode Gaussian state have the same sign. They are negative (positive) iff the single-mode state is sub(super)-Poissonian. Proof: The mode function f (ω k ) defines the first row of the unitary transformation matrix U . Let the parameters of the single-mode Gaussian state beX,P , V XX , V P P , and V XP . The multi-mode parametersx k , p k , V xk,xl , V xk,pl , etc., can be calculated by means of the simplified summation as in (A3) and (A5). The sign of the off diagonal element of the correlation matrix (36) is determined by the sign of the covariance (33). After some algebra with applying Eqs. (A3), (A5), and (34) we arrive at Thus, we can see that the sign of the non-zero elements does not depend on the arguments ω k and ω l , but is fully determined by the expression in the square brackets. If we denote the photon number in the single mode byn, we find that the quantity ∆n 2 − n is equal to the expression in the square brackets. This quantity is negative for sub-Poissonian states and positive for super-Poissonian states by definition, QED.

Showing that
where c is a matrix with the elements c kl . Eq. (D3) could only be valid if there exists some eigenvalue c λ of c such that c λ > 1. Let us assume that such an eigenvalue does exist. Let u k be the elements of the corresponding eigenvector, i.e., i.e., Let us define a function q(ω) as Since u k = q(ω)f k (ω)dω, one has i.e., which cannot happen for any function q(ω) since |c 2 (ω)| ≤ 1. Thus, A 2 ≤ 1, QED.

APPENDIX E: ELIMINATION OF COHERENT AMPLITUDES OF A MULTI-MODE GAUSSIAN STATE
Let us assume that N modes of the b-system are excited and the rest is in vacuum. Then, one can redefine the N modes such that only one of them (say modeb 1 ) has nonzero values ofX 1 ,P 1 , whereasX m = 0,P m = 0 for m > 1 (the variance matrix elements corresponding to these modes are, however, generally non-zero).
Proof In [7] (see also [12,13]) it is suggested to describe the quantum fluctuations of the Nonlinear Schrödinger equation (NSE) soliton by writinĝ a(x, t) = a o (x, t) + ∆â(x, t), with k, l = n, θ, x, p (note that f k themselves do not form an orthogonal set; explicit form of functions f k and f k can be found in [12]). Writing ∆â as a sum of hermitian operators ∆â = ∆â (1) + i∆â (2) one finds that ∆n = f * n (x)∆â (1) ∆θ = i f * θ (x)∆â (2) The relationship of this approach to our scheme can be easily examined if one assumes the first four mode functions f k (ω) in the form These functions were obtained by an orthogonalization procedure from the Fourier transformed functions f n,θ,x,p with ω o = 2c/(πξ). With this choice of the mode functions one finds that the soliton perturbation operators can be expressed by means of the quadraturesX 1 ,X 2 , P 1 ,P 2 ,P 3 , andP 4 as Assuming vacuum fluctuations ofX k andP k , X 2 k = 1/2, P 2 k = 1/2, one finds that the uncertainty products of the soliton perturbation operators are ∆n 2 ∆ θ 2 = 1 3 + π 2 36 ≈ 0.675, which are larger than the minimum uncertainty value 1/4 following from the commutation relations [∆n, ∆θ] = i and [∆x, n o ∆p] = i. This result corresponds exactly to that of [7,12]. Our interpretation of the result is that the operator ∆θ is not purely a conjugate of ∆n, but it contains an admixture of an operator commuting with ∆n. This admixture [in (F16) proportional toP 3 ] increases the noise above the minimum uncertainty limit.
Similar interpretation holds for the pair ∆x, n o ∆p.
From Eqs. (F15)-(F18) one sees that from the statistics ofX 1,2 andP 1,2,3,4 one can determine the statistics of the soliton perturbation operators ∆n, ∆θ, ∆x and ∆p. This does not hold vice versa: knowledge of the four soliton operators is not enough to determine the six quadratures.