1. Introduction

In physical systems governed by nonlinear dynamics, any tiny variation of the input may lead to a significant variation at the output. This is particularly relevant in strong nonlinear fibers and waveguides, as pointed out in recent studies on the stability of the supercontinuum generation with long (picosecond) or even continuous-wave light [1], or the generation of optical rogue waves [2], for example. A proper description of the stochastic nature of the input light and its nonlinear propagation is a topic of increasing relevance. Among the different existing approaches, the theory of optical coherence is a well-established formalism whose aim is to study the ensemble statistics of the stochastic light in terms of correlation functions [3].

Since the discovery of incoherent spatial solitons [4], the propagation of spatially partially coherent light in nonlinear media has been widely studied using different theoretical methods [5–7]. In these studies, the nonlinear response of the medium is assumed to be noninstantaneous, in which case a closed-form propagation equation for the second-order coherence function is found. However, in standard optical fibers and other highly nonlinear waveguides of interest, the basic Kerr nonlinearity causing self-phase modulation is practically an instantaneous effect and the same assumptions as in the spatial domain studies cannot be normally used for modeling partially coherent pulses.

While the propagation of arbitrary partially coherent pulses in linear dispersive media is well understood [8–11], the extension of such studies to the nonlinear regime remains a challenge owing to the fast response of the Kerr effect. Approximate analytical expressions for the average properties or correlation functions of partially coherent pulses propagating in nonlinear media have been derived under some restrictions [12, 13], and the spectral changes caused by selfphase modulation have been studied in the special cases of stationary [14], quasi-stationary [15] and cyclo-stationary [16] fields. In the general case, the propagation of the coherence properties of pulses in instantaneous nonlinear media is governed by an infinite hierarchy of moment equations for which a closure can be found if the fields obey Gaussian statistics [6]. This kind of statistics is commonly encountered in light sources obeying a chaotic behaviour, such as flash lamps, amplified spontaneous emission sources, or blackbody radiation. However, the statistical properties are generally not conserved when the light propagates in strongly nonlinear media and thus the Gaussian assumption used to achieve a closed-form propagation equation is not physically valid [15, 17]. Moreover, most of the common ultrashort laser pulses deviate from the Gaussian statistical distribution as they are more stable in amplitude than thermal sources.

In this paper we present an alternative modeling approach based on numerical Monte Carlo simulations, where a set of field realizations are propagated individually through the system and the statistical properties of the ensemble are studied in the output. The use of Monte Carlo method has been suggested for studying the propagation of partially coherent radiation in nonlinear media [17], but in the case of pulsed fields it has been applied only in some earlier works that assumed chaotic fields with specific Gaussian-shaped intensity and correlation functions and Gaussian statistics [12, 18]. Here we propose a general method for the construction of an ensemble of pulses that represent an arbitrary input second-order coherence function and present a procedure to create ensembles that follow non-Gaussian statistics, which allows a more realistic description of non-thermal fields such as laser pulses. Furthermore, we study numerically the convergence of the results in order to estimate the required number of samples in this kind of simulations. By propagating the pulses of the ensemble individually, i.e., coherently, through the nonlinear fiber we can construct the corresponding temporal and spectral coherence functions in the output without restricting assumptions. This allows us to study the average intensities, spectra, and correlation properties of partially coherent pulses in strongly nonlinear media having an instantaneous response.

2. Theory

2.1. Basic coherence definitions

In this section we review the basic coherence concepts used throughout the paper. Realistic optical pulses contain stochastic variations, e.g., in the amplitude, phase, shape or temporal duration. The second-order statistical properties of an ensemble of plane-wave pulses are characterized by the mutual coherence function,

that measures the second-order correlation of the complex fields [3]. Here the angle brackets denote an ensemble average over the random field complex amplitudes U at time instants t ₁ and t ₂ in the plane z. The average intensity of the pulses is obtained from the mutual coherence function evaluated at a single space-time point as I(t,z) = Γ(t,t;z) = 〈∣U(t,z)∣²〉. When the mutual coherence function is normalized by the intensity as

we get the complex degree of coherence. When ∣γ(t ₁,t ₂;z)∣ = 0, the field realizations at the temporal points t ₁ and t ₂ are uncorrelated and the field is said to be fully incoherent. On the other hand, ∣γ(t ₁,t ₂;z)∣ = 1 indicates that the field is fully coherent and therefore complete correlation between the spatiotemporal points appears. Realistic optical fields, however, lie in between those extreme cases and are called partially coherent.

In the frequency domain, the corresponding quantity measuring the second-order correlation between angular frequencies ω is named the cross-spectral density,

where Ũ is the Fourier transform of the field complex amplitude. This function is related to Γ through the generalized Wiener-Khintchine theorem,

and to the average energy spectrum of the pulses by S(ω,z) =W(ω,ω;z). The complex degree of spectral coherence, μ(ω ₁,ω ₂;z), can be defined analogously to Eq. (2).

Other physical quantities, such as the coherence time and the spectral coherence width, can be derived from the above correlation functions. For instance, the coherence time is related to the counter-diagonal width of the complex degree of coherence γ [19], and it determines the time duration over which two pulses can be added coherently. The corresponding frequency-domain quantity can be used, e.g., to estimate the useful optical bandwidth to achieve phasestable superposition, which is relevant for pulse compression applications or generation of terahertz and microwave radiation through beating of the optical components [20]. This and further information of the properties of the pulses could also be obtained by studying directly the structure of their coherence functions, including the higher-order correlation terms. The spectral degree of coherence defined in this work accounts for the cross-correlation between any two spectral components of the signal. It is worth noting that this definition is different from the concept of “degree of coherence” typically used in supercontinuum studies [1], which accounts for the standard phase deviation for individual wavelength components. The coherence functions expressed in this work are defined in the framework of the formalism of the coherence theory [3] and complement the information available from these other studies.

2.2. Propagation in nonlinear media

Let us provide a brief description of the statistical evolution of partially coherent radiation in nonlinear instantaneous Kerr media. On one side, the propagation of coherent pulses in nonlinear media is governed by the well-known nonlinear Schrödinger equation (NLSE) [21],

where U = U(t,z) denotes in this case the slowly-varying pulse envelope, β ₂ is the groupvelocity-dispersion (GVD) parameter, and γ is the nonlinear coefficient of the material. Moreover, as usual, the time coordinate t is measured in the reference frame moving at the group velocity of the pulse. This form of the NLSE is valid for pulse widths T ₀ ≥ 1 ps. The accurate description of the propagation of shorter pulses in nonlinear media usually needs to take into account the effects of higher-order dispersion parameters and nonlinearities.

On the other side, the corresponding propagation equation for the coherence functions of partially coherent pulses can be derived by using the standard procedure [3]. Accordingly, taking the ensemble average of the difference of the two equations obtained by considering, first, the complex conjugate of Eq. (5) for U ₁ at t ₁ multiplied by U ₂ = U(t ₂, z) and, second, multiplying Eq. (5) written for U ₂ at t ₂ by U ^* ₁ = U ^*(t ₁, z), we get

where now U describes the random complex amplitude. Equation (6) determines that the evolution of the coherence function upon propagation depends on the fourth-order correlations of the field. In the same way, it is possible to derive a propagation equation for the fourth-order correlations, which in turn depends on the sixth-order correlations of the field. Continuing on this process, an infinite hierarchy of propagation equations could be obtained.

The common approach to find a closure for the hierarchy of the propagation equations is to assume that the fields obey Gaussian statistics. In this case all the higher-order coherence properties can be expressed in terms of the second-order correlation function [3]. Unfortunately, the Gaussian property is not generally preserved during the propagation in nonlinear media [15,17], since the above assumption is valid only when dispersive effects are dominant. Thus, simple and general propagation laws are not easily found for the statistical description of nonstationary partially coherent pulses.

3. Formulation of the method

Instead of finding a propagation equation for the mutual coherence function, statistical properties of partially coherent pulses in nonlinear media can be studied numerically using waveoptical Monte Carlo method [17]. The numerical procedure is based on describing the partially coherent field by random sample pulses. A similar approach has been used earlier to model the nonlinear propagation of partially coherent input pulses with Gaussian-shaped intensity and correlation functions and Gaussian statistics [12,18]. Here we provide the general guidelines for applying the method to more realistic situations, i.e., to study pulses with arbitrary coherence properties and to create ensembles whose statistics are non-Gaussian. We thus pave the way to study the evolution of cross-correlation functions of noisy (although not necessarily chaotic) radiation, such as realistic laser pulses, in strong and instantaneous nonlinear media.

3.1. Simulation steps

In the current problem the nonlinear material is assumed to be regular and only the partially coherent input field features random characteristics. Hence, the simple simulation scheme consists of three main steps:

In the first step the coherence properties of the input fields must be known in advance either by measurements [22], or by choosing some appropriate model field to describe them properly. The construction of the corresponding random pulse ensemble can be carried out in different ways. In the next subsections we present one straightforward option for doing this.

The propagation of the pulses through the optical system, such as a nonlinear fiber, can be done by using any standard algorithm to solve Eq. (5), since despite there are random fluctuations from pulse to pulse in the ensemble, every single pulse is completely coherent. For the numerical examples of this article, we have used the basic split-step Fourier method [21].

Finally, in the last step the functions describing the coherence properties of the propagated pulse ensemble are reconstructed. In accordance with Eq. (1), the mutual coherence function of an ensemble of L random pulses U_i(t) is given by

The rest of coherence functions can be easily obtained using the definitions presented in Section 2.

3.2. Ensemble construction

The construction of the ensemble of random pulses is based on the coherent mode representation of the input coherence function [3, 23, 24]. Since the mutual coherence function is Hermitian and non-negative definite, it can be represented as an absolutely and uniformly converging series,

where the coefficients λ_n and functions ϕ_n(t) are the eigenvalues and eigenfunctions, respectively, of the Fredholm integral equation,

For some special coherence functions, such as Gaussian Schell-model fields [3], there exist analytical solutions for λ_n and ϕ_n(t), but in general they must be worked out numerically.

Due to the properties of the second-order coherence function, the eigenvalues are real and non-negative and the eigenfunctions orthonormal. Since the eigenfunctions form an orthonormal basis, any field realization can be represented as a series such that

This representation is also known as the Karhunen-Lóeve expansion of the random pulse U(t) [3]. Now if a_n are uncorrelated random variables satisfying the expression

it can be shown that an ensemble of such pulses corresponds to the initial mutual coherence function. However, it should be noted that the choice of the coefficients a_n fixes the higherorder statistics of the pulses. In other words, it is possible to have different ensembles providing the same Γ(t ₁,t ₂) but different higher-order correlations. In the next subsection we will provide further discussion on one possible way to define suitably the coefficients a_n.

The number of λ_n with significant value depends on the degree of coherence of the field [23]. Thus, when the coherence of the input field decreases, an increasing number of terms must be included in the sum of Eq. (10) in order to generate reliable samples. Alternatively, one could create a set of sample pulses in the frequency domain based on the eigenmodes obtained by representing the cross-spectral density in a similar manner to the sum in Eq. (8). With the same random coefficients a_n, these frequency-domain realizations would be equal to the Fourier transforms of the time-domain samples given by Eq. (10).

The above way to construct random field realizations is based on the elementary properties of the coherent mode expansion [3, 24]. However, we stress that our approach is not equivalent to the propagation of the coherent modes of the field separately through the nonlinear media, which is in fact the conventional approach when studying both linear optical systems [3] and, related to the self-consistent multimode theory, partially coherent fields in non-instantaneous nonlinear media [5]. Unlike the well-defined mode functions, each random realization contains individual fluctuations affecting the nonlinear propagation. The average effects of these fluctuations can be related to the coherence functions only when a sufficiently large number of random realizations is propagated throughout the system.

A similar approach has been used before in the particular case of Gaussian Schell-model fields based on the analytical solution for the eigenfunctions ϕ_n [12]. At this point it is worth noting that numerically the method can be extended to a wider class of fields with different average intensity and correlation profiles. At the same time, other ways to simulate partially coherent field realizations have been reported in some recent papers [25, 26]. In all these earlier works, only pulse ensembles obeying Gaussian statistics have been considered. We again emphasize that this choice corresponds to chaotic light produced by conventional light sources, such as gas discharge or filament lamps and thermal cavities [27]. However, the statistical properties of pulsed laser light are quite different and thus such ensembles do not provide a proper model for laser pulses. In the following, we present an option for the generation of an ensemble that follows non-Gaussian statistics, which could be used to model the coherence properties of noisy laser pulses.

3.3. An ensemble with non-Gaussian statistics

Not even the pulses emitted by laser sources have ideal and well-defined characteristics, since different noise sources may cause fluctuations in the phase, amplitude, central frequency, pulse width, etc. Small variations in these parameters can have a significant effect when pulses propagate in strong nonlinear media. One of the sources of such variations is the technical noise that results from the fluctuations of the laser input power. However, this kind of noise can be reduced using stabilization techniques and therefore it is often neglected in studies of nonlinear pulse propagation [1, 28]. On the other hand, in modern mode-locked fiber lasers, the fluctuations in the pulses may be stronger than predicted by simplified models that only consider quantum-limited shot noise [29]. While the exact structural form of each noisy pulse is difficult to predict, their average and correlation properties can be studied by the coherence functions that easily allow us to model pulses with different amounts of noise, or, in other words, with different degrees of coherence.

If the coefficients a_n in Eq. (10) are chosen to be complex Gaussian random variables, the resulting pulse ensemble obeys Gaussian statistics [12, 30]. However, we note that with a different choice of the coefficients, an ensemble with different higher-order coherence properties is constructed. For this purpose, random coefficients that satisfy the condition in Eq. (11) can be defined as

where φ_n are random phases uniformly distributed in the range [−π,π]. If we interpret the modes of the coherent-mode presentation in Eq. (8) as natural modes of oscillations of the source [3], the above choice of a_n leads to random sample pulses that are linear combinations of such elementary fields with random phases. It can be shown that the resulting samples follow non-Gaussian statistics and consequently provide a way to represent coherence functions of non-thermal fields. The field realizations constructed following the above prescription have in general random fluctuations both in the phase and amplitude distributions, but each random pulse realization U_i(t) carries the same amount of energy. Thus, this model could be used to describe a set of laser pulses without technical noise.

As an example of the characteristics of the random pulses generated in such a manner, Fig. 1 shows the temporal amplitude of different samples corresponding to Gaussian Schell-model pulses (GSMPs) with two coherence times. In all our calculations, the number of the eigenfunction terms included in the sum of Eq. (10) depends on the number of the corresponding λ_n with numerically significant values, which varies with the degree of coherence. The average amplitude and envelopes of a large set of realizations are also shown in the figure. The phase of each pulse also shows random, continuous variations that are more significant the less coherent the pulses are. Similar characteristics can be recognized when pulses are plotted in the frequency domain.

 

Fig. 1. Amplitude of 30 random realizations (red, thin lines) corresponding to 10 ps GSMPs with non-Gaussian statistics fixed by Eq. (12) and coherence times: (a) T_c = 50 ps, and (b) T_c = 10 ps. The average amplitude (black, thick line) and the envelopes (blue, dashed line) are shown for a set of 20000 random realizations. The absolute value of the corresponding mutual coherence functions are shown in the insets. As expected, each individual realization has the same energy.

Download Full Size | PDF

4. Accuracy of the simulations

The above Monte Carlo simulations try to reproduce the partially coherent field by means of a set of random pulse samples properly selected. The accuracy of such a replica depends on the number of the samples included in the ensemble. The higher the number of sample pulses, the closer their average and correlation properties match those corresponding to the input field. The error due to the finite number of samples in the ensemble can be estimated by studying the convergence of the results as the number of samples increases.

Let us consider an input field with a certain mutual coherence function Γ₀(t ₁,t ₂). A set of L random pulses defined by Eq. (10) provides, according to Eq. (7), a coherence function Γ_L(t ₁,t ₂). The mean square error ε_L between these functions is given by

Clearly, ε ² _L depends on the number of random pulses included in the ensemble. Due to the properties of the Fourier transformation, the error is the same when it is calculated in the spectral domain. An example of the convergence of the procedure is shown in Fig. 2. Again the characteristics of the input fields are the same as those considered in Fig. 1. It is apparent that a larger number of samples is needed to reproduce the coherence function accurately when we deal with less coherent fields.

Typically in Monte Carlo simulations the error decays at the rate L ^−1/2 for large enough values of L. The test ensembles shown in Fig. 2 also follow this behavior. The deviations from the above decay behavior are more significant when the number of samples is relatively small and we deal with less coherent pulsed light. In the latter case, due to the statistical nature of the process, every random realization may suffer strong fluctuations from the average and the resulting error is different for distinct ensembles even with the same number of samples. In other words, it is clear that the mean value of a few sets of realizations with a fixed value of L lies, in principle, closer to the fitted curve than a single one. Moreover, note that all numerical calculations are worked out with standard double precision.

 

Fig. 2. Variation of the error function for random samples fixed by Eq. (12) corresponding to 10 ps GSMP input fields with coherence time equal to T_c = 10 ps and T_c = 50 ps. The curves fitted to the computed data, which follow a L ^−1/2 decay rate, are plotted in solid line.

Download Full Size | PDF

When light propagates in nonlinear media, the coherence properties of the input field are typically changed. Thus it is also necessary to find out the required size of the ensemble of random pulses to represent accurately the propagated field at the output. However, we do not know the exact form for the propagated coherence function in advance and the error must be estimated against the coherence function determined from a chosen large number of samples. In mathematical terms, the Γ₀ in Eq. (13) must be replaced by certain Γ_M where M is a sufficiently large number.

In order to study the convergence of the error function after propagation, we again consider GSMPs with temporal duration T ₀ = 10 ps, coherence time T_c = 10 ps, and central wavelength λ ₀ = 1550 nm, but in two different nonlinear Kerr media. In the first case (fiber 1), the GVD parameter is chosen to be β ₂ = 50 ps²/km, the nonlinear coefficient γ = 0.5 (Wkm)⁻¹, and the average peak power of the pulses P ₀ = 1 W. Thus the propagation lengths where the dispersive and nonlinear effects become important, defined correspondingly as L_D = T ² ₀/∣β ₂∣ and L_NL = 1/(γP ₀), have the same value. Consequently the parameter N ², or the soliton number, becomes N ² = L_D/L_NL = 1 [21]. The convergence of the error function when the pulses are propagated distances z = 2,10, and 20 km is shown in Fig. 3(a). In the other case (fiber 2), the peak power is changed to P ₀ = 50 W and the fiber parameters are assumed to be β ₂ = 20 ps²/km and γ = 10 (Wkm)⁻¹ corresponding to a highly nonlinear fiber. Now the nonlinear effects are severe and the soliton number is N ² = 2500. The pulses feature strong changes after the propagation in such a medium, and the convergence of the error function is shown for values z = 5, 15, and 30 m in Fig. 3(b). According to the previous paragraph, in all the cases the different coherence functions are compared to Γ_M with M = 20000.

The results in Fig. 3 show that the convergence of the method depends strongly on the properties of the nonlinear media. In fiber 1, the changes in the output coherence function are relatively small and the propagation distance does not substantially affect the required number of samples. In contrast, in fiber 2, with higher nonlinearities, we require an increasing number of samples to describe the coherence function accurately as the pulses propagate. Furthermore, the error is more significant in Fig. 3(b) for a given number L due to the intricate structure of the coherence function at the output. We recognize that in both fibers the decay of the error agrees closely with the predicted L ^−1/2 rate. We can reach the same conclusions when we deal with pulses with a higher degree of coherence, although smaller sets of samples are sufficient to achieve the same reliability level for the results. For comparison, we note that in Ref. [12] the numerical simulations were carried out using only 100 samples, which clearly would not give accurate results for the situations considered here. Broadly speaking, the convergence of the error function should be checked individually at each particular situation where the simulation method is used.

 

Fig. 3. Variation of the error function for random samples fixed by Eq. (12) corresponding to 10 ps GSMPs with coherence time T_c = 10 ps propagated in a nonlinear medium with soliton number: (a) N ² = 1 (fiber 1), and (b) N2 = 2500 (fiber 2). The curves fitted to the computed data, which follow a L ^−1/2 decay rate, are plotted in solid line.

Download Full Size | PDF

5. Numerical examples

Finally, we show some illustrative examples corresponding to the modifications suffered by the output coherence functions when partially coherent pulses propagate through a highly nonlinear fiber. In this section, the fiber and pulse parameters are assumed to be the same ones as in Fig. 3(b).

First, Fig. 4 illustrates the characteristics of a set of propagated random pulse realizations in the frequency domain for different propagation distances. When the pulses propagate, new frequencies are generated by self-phase modulation and the broadened spectrum of each random pulse shows also the oscillations typical to this phenomenon. As expected, these oscillations are flattened out in the average spectral amplitude of the partially coherent pulse ensemble.

 

Fig. 4. Spectral amplitudes of 30 random realizations (red, thin lines) fixed by Eq. (12) corresponding to 10 ps GSMPs with coherence time T_c = 10 ps propagated in fiber 2 distances: (a) z = 0, (b) z = 15, and (c) z = 30 m. The average amplitude (black, thick line) and the envelopes (blue, dashed line) are shown for a set of 20000 random realizations.

Download Full Size | PDF

Assuming higher coherence for the input pulsed light, more oscillations are observed in the average spectrum. In the time domain, the changes in the propagated random pulses are less pronounced.

Based on the output random sample pulses, we also study the evolution of the coherence functions corresponding to the pulse ensemble. In Fig. 5 (Media 1), the evolution of the mutual coherence function, cross-spectral density, and temporal and spectral degrees of coherence are shown for the same fiber and pulse parameters as above. The diagonal of the plots corresponding to Γ and W show the variation of the average temporal and spectral intensities, while the counter-diagonal of γ and μ are related to the central temporal and spectral coherence widths of the pulses. At the beginning of the propagation, intricate changes in the coherence functions are noticed and the coherence decreases strongly as the pulses propagate further. Note that some apparent coherence is finally preserved only in the edges. This feature corresponds to temporal instants or frequencies where the pulses carry no significant energy and, then, such correlations are not relevant. For this simulation 5120 sample pulses were used, so the numerical error in the results can be estimated from Fig. 3(b) to be less than 5%.

 

Fig. 5. Coherence functions corresponding to 10 ps GSMPs with non-Gaussian statistics fixed by Eq. (12) and coherence time T_c = 10 ps after propagation in fiber 2 from z = 0 to z = 30 m. (Media 1).

Download Full Size | PDF

For comparison, we illustrate the influence of the input degree of coherence on the changes of the coherence functions. In Fig. 6 (Media 2) we show the evolution of similar input pulses but with longer coherence time, T_c = 50 ps, in the same fiber 2. The propagation in the nonlinear medium shows that the coherence strongly decreases also in this case, even though the stochastic fluctuations in the input pulses are relatively small. Using similar convergence analysis as discussed in the previous section, we find that the chosen number of 2048 sample pulses provides smaller error, less than 2%, than in the previous example, as could be expected also based on Fig. 2.

 

Fig. 6. Coherence functions corresponding to 10 ps GSMPs with non-Gaussian statistics fixed by Eq. (12) and coherence time T_c = 50 ps after propagation in fiber 2 from z = 0 to z = 30 m. (Media 2).

Download Full Size | PDF

Similar modifications in the characteristics of partially coherent fields propagating in nonlinear media have also been reported in some earlier studies. For instance, the smoothing of the characteristic peaks of the self-phase modulated spectrum that can be seen in the average spectral amplitudes of Fig. 3 has also been reported in the case of thermal fields within the framework of cyclostationary processes [16]. The compression of the coherence time that is shown in Figs. 5 and 6 (Media 1 and Media 2) has been noticed for incoherent Gaussian fields using analytical and numerical studies as well [14,17,18]. By studying the changes of the complete coherence functions we can get a wider view on the possible effects of the propagation in nonlinear media, such as the local variations in the degree of coherence that are especially evident in Fig. 6 (Media 2).

6. Conclusions and outlook

We have presented a numerical technique for describing the evolution of the coherence functions when partially coherent pulses propagate in instantaneous nonlinear media. The method is based on Monte Carlo simulations where an ensemble of random sample pulses that represents the partially coherent field is propagated pulse by pulse through the medium. The output coherence functions are then reconstructed based on the statistical properties of the propagated ensemble. We have carried out a systematic study of the numerical accuracy of the method that depends on the number of random pulses included in the ensemble and the nonlinear properties of the medium at issue.

Our method can be applied for the description of the propagation of arbitrary input coherence functions, including those corresponding to non-thermal light such as noisy laser pulses, in media with strong nonlinearities and instantaneous response, which hitherto has remained challenge. The numerical examples show that the nonlinear propagation in such media generates intricate changes in the coherence functions and the degree of coherence of the input light strongly affects the coherence characteristics of the output. We would like to emphasize that the propagated coherence functions contain information on a large variety of the basic characteristics of the pulses, such as average intensity, average spectrum, coherence time, spectral coherence width and other more detailed statistical properties of the ensemble, like correlations between two particular frequencies.

In this paper, we have discussed one possible way to create partially coherent pulse ensembles that correspond to a specified second-order mutual coherence function and do not necessarily follow Gaussian statistics. If all the higher-order coherence properties of the pulses are known, more rigorous rules for defining the ensembles could be derived. In future work we aim to study further the influence of the higher-order correlations in the nonlinear propagation of partially coherent pulses. In addition, even though in the examples of this paper we have only considered media with second-order dispersion and the basic Kerr-type nonlinearity, the same analysis method is directly applicable to more complicated nonlinear propagation problems, such as supercontinuum generation. Furthermore, our method can easily include, if needed, the spatial transverse coordinates of partially coherent fields.