Numerical Study of Mid-IR Ultrashort Pulse Reconstruction Based on Processing of Spectra Converted in Chalcogenide Fibers with High Kerr Nonlinearity

: Ultrashort optical pulses play an important role in fundamental research and applications. It is important to have reliable information about pulse parameters such as duration, intensity proﬁle, and phase. Numerous methods for characterizing pulses in the near-IR range have been well developed by now. However, there is a challenge with pulse measurement in the mid-IR, which is largely related to the underdeveloped component base in this spectral range. We investigate by means of numerical simulations a simple method of pulse reconstruction applicable in the mid-IR. The method is based on measuring and processing only the initial pulse spectrum and two converted spectra in elements with Kerr nonlinearity for different B -integrals characterizing nonlinear phase accumulation. The hardware implementation of the proposed method is very simple. This method requires only a one-dimensional data set, has no moving parts in the optical scheme, and allows for working with high-energy as well as low-energy pulses. We propose a novel simple, efﬁcient, noise-tolerant algorithm for data processing that assumes spectral phase approximation by a polynomial function. We demonstrate numerically the reconstruction of mid-IR ultrashort pulses, namely 3 µ m wavelength pulses, using commercial chalcogenide As 2 S 3 -based glass ﬁbers as nonlinear elements.


Introduction
Ultrashort optical pulses with durations ranging from a few fs to one hundred ps play an important role in fundamental research and various applications. Reliable information about pulse parameters such as duration, intensity profile, and phase is needed in many tasks. At present, numerous methods for characterizing pulses in the near-IR range are quite well developed. Note that the widely used methods with commercial hardware and software implementation are based on second harmonic generation. These are ACF (Autocorrelation Function), SHG FROG (Second-Harmonic Generation Frequency-Resolved Optical Gating), and SPIDER (Spectral Phase Interferometry for Direct Electric-field Reconstruction) [1][2][3]. However, for these methods, there are ambiguities associated with the direction of the time axis and the sign of the phase and limitations associated, for example, with the phase-matching bandwidth of the used nonlinear crystals. Therefore, quite a number of modifications and original methods were proposed with various advantages for specific purposes [1,2,[4][5][6][7][8][9][10][11][12], for example, to avoid ambiguities, increase the temporal and/or spectral resolution and expand the characterized range, measure signals with a complex spectral-temporal structure, and avoid information loss due to averaging over many (sometimes unstable) pulses.
But there exists a challenge of pulse measurement in the mid-IR, which is related to the underdeveloped component base in this spectral range. Thus, the development of simple, reliable, and inexpensive methods for measuring the intensity profile and phase of mid-IR ultrashort pulses is of great importance, especially for low pulse energies. The present study aims to address this issue.
Recently, a method for measuring the shape and phase of ultrashort pulses based on the processing of self-modulated spectra was proposed, which makes it possible to reconstruct signals without ambiguities and without limitation on their spectral width [13,14]. In the current contribution, we continue the development and numerical study of the method of pulse reconstruction based on fibers with high third-order Kerr nonlinearity. We propose a novel efficient algorithm for data processing and demonstrate the applicability of the method for measuring mid-IR ultrashort pulses, namely 3 µm wavelength pulses, using commercial chalcogenide As 2 S 3 -based glass fibers as nonlinear elements. Such fibers are transparent up to 6.5 µm and have large nonlinear Kerr coefficients, therefore, they are very suitable for the considered problem. The method is based on spectral measurements. Only the initial pulse spectrum and two converted spectra in elements with Kerr nonlinearity are required. These two converted spectra must have different B-integrals characterizing nonlinear phase accumulation. Moreover, there is no need to know the exact values of the B-integrals, only their ratio should be known [13,14]. The hardware implementation of the method is simple and involves the use of elements with third-order rather than second-order nonlinearity. This method requires only a one-dimensional data set (unlike the FROG family of methods), has no moving parts in the optical scheme (unlike most correlators and FROG implementations), and allows working with high-energy as well as low-energy pulses, which makes it promising in terms of potential implementation in devices of integrated photonics and fiber optics [14]. The earlier proposed reconstruction algorithms are based on the Gerchberg-Saxton algorithm [13,14]. However, these Gerchberg-Saxton-like algorithms work correctly if the pulse phase is not large (does not exceed π) [13]. In the case of a large chirp, the global minimum may not be reached; instead, stagnation may occur at a local minimum of an error function characterizing the difference between the original and reconstructed spectra (iterative algorithms often loop around a local minimum). Previously, we applied elements of genetic algorithms to solve this problem [13], which significantly increased the calculation time but did not guarantee to find the global minimum. Here, we propose a very simple and efficient algorithm to overcome this challenge. The algorithm is based on exploiting a polynomial spectral phase and finding polynomial coefficients using an exhaustive search on a reasonable grid of parameters, which allows finding global and all local minima on the preset grid. Moreover, the found polynomial phase corresponding to the global minimum can be used to seed the iterative algorithm. This allows refining the "right" solution and attaining an almost perfect agreement between the original and reconstructed data, which is demonstrated through numerous simulations. Next, we calculate parameters for commercial chalcogenide As 2 S 3 -based glass fibers and demonstrate their application as nonlinear elements in the analyzed method with the developed algorithm for measuring mid-IR pulses through representative numerical examples.

Methods
The method for the spectral phase reconstruction investigated in this work is based on measuring spectra converted during the pulse propagation in nonlinear elements with Kerr nonlinearity. Such elements are short segments of a highly nonlinear fiber of length L. If the dispersion length L D is much larger than the fiber length (L D = T 0 2 /|β 2 | >> L, where T 0 is the characteristic pulse duration, β 2 is the quadratic dispersion coefficient [15]), the dispersionless approximation can be used to describe pulse propagation, and the evolution of the complex electric field amplitude E is described by a very simple equation [15]: where γ is the nonlinear Kerr coefficient of the fiber; t is the time; and z is the coordinate along the fiber. This equation can be solved analytically [15]. The pulse intensity where P is peak power. The complex field envelopes of the output pulses during propagation in fiber for two different values of the B-integral B 1 and B 2 are where E 0 (t) is the field at the input and E 1 (t) and E 2 (t) are the fields at the fiber output. Let us also introduce the spectral complex amplitudes: where f is the frequency counted from the central pulse frequency f 0 andF is the Fourier transform operator. The corresponding spectral phases are ϕ 0 , ϕ 1 , and ϕ 2 .
The main idea of the method is illustrated in Figure 1a. There is a pulse with known spectral intensity I 0 (f ) and spectral phase (ϕ real ) that should be reconstructed. Let us denote the converted spectral intensities at the output of the nonlinear fiber as I 1 (f ) and I 2 (f ) measured for two different values of the B-integral (B 1 and B 2 , respectively). We guess with some accuracy the spectral phase ϕ 0 is close to ϕ real (Figure 1a). Using the inverse Fourier transform operatorF −1 , we can calculate the corresponding field E 0 (t) in the time domain: Next, we can calculate the complex fields E 1 (t) and E 2 (t) using expressions (3) and (4).
Note that E 1 ( f ) 2 and E 2 ( f ) 2 will be close to I 1 (f ) and I 2 (f ), respectively. Moreover, the smaller the difference between ϕ 0 and ϕ real , the better match E j ( f ) 2 with I j (f ) (j = 1,2). To quantitatively characterize the coincidence of the real and reconstructed spectra, we introduce the error function: where N is the number of sampling points of the frequency interval and f q is the frequency at the q-th point. plitude ( ) and corresponding field amplitude E0(t) is calculated. Next, field amplitudes E1(t) and E2(t) are simulated using Equations (3) and (4) and the error function is estimated using Equation (7). An exhaustive search allows for finding quadratic and cubic coefficients giving global and all local minima on the preset grid. If necessary, it will not be difficult to implement this algorithm taking into account the search for a polynomial phase of a higher order. In this work, we purposely study representative cases in which, in addition to the global minimum corresponding to the solution, there are local minima.  Further, we can formulate the problem formally. We need to find a spectral phase that minimizes the error ∆. Note that in the described procedure it is necessary to set B-integrals which can be pre-assessed but not exactly known as a rule; however, the ratio B 2 /B 1 is known. We run the developed algorithms for various B 1 and select the value minimizing the error function. It will be demonstrated in Section 3.4 that the found B-integrals coincide with the original B-integrals. The same approach was successfully applied for the processing of experimental data using a Gerchberg-Saxton-like algorithm with the genetic algorithm elements in [13,14].
Here, to prevent possible looping of iterative algorithms around a local minimum of the error function and guarantee to find the solution near the global minimum, we propose a very simple, efficient, and fast algorithm based on approximating the spectral phase by a polynomial function and optimizing the corresponding polynomial coefficients by exhaustive search on a reasonable grid so as to minimize the error and find the optimal polynomial phase ϕ 0 opt . We consider the case of spectral phase approximation taking into account the quadratic and cubic contributions, which is the most important case from the practical point of view and is usually sufficient in the common case of pulse propagation in a dispersive medium [15]. In addition, this case allows a simple and intuitive graphical interpretation, providing visual control of seeking a solution. The scheme of the proposed algorithm is shown in Figure 1b. A polynomial spectral phase is constructed on a reasonable grid of the quadratic and cubic coefficients; spectral complex amplitude E 0 ( f ) and corresponding field amplitude E 0 (t) is calculated. Next, field amplitudes E 1 (t) and E 2 (t) are simulated using Equations (3) and (4) and the error function is estimated using Equation (7). An exhaustive search allows for finding quadratic and cubic coefficients giving global and all local minima on the preset grid. If necessary, it will not be difficult to implement this algorithm taking into account the search for a polynomial phase of a higher order. In this work, we purposely study representative cases in which, in addition to the global minimum corresponding to the solution, there are local minima.
Next, we propose to use the found approximate solution with a polynomial phase to seed the previously developed iterative algorithm. In this case, the iterative search for an improved solution occurs near the global minimum, so this hybrid algorithm is well suited for phase refining (finding the phase difference from the polynomial one). The scheme of one cycle of the iterative algorithm is shown in Figure 1c (a more detailed description can be found in [14]).

Reconstruction of "Ideal" Signal
First, we considered an "ideal" case. An initial pulse without any noise was set with a spectral phase containing only quadratic and cubic terms: where C 2 (0) = 6.5 ps 2 and C 3 (0) = −3.6 ps 3 . The spectral intensity profile was sech 2 -shaped; the spectral width was FWHM = 0.6 THz (full width at half maximum). We applied the proposed algorithm shown in Figure 1b. We performed an exhaustive search over the range of parameters −10 ps 2 ≤ C 2 ≤ 10 ps 2 , −10 ps 3 ≤ C 3 ≤ 10 ps 3 (with steps dC 2 = 0.1 ps 2 and dC 3 = 0.1 ps 3 ) in which, in addition to the global minimum corresponding to the true solution, a local minimum was also found. The required range of C 2 and C 3 depends on the pulse spectral width and the maximum expected pulse chirp (which depends on the net dispersion). Figure 2a-d display the numerically found solution (the computational time is a few seconds). As expected for such a case, the numerically found spectral phase ideally coincides with the initially specified one; therefore, the temporal intensity profile and temporal phase as well as the output spectra after conversion in the Kerr fiber with the given B 1 and B 2 for the reconstructed pulse also perfectly coincide with the original ones. The calculated error function ∆ is shown in Figure 2e, where the global minimum is labeled "I" and the local minimum is labeled "II". We also specially constructed a pulse with a spectral phase corresponding to the local minimum II. In this case, the error is much higher and there is a significant difference in the converted spectra plotted for the initial pulse and for the "wrong" reconstructed one.

Reconstruction of Signals with Spectral Noise
Next, we investigated the effect of spectral noise on the results of the algorithm execution. We assumed that the original spectral phase was still determined by expression (8) with the same polynomial coefficients as in the previous case (C 2 (0) = 6.5 ps 2 and C 3 (0) = −3.6 ps 3 ). Then, we constructed the original spectrum and the nonlinearly converted ones but random noise was added to all three spectra at a level of 5% of the maximum of I 0 (f ) and we applied the numerical algorithm for finding the polynomial phase (Figure 3a-d). In this case, the results of the reconstruction of the polynomial coefficient C 2 coincided with C 2 (0) and the found value C 3 = −3.8 ps 3 slightly differed from the precisely specified value C 3 (0) (Figure 3a-d). The error was ∆ = 0.0072. For comparison, we constructed the pulse for

Reconstruction of Signals with Spectral Noise
Next, we investigated the effect of spectral noise on the results of the algorithm execution. We assumed that the original spectral phase was still determined by expression (8) with the same polynomial coefficients as in the previous case (C2 (0) = 6.5 ps 2 and C3 (0) = −3.6 ps 3 ). Then, we constructed the original spectrum and the nonlinearly converted ones but random noise was added to all three spectra at a level of 5% of the maximum of I0(f) and we applied the numerical algorithm for finding the polynomial phase (Figure 3a-d).
In this case, the results of the reconstruction of the polynomial coefficient C2 coincided with C2 (0) and the found value C3 = −3.8 ps 3 slightly differed from the precisely specified value C3 (0) (Figure 3a-d). The error was Δ = 0.0072. For comparison, we constructed the pulse for the polynomial phase corresponding to the local minimum II (Figure 3f-i) with a higher error Δ = 0.0201. ((f-i), right column) "Wrong" solution constructed for the local minimum II of the error function ∆ (C 2 = 2.2 ps 2 and C 3 = 4 ps 3 ) shown in (e). (a,f) Intensity and phase in the time domain; (b,g) input pulse spectrum and spectral phase; (c,h) converted spectra at the output for B 1 = 1; (d,i) converted spectra at the output for B 2 = 2. Black curves calculated for the initial pulse and color (magenta and blue) curves correspond to numerically reconstructed pulses for subplots (a-d,f-i).
Next, we analyzed a large number of different pulses generated for polynomial spectral phases with randomly set values of the coefficients C 2 (0) and C 3 (0) . The spectral intensity profile was also sech 2 -shaped and the spectral width was FWHM = 0.6 THz. For each pulse, random spectral noise at a level of 5% of the maximum of I 0 (f ) was added to the original spectrum and to the converted spectra (for B 1 = 1 and B 2 = 2), and then the polynomial phase search algorithm was applied (with steps dC 2 = 0.1 ps 2 , dC 3 = 0.1 ps 3 ). The values of the randomly generated coefficients corresponding to each specific test signal (100 pulses in total) are shown by dots in Figure 4a. Figure 4b shows the error function ∆ calculated for each specific test signal, and Figure 4c demonstrates the difference between the reconstructed and original polynomial coefficients. Note that the error function is practically independent of the value of the coefficients C 2 and C 3 . We also visually checked that in all cases the solution corresponding to the original signal was found.  Next, we analyzed a large number of different pulses generated for polynomial spectral phases with randomly set values of the coefficients C2 (0) and C3 (0) . The spectral intensity profile was also sech 2 -shaped and the spectral width was FWHM = 0.6 THz. For each pulse, random spectral noise at a level of 5% of the maximum of I0(f) was added to the original spectrum and to the converted spectra (for B1 = 1 and B2 = 2), and then the polynomial phase search algorithm was applied (with steps dC2 = 0.1 ps 2 , dC3 = 0.1 ps 3 ). The values of the randomly generated coefficients corresponding to each specific test signal (100 pulses in total) are shown by dots in Figure 4a. Figure 4b shows the error function Δ calculated for each specific test signal, and Figure 4c demonstrates the difference between the reconstructed and original polynomial coefficients. Note that the error function is practically independent of the value of the coefficients C2 and C3. We also visually checked that in all cases the solution corresponding to the original signal was found.

Reconstruction of Signal with Spectral Phase Perturbations
Here, we demonstrate the use of a hybrid algorithm (Figure 1a,b) for the reconstruction of a pulse whose phase differs from the polynomial one. We considered the case in which a uniformly distributed random value (0.1π × random[−1, 1]) was added at each point fq to the polynomial phase defined by expression (8). Therefore, the spectra I1(f) and I2(f) are strongly modulated. Spectral noise was not added for this test case. Using an ex-

Reconstruction of Signal with Spectral Phase Perturbations
Here, we demonstrate the use of a hybrid algorithm (Figure 1a,b) for the reconstruction of a pulse whose phase differs from the polynomial one. We considered the case in which a uniformly distributed random value (0.1π × random[−1, 1]) was added at each point f q to the polynomial phase defined by expression (8). Therefore, the spectra I 1 (f ) and I 2 (f ) are strongly modulated. Spectral noise was not added for this test case. Using an exhaustive search for the spectral polynomial coefficients, we reconstruct the pulse (Figure 5a-d, left  column). Further, we used this pulse with the polynomial spectral phase to seed the iterative algorithm. The results of its execution are shown in Figure 5e-h (right column). The phase peculiarities were reconstructed very well. The error ∆ was reduced by a factor of 27 (from ∆ = 0.0269 to ∆ = 0.0010) and the original and found converted spectra coincided very well. For comparison, we run the iterative algorithm under other initial conditions. Sometimes we observed that the iterative algorithm could loop near a local minimum (the error functions were similar to the ones shown in Figures 2e and 3e). In these cases, the original and found converted spectra were not in good agreement; the errors ∆ were significantly higher than for the solution near the global minimum. Therefore, the use of the hybrid algorithm guaranteed finding the "right" solution.  ; (b,f) input pulse spectrum and spectral phase; (c,g); converted spectra at the output for B1 = 1; (d,h) converted spectra at the output for B2 = 2. Black curves calculated for the initial pulse and color (magenta and blue) curves correspond to numerically reconstructed pulses.

Reconstruction of Mid-IR Signals Using Chalcogenide Fiber
Next, we considered the possibility of implementing the method for measuring ultrashort pulses in the mid-IR range using chalcogenide fibers with high Kerr nonlinearity. As an example, we numerically analyzed the reconstruction of the pulses produced by 3 μm wavelength lasers (with a central wavelength in the 2.8-3.9 μm range). Such pulses are routinely generated, for example, in rare-earth ion-doped fluoride fiber lasers [16]. We proposed to use step-index arsenic sulfide (As2S3-based) glass fibers as nonlinear elements. Note that such fibers are produced by different companies and are commercially available, which will make it quite easy to implement this method in research laboratories.  ; (b,f) input pulse spectrum and spectral phase; (c,g); converted spectra at the output for B 1 = 1; (d,h) converted spectra at the output for B 2 = 2. Black curves calculated for the initial pulse and color (magenta and blue) curves correspond to numerically reconstructed pulses.

Reconstruction of Mid-IR Signals Using Chalcogenide Fiber
Next, we considered the possibility of implementing the method for measuring ultrashort pulses in the mid-IR range using chalcogenide fibers with high Kerr nonlinearity. As an example, we numerically analyzed the reconstruction of the pulses produced by 3 µm wavelength lasers (with a central wavelength in the 2.8-3.9 µm range). Such pulses are Fibers 2022, 10, 81 9 of 13 routinely generated, for example, in rare-earth ion-doped fluoride fiber lasers [16]. We proposed to use step-index arsenic sulfide (As 2 S 3 -based) glass fibers as nonlinear elements. Note that such fibers are produced by different companies and are commercially available, which will make it quite easy to implement this method in research laboratories. As an example, we considered IRflex Corporation fibers [17]. Based on technical specifications [17], we estimated the fiber parameters and then simulated pulse reconstruction using these parameters.
To calculate the nonlinear Kerr coefficient γ and the quadratic dispersion coefficient β 2 of arsenic sulfide glass fibers, we used the well-known method of finding the propagation constant β and electric fields of fundamental modes for the linear polarization by solving numerically the characteristic equation at different frequencies [15]. The refractive index was given by the Sellmeier equation based on measurements performed by Amorphous Materials Inc. ("AMTIR-6" [18]) and the numerical aperture was NA = 0.3 (corresponding to the difference between the core and cladding refractive indices of about 0.019). The calculated nonlinear coefficients γ and dispersion coefficients β 2 are plotted in Figure 6a,b, respectively, for fiber diameters of 5, 6.5, 7, and 9 µm. All the considered fibers had a large γ of about 100 (W · km) −1 (Figure 6a). The dispersion curves were almost flat functions of wavelength in the 2.8-4 µm range (Figure 6b). Their calculations were required for estimating L D .  Next, we simulated diverse 3 μm pulses with realistic parameters and reconstructed them using the hybrid algorithm, assuming that an As2S3 fiber with a core diameter of 6.5 μm serves for spectrum conversions. We purposely set spectral phases with leading polynomial contributions due to quadratic and cubic terms and added different regular perturbations. Representative examples are shown in Figure 7. Each column corresponds to an individual pulse with original data (black), data reconstructed with the proposed algorithm for searching an optimal polynomial spectral phase (magenta), and data reconstructed with the iterative algorithm near the global minimum of Δ (blue). We added spectral noise for I0(f), I1(f), and I2(f) at a level of 2% of the maximum of I0(f). We assumed that B1 is not well known and run the algorithms for various B1 values in the 0.5-2 range. The minima of error functions corresponded to the original B1 values (Figure 7, bottom row).
The first representative example is a pulse at 2.8 μm with the energy of order 1 nJ, an FWHM duration of 340 fs and an FWHM Fourier transform limited duration TFTL = 260 fs (Figure 7, left column). Pulses with similar characteristics can be routinely generated with mode-locked fluoride fiber lasers [16]. We set φreal = −0.4[ps 2 ]f 2 + 0.15[ps 4 ]f 4 , B1 = 0.9, and B1 = 1.8. Such B-integrals can be achieved for a sub-cm fiber length. For example, B2 = 1.8, if the launched energy is only 0.85 nJ and L = 0.5 cm (γ = 150 (W km) −1 @ 2.8 μm). For a sech 2shaped spectrum, T0 ≈ TFTL/1.763 = 150 fs, therefore LD = 9 cm is much longer than 0.5 cm and the dispersionless approximation is correct. Since the proposed algorithm for searching the polynomial phase does not take into account the fourth-order contribution, the mismatch with the original data is visible but not very high. Using the iterative algorithm after the reconstruction of the quadratic phase allows for refining the pulse characteristics (Figure 7, left column). Error Δ was decreased from 0.0059 to 0.0017 (this took eight iterations).
The second example is a pulse at 3.5 μm with an FWHM duration of 4 ps and FWHM Fourier transform limited duration TFTL = 2 ps (T0 ≈ 1.1 ps) (Figure 7, middle column). For such a pulse LD = 2.8 m, therefore, the L maximum can be about 15-20 cm. We set φreal = −150[ps 2 ]f 2 + 450[ps 3 ]f 3 + 0.05π⸱sin(f/0.01[THz]). We simulated that for B2 = 2.2, L = 20 cm, Next, we simulated diverse 3 µm pulses with realistic parameters and reconstructed them using the hybrid algorithm, assuming that an As 2 S 3 fiber with a core diameter of 6.5 µm serves for spectrum conversions. We purposely set spectral phases with leading polynomial contributions due to quadratic and cubic terms and added different regular perturbations. Representative examples are shown in Figure 7. Each column corresponds to an individual pulse with original data (black), data reconstructed with the proposed algorithm for searching an optimal polynomial spectral phase (magenta), and data reconstructed with the iterative algorithm near the global minimum of ∆ (blue). We added spectral noise for I 0 (f ), I 1 (f ), and I 2 (f ) at a level of 2% of the maximum of I 0 (f ). We assumed that B 1 is not well known and run the algorithms for various B 1 values in the 0.5-2 range. The minima of error functions corresponded to the original B 1 values (Figure 7, bottom row).
The first representative example is a pulse at 2.8 µm with the energy of order 1 nJ, an FWHM duration of 340 fs and an FWHM Fourier transform limited duration T FTL = 260 fs (Figure 7, left column). Pulses with similar characteristics can be routinely generated with mode-locked fluoride fiber lasers [16]. We set ϕ real = −0.4[ps 2 ]f 2 + 0.15[ps 4 ]f 4 , B 1 = 0.9, and B 1 = 1.8. Such B-integrals can be achieved for a sub-cm fiber length. For example, B 2 = 1.8, if the launched energy is only 0.85 nJ and L = 0.5 cm (γ = 150 (W km) −1 @ 2.8 µm). For a sech 2 -shaped spectrum, T 0 ≈ T FTL /1.763 = 150 fs, therefore L D = 9 cm is much longer than 0.5 cm and the dispersionless approximation is correct. Since the proposed algorithm for searching the polynomial phase does not take into account the fourth-order contribution, the mismatch with the original data is visible but not very high. Using the iterative algorithm after the reconstruction of the quadratic phase allows for are close to each other). Using the iterative algorithm near the global minima of the error functions allowed for refining the features of the reconstructed signals and obtaining an almost perfect agreement with the original signals (black and blue lines almost coincide in Figure 7). The error function was reduced by an order of magnitude after applying the iterative procedure. Figure 7. Reconstruction of mid-IR 3 μm class laser pulses with realistic parameters for commercially available arsenic sulfide glass fibers as nonlinear elements for spectrum conversion. Each column corresponds to an individual pulse. The upper row shows intensity profiles and phases; the second row shows the initial spectra and spectral phases; the third and fourth rows show converted spectra for different B-integrals. Solid black curves correspond to original signals, magenta curves correspond to pulses obtained by the proposed algorithm for polynomial spectral phase searching (dashed-temporal or spectral intensities, dash-dotted-temporal or spectral phases), and blue curves are obtained by the iterative algorithm near the global minimum of error Δ (dashed-temporal or spectral intensities, dash-dotted-temporal or spectral phases). Bottom row shows error functions calculated for different B1 values (the minima correspond to the original B1 values). Magenta graphs (left axes) are for algorithm optimizing polynomial phase and blue graphs (right axes) are for iterative algorithm.

Discussion
In the presented work, we paid attention to the development of a method of ultrashort pulse reconstruction based on processing spectra converted into fibers with high third-order Kerr nonlinearity. The method without ambiguities requires only three spectra: the initial spectrum of an ultrashort pulse and two converted spectra for the different values of B-integrals [13,14]. It was shown previously that exact knowledge of the B-integral values is not necessary, their ratio is quite enough [13]. Note that for measuring pulse Figure 7. Reconstruction of mid-IR 3 µm class laser pulses with realistic parameters for commercially available arsenic sulfide glass fibers as nonlinear elements for spectrum conversion. Each column corresponds to an individual pulse. The upper row shows intensity profiles and phases; the second row shows the initial spectra and spectral phases; the third and fourth rows show converted spectra for different B-integrals. Solid black curves correspond to original signals, magenta curves correspond to pulses obtained by the proposed algorithm for polynomial spectral phase searching (dashedtemporal or spectral intensities, dash-dotted-temporal or spectral phases), and blue curves are obtained by the iterative algorithm near the global minimum of error ∆ (dashed-temporal or spectral intensities, dash-dotted-temporal or spectral phases). Bottom row shows error functions calculated for different B 1 values (the minima correspond to the original B 1 values). Magenta graphs (left axes) are for algorithm optimizing polynomial phase and blue graphs (right axes) are for iterative algorithm.
The second example is a pulse at 3.5 µm with an FWHM duration of 4 ps and FWHM Fourier transform limited duration T FTL = 2 ps (T 0 ≈ 1.1 ps) (Figure 7, middle column). For such a pulse L D = 2.8 m, therefore, the L maximum can be about 15-20 cm. We set ϕ real = −150[ps 2 ]f 2 + 450[ps 3 ]f 3 + 0.05π·sin(f /0.01[THz]). We simulated that for B 2 = 2.2, L = 20 cm, and γ = 85 (W km) −1 @ 3.5 µm, the peak power is only 130 W, which corresponds to the launched pulse energy of 0.8 nJ (with allowance for the temporal shape of the pulse). Error ∆ was decreased from 0.0151 to 0.0012 after running the iterative algorithm (16 iterations). The final example is a pulse at 3.9 µm with an FWHM duration of 18 ps and T FTL = 10 ps (T 0 ≈ 5.7 ps) (Figure 7, right column). For such a pulse, the dispersion length is very large, L D = 140 m, therefore the maximum fiber length can be about a few meters. We set ϕ real = −900[ps 2 ]f 2 + 600[ps 3 ]f 3 + 0.05π·cos(f /0.005 [THz]) and obtain that for B 2 = 2.6, γ = 60 (W km) −1 @ 3.9 µm, and L = 5 m, the peak power is as low as 9 W and energy can be lower than 0.2 nJ. Error ∆ is decreased from 0.0256 to 0.0030 after running the iterative algorithm (23 iterations).
Thus, for all tested pulses shown in Figure 7, the algorithm for searching the polynomial phase gave quite a good agreement with the original data (black and magenta lines are close to each other). Using the iterative algorithm near the global minima of the error functions allowed for refining the features of the reconstructed signals and obtaining an almost perfect agreement with the original signals (black and blue lines almost coincide in Figure 7). The error function was reduced by an order of magnitude after applying the iterative procedure.

Discussion
In the presented work, we paid attention to the development of a method of ultrashort pulse reconstruction based on processing spectra converted into fibers with high thirdorder Kerr nonlinearity. The method without ambiguities requires only three spectra: the initial spectrum of an ultrashort pulse and two converted spectra for the different values of B-integrals [13,14]. It was shown previously that exact knowledge of the B-integral values is not necessary, their ratio is quite enough [13]. Note that for measuring pulse trains delivered by mode-locked laser systems, one can use only a single piece of a nonlinear fiber. The exact ratio of the B-integrals can be obtained by controllable attenuation of a signal before the fiber by measuring the output average power. For differently attenuated pulses, the ratio of B-integrals is equal to the ratio of average powers.
In this work, we proposed and implemented a novel, very simple and efficient algorithm for data processing which is based on searching for the optimal polynomial spectral phase. Through numerous simulations, we showed that the algorithm is noise-tolerant. Moreover, the results of its execution can be used to seed the iterative Gerchberg-Saxton-like algorithm for refining the "right" solution.
The considered method of spectral phase reconstruction using the proposed and implemented numerical algorithm for searching the polynomial phase is applicable for measuring ultrashort pulses in a wide range of durations. The upper duration limit is determined primarily by spectrometer resolution and peak powers achievable for long pulses (about a few tens ps or may be up to~100 ps). The lower duration limit is determined by dispersion effects in nonlinear fibers, which lead to pulse distortions at relatively short durations (shorter than a few tens of fs). However, earlier we proposed a modified iterative algorithm with allowance for fiber dispersion to mitigate the requirements for a lower duration limit and demonstrated its applicability in experiments [14]. Now let us address the limitations related to a pulse peak power. The minimum value of B-integral at which the method is applicable is about 1. At noticeably smaller values of B-integrals, the pulse spectrum is converted insignificantly, which, taking into account noise and imperfection of the spectral measurements, does not allow for reliable reconstruction of the phase. Therefore, the minimum characterizable signal peak power depends on γ and L. The length of the fiber segment L can be limited by dispersion effects. Therefore, for estimation, we can take L << L D . Thus, the peak power of the measured signals should be P >> 1/(γL D ). Note that for sufficiently long pulses with a duration of several ps or tens of ps with a narrowband spectrum, the limiting factor may be not dispersion but optical losses.
We demonstrated that the method with the developed algorithms can be applied for the reconstruction of pulses in the mid-IR range, namely 3 µm wavelength pulses, using commercially available chalcogenide As 2 S 3 -based glass fibers. Such fibers are transparent in this range and have a huge Kerr nonlinearity. We demonstrated numerically by various representative examples with realistic parameters that the proposed algorithm for searching the polynomial spectral phase gives a fairly good match, and the subsequent application of the iterative algorithm near the global minimum of the error function allows achieving an almost perfect match between the original and reconstructed signals. The next step of the work can be the experimental implementation of the method for measuring mid-IR pulses using chalcogenide fibers. Additionally, note that As 2 S 3 -based glass fibers can be used to characterize ultrashort pulses in a wider spectral range of up to 5-6 µm. The use of other types of commercial chalcogenide fibers, such as As 2 Se 3 -based glass fibers or customized chalcogenide fibers based on glasses with special compositions, can provide even more opportunities for characterization of mid-IR ultrashort optical pulses since they can have higher nonlinear coefficients and a wider transparency range [19][20][21][22][23].

Conclusions
To conclude, we developed and investigated a novel numerical algorithm for ultrashort pulse reconstruction based on processing the initial pulse spectrum and two spectra converted into fibers with high third-order Kerr nonlinearity. This simple, efficient, noisetolerant algorithm assumes spectral phase approximation by a polynomial function and finds the global minimum of error function using an exhaustive search on a reasonable grid of polynomial coefficients. Moreover, the results of the algorithm execution can be used to seed the iterative Gerchberg-Saxton-like algorithm for refining the solution (finding non-polynomial addition to the spectral phase). We demonstrated that the method can be applied for the reconstruction of pulses in the mid-IR range, namely 3 µm wavelength pulses, using commercially available chalcogenide As 2 S 3 -based glass fibers.