Intensity-only-measurement mode decomposition in few-mode fibers

: Recovery of optical phases using direct intensity detection methods is an ill-posed problem and some prior information is required to regularize it. In the case of multi-mode fibers, the known structure of eigenmodes is used to recover optical field and find mode decomposition by measuring intensity distribution. Here we demonstrate numerically and experimentally a mode decomposition technique that outperforms the fastest previously published method in terms of the number of modes while showing the same decomposition speed. This technique improves signal-to-noise ratio by 10 dB for a 3-mode fiber and by 7.5 dB for a 5-mode fiber.


Introduction
Recent renewed interest in multi-mode fibers (MMFs) and, in particular, in few-mode fibers (FMFs) is dominated by their potential in telecommunications, however, their features are equally important for a wide range of applications. In contrast to conventional single-mode fiber (SMF), FMFs and MMFs provide larger mode areas allowing nonlinearity suppression and damage threshold improvement for high-power laser systems and amplifiers [1,2]. Such fibers open up new ways of imaging [3], including starlight [4], sensing [5], and fundamental studies of novel nonlinear dynamics [6] including optical beam self-cleaning [7], spatiotemporal solitons [8] and ultrafast characterization methods [9,10]. Moreover, an unrelenting capacity increase of the contemporary telecommunication systems dictate using novel tools and principles of signal transmission to seize this demand [11]. Spatial mode-division multiplexing is one of the popular approaches to improve the capacity of optical links especially by using FMFs that can significantly increase communication capacity compared to SMF without a drastic complexity of signal processing [12]. However, all of these applications require beam characterization not only in terms of beam divergence (M 2 parameter), but also in terms of waveguide eigenmodes amplitudes and phases. This characterization is also known as mode decomposition (MD) problem. Recovery of a signal from the measured intensity is an attractive approach with various applications ranging from imaging to telecommunications. However, direct detection-based MD methods are usually computationally-intensive and/or capable of recovering few modes only.
Some of the existing MD methods require a reference beam. Digital holography [13,14] and multi-lane light conversion [15] methods are both in this method group. The main drawback of these methods is the requirement of a coherent laser source in the receiver part of the system that limits potential applicability. Those methods that do not require a reference beam can be divided into iterative and non-iterative algorithms. Iterative methods like Gerchberg-Saxton technique [16], line-search [17], and stochastic parallel gradient descent [18] show high accuracy and speed.
However, these methods are vitally sensitive to initial conditions and can be looped at a local minimum.
On the other hand, non-iterative MD methods include machine learning methods [19][20][21][22][23][24] and fractional Fourier system [25]. Machine learning MD algorithms outperform iterative methods in the decomposition speed with a price of high-performance computer requirements, long training times, and poor performance in cases with more than five spatial modes. Recently, a novel non-iterative method, based on dividing inherently non-linear MD problem into a cumbersome linear part and a simple non-linear part, was proposed [26]. This algorithm was based on near-field (NF) MD and allowed to decompose up to 100,000 multi-mode distributions per second for three-, five-, and eight-mode fibers; along with maximum 27-mode decomposition in noiseless case.
Here we propose a significant improvement of the algorithm presented in [26] in terms of the decomposition accuracy by using far-field (FF) and near-far-field (NFF) MD techniques. In addition to the accuracy boost, both techniques allowed a boost in the maximum number of modes decomposed to up to 46 modes for FF method and up to 49 modes for NFF method. In addition, all three possible MD technique variants based on the NF, FF, and NFF decomposition techniques are compared in terms of their accuracy, number of decomposed modes and computational complexity. Moreover, we experimentally verify proposed FF method using a 5-mode fiber and a 21 dB signal-to-noise ratio. We expect the developed technique to find many applications from ultra-fast fiber imaging to laser beam characterization.

Methods
To solve the MD problem, it is required to determine the weighting factors C k = A k exp(iϕ k ) at which a given intensity distribution is observed at the fiber end facet (near-field) or at the infinite distance from the fiber end facet (far-field). The inference of amplitudes and phases based on the intensity distribution is a non-linear problem that can be solved with different approaches. In the previous article, we presented a fast MD method in few-mode fibers, which is based on dividing of the problem into a cumbersome linear part and a simple non-linear part [26]. The presented method makes it possible to determine the coefficients based on the intensity distribution in the near field: where P NF (x, y) is the intensity distribution in the near field, and Ψ k is the near field distribution of k th eigenmode of the fiber. However, the presented method makes it possible to perform MD only up to the complex conjugation of the weight coefficients, since it is based on the analysis of the intensity distribution in the near field. The near-field (NF) technique requires the use of a rather bulky 4-f system ( Fig. 1(a)) to magnify and project the NF distribution onto the image sensor. In this paper, we present a development of the method that allows either the use of a much more compact experimental FF image capture system ( Fig. 1(b)) or the use of both near and far-field (NFF) methods ( Fig. 1(c)) to determine the true values of the mode coefficients, rather than only the complex conjugations.
Far-field measurements are more convenient both because of the compactness of the optical image capturing system and the possibility to easily scale the image resolution. In the case of NF intensity distribution capturing, an increase in the resolution requires an increase in the magnification of the 4-f lens system, which leads to the proportional scaling of the experimental scheme. Thus, to obtain a magnification factor of 67 for the NF capturing scheme, we used a 4-f system with focal lengths of 4.5 mm and 300 mm, which led to an experimental setup overall length of more than 600 mm. To obtain a higher magnification and, accordingly, a larger image size, it is necessary to further increase the focal length of a long-focus lens, which leads to an almost inapplicable length of the experimental setup. In the case of the FF capturing, the magnification coefficient of 67 can be obtained with the length of the free-space beam path of 5 mm, and to obtain a magnification factor of 200, only a length of 15 mm is required. To solve the problem of the FF MD, it is required, as in [26], to calculate the matrix of pairwise products of the fiber eigenmodes. First of all, we calculate the FF distribution for each of the fiber modes and their pairwise products. The FF distributionΦ(ξ, η) can be calculated as the 2D Fourier transform of the NF distribution: We can rewrite the Eq. (1) for the FF distribution: where P(x, y) is the FF intensity distribution. In FF case, either a purely real distribution (for centrally symmetric modes) or a purely imaginary (for centrally antisymmetric modes) is observed in the FF. When calculating the matrix of pairwise products of eigenmodes, we "turn" the phase of the purely imaginary FF modes by an angle of π/2 to bring them to purely real distributions: In what follows, we will take this rotation of π/2 into account when calculating the phases of the eigenmodes. We also need to calculate the matrix of pairwise products of eigenmodes for the FF distribution T FF : Here the upper index (m) corresponds to the m th pixel of the image of eigenmode. Figure 2 shows the distributions of the T FF matrix components for 5-mode fiber. Each image corresponds to one column of matrix T FF . Each image is a point-wise product of the corresponding eigenmodes, and the column of the T FF matrix is obtained by rearranging these M by M distributions into M 2 by 1 vector. The same rearranging is then applied to the intensity distribution to obtain vector P. The dependency of matrix T FF condition number on the number of fiber eigenmodes is shown in Fig. 3.
In comparison to the T matrix for the near field distribution (see [26]) there is no sharp increment in the region around 10 modes. This is due to the fact that unlike the near-field matrix T the columns of the T FF matrix which are composed of LP 11 and LP 12 modes do not match exactly (see Fig. 4). Such a lack of symmetry positively affects the problem numerical solution, as the T FF is well-conditioned. The phase recovery method is generally similar to the technique described in [26]. We first rewrite the Eq. (3) in the linear form using Eq. (5) and solve the linear problem: Here P is the intensity distribution in the FF rearranged into a M 2 by 1 vector, T FF is the matrix of pairwise products of the far-field of eigenmodes.  For combinations of eigenmodes including either only centrally symmetric or only centrally antisymmetric modes, we obtain: or z n = A k A j cos(ϕ k − ϕ j ).
For combinations of eigenmodes including one centrally antisymmetric mode, we get: or z n = A k A j sin(ϕ k − ϕ j ), since we have multiplied the central antisymmetric modes by exp(iπ/2). Including the above mentioned relations, one can represent the vector z as a matrix Z. To do this, we numerate the vector z along the main diagonal first, and then along the columns of the lower triangular matrix Z: And, for example, for 5-mode fiber we can obtain the matrix Z as following: We can also obtain the general relations between the amplitudes, phases and the matrix eigenvalues as following: To determine the sign of the phase when taking the arcsine or arccosine, we use the second column of the matrix Z . It should be noted that we use only the first 3N − 3 components of the vector z for calculating mode weights. For a number of modes N>3, it means that we do not need to calculate the product of full matrix T −1 with the intensity vector P. We need to calculate only the product of the first 3N − 3 rows of this matrix with the intensity vector. Thus, based on the FF intensity distribution, we calculate the vector z, using the Eq. (6) by multiplying the pseudo inverse matrix of pairwise products of eigenmodes by the vector of intensity. Then using this vector z we calculate the amplitudes and phases of the modes using Eqs. (9) and (11).
The analysis of only the NF or only the FF does not allow one to unambiguously reconstruct the phases of the modes, since there is phase transformation that leaves the intensity distribution unchanged. In the case of the near field, this is a complex conjugation of the weights (changing the signs of all phases to opposite ones): ϕ new k = −ϕ k , and for the FF, this is the replacement of the phase sign with the opposite one for centrally symmetric modes (for example, LP 21 or LP 02 ): ϕ new k = −ϕ k and transformation ϕ new k = π − ϕ k for non-centrally antisymmetric (for example LP 11 or LP 31 ): The phase transformations, which conserve the intensity distribution, are different for the NF and FF. This means that knowing both NF and FF distributions allow to calculate the phases exactly, and not only up to complex conjugation (as in the case of the near field) or transformation (12), as in the case of the FF. Therefore, if both FF and NF intensity distributions are measured, then one can decompose mode weights up to a constant phase shift, not only to an intensity-conserving transformation. Figure 5 shows how NF and FF intensity distributions look like for different sets of weights that show either the same NF intensity distribution or the same FF intensity distribution.
A complete MD can be approached either by calculating the cosine value of the second mode by using NF decomposition and then applying FF MD or it can be done by combining NF and FF matrices in one system of equations. The first approach is just a simple extension of the presented FF method and has the same accuracy as the FF MD method. The second approach requires complex-value calculations which increases the complexity of the method, leading to higher time consumption and accuracy. To do this we simply combine both NF and FF systems of equations and just make a change of variables as simple as: In this case, we cannot combine C k C * j and C * k C j as there is no symmetry in matrix T anymore. In the NFF case we use complex-valued pseudo inverse T −1 matrix to calculate the complexvalued z vector. Mode amplitudes and phases are then calculated similarly to (11), but cosine values are proportional to real parts of the corresponding values of matrix Z and sine values are proportional to imaginary parts of the corresponding values of matrix Z.

Results
The MD based on the FF analysis shows better results in the noiseless case compared to the NF case in terms of the decomposition accuracy. In addition, the method works up to 46 modes in the noiseless decomposition (Fig. 6) leading to a significant increase in comparison to 27 mode decomposition for the NF case [26]. For 48-mode fiber noiseless decomposition is performed with some errors, but the most of the weights are still determined with a high accuracy (Fig. 7). So, in principle, the proposed method can be applied for the FF MD in up to 46-mode fiber. The effect of the noise presence on the accuracy of the proposed MD algorithm has also been studied using following error metrics: The amplitude error ϵ A represents the normalized amplitude error, the phase error ϵ φ shows the maximum phase error across phases of all modes normalized by 2π, and the total error ϵ C represents the normalized complex-valued weights error. The latter error includes both amplitude and phase errors. To calculate the accuracy of the presented algorithm we performed multiple decompositions for random sets of mode weights. For FF MD we chose amplitudes from a continuous uniform distribution: A true 1..N U(0, 1). To avoid the ambiguity in the MD we chose the phase of the second mode from a continuous uniform distribution ϕ true 2 U(−π/2, π/2) and all other phase coefficients ϕ true 3..N U(0, 2π). Then we calculated the true FF intensity corresponding to each set of weights. We added Gaussian noise, and used this noisy intensity distribution to recover mode weights using our algorithm. Based on the recovered weights we calculated error metrics. Figure 8 shows examples of noisy MD using our algorithm. We consider a total error ϵ C of less than 0.1 as acceptable. So later we refer to this value when we calculate the SNR value necessary to achieve such an accuracy. However, even if ϵ C >0.1 but less than 0.2-0.3, the recovered weights lay close to their true values. Figure 9 shows how these errors depend on the signal to noise ratio and the number of modes. Each point is averaged over 5000 decompositions. Decomposition is performed by processing 100x100 FF intensity distributions. The method can be applicable for the MD with an accuracy of 0.1 in 3-mode fiber when SNR is better than 15 dB, in 5-mode fiber for the SNR>22.5 dB and in the 6-mode fiber when the SNR is better than 32 dB. For the 8-mode fiber only amplitudes can be determined with an accuracy of 0.1 when the SNR is better than 25 dB, and for the 10-mode fiber the SNR better than 33 dB is required to determine amplitudes with the same accuracy. This results in a superior performance compared to our previously published MD algorithm based on the decomposition of the NF intensity distribution.
To calculate the accuracy of the NFF method we chose amplitudes from A true 1..N U(0, 1) and phases from ϕ true 3..N U(0, 2π) as this method allows decomposing mode weights without any ambiguity. In a noiseless case, the NFF method can perform MD for 49-mode fibers (Fig. 10). This method is also better in the presence of noise. Figure 11 shows how decomposition errors depend on the SNR and the number of modes. This graph was obtained with the same metrics as theFig. 9. Fig. 11. Amplitude, phase, and net errors for NFF decomposition depending on the noise level in the input image and the number of modes. Error metrics are calculated for an image size of 100x100 pixels.
This method can be applicable for the MD with an accuracy of 0.1 for 3-mode fiber when the SNR is better than 10 dB, in 5-mode fiber for the SNR>19 dB, and for the 6-mode fiber when the SNR is better than 31 dB. For 8-mode fiber only amplitudes can be determined with an accuracy of 0.1 when the SNR is better than 25 dB, and for the 10-mode fiber the SNR better than 32 dB is required to determine amplitudes with the same accuracy. Thus, the desired accuracy of 0.1 is achieved with 3 dB lower SNR for 3-and 5-mode fibers when using NFF intensity technique than when using only FF intensity distribution.

Computational complexity
One of the universal metrics of the computational complexity of algorithms is the number of multiplications used. In this section, we calculate the number of multiplications for our algorithms. The MD method based on the FF analysis consists of the multiplication of pseudo inverse of the matrix of pairwise products of eigenmodes with the vector of intensity. For intensity image of size M by M and the number of modes N, one needs to calculate the product of the first 3N − 3 rows of the size of matrix T −1 FF with the intensity vector. So, the linear part of our algorithm requires Q lin = 3M 2 (N − 1) multiplications.
The nonlinear part consists of calculating square roots, arc-cosine, and cosine functions which are not taken into account as these operations can be performed by using lookup tables. The only operation that requires multiplication is for each of the N-1 phase coefficients. Here we have 2 multiplications and 2 divisions, so for the nonlinear part of our algorithm we have Q nln (N, M) = 4(N − 1) multiplications. The total number of multiplications for the FF method is, therefore: As 3M 2 >>4 in every practical case, Q FF ≈ 3M 2 (N − 1) a good approximation.
To compare the computational complexity with the best CNN-based MD algorithm presented before [24] we need to calculate the complexity of the VGG-16-like CNN that was used for MD. The adaptation of the VGG-16 convolutional neural network that was used for the MD problem consists of 13 convolutional layers (CL) and 2 fully connected (FC) layers. The number of convolutional filters and sizes of layers is shown in Table 1. The VGG-16-like CNN consists of N c = 13 convolutional layers. Each convolutional filter is of 3x3 pixels size. So, the number of multiplications per 1 convolution is n = 9. Padding: same was used on each convolutional layer. So, the number of convolutions per layer is the number of pixels on this layer. For a layer of size M by M pixels, it is M 2 convolutions.
The total number of multiplications needed during propagation through N c = 13 convolutional layers is where M i is the size of i-th layer and C i is the number of convolutional filters applied at i-th layer. The total number of multiplications per 1 FC layer is the product of the number of nodes in the previous layer and the number of nodes in the current layer. So, for the first FC layer with F 1 elements, it is defined as The second FC layer is of size F 2 = 2N − 1, so the total number of multiplications for the second layer is Then to recover phases from the cosine values the CNN-based method requires the calculation of intensity distributions for each set of possible phase values and then compare them to the true intensity distribution. For each of the N − 1 cosine values, there are 2 possible phase values that correspond to this cosine value, so the total 2 N−1 intensity calculations are needed. Each intensity calculation consists of (2N − 1)M 2 multiplications as each of M 2 pixels for each of N modes must be multiplied by its amplitude (N multiplications) and phase (N − 1 multiplications) values. So, the total number of multiplications for deriving phase coefficients from their cosine values is This leads to the total number of multiplications per 1 MD using CNN: An approximation of the above equation can be given by The NFF MD method requires multiplications of complex values. Each multiplication of complex-valued numbers requires 4 multiplications of real-value numbers. Moreover, the number of pixels is doubled as we analyze both NF and FF intensity distributions. Therefore, for the complete MD algorithm the number of multiplications is (23) Figure 12 shows how the performance gain depends on the number of modes. The NFF MD algorithm requires 8 times more multiplications than the FF only method, but it benefits from both higher accuracy and ability to recover phase coefficients unambiguously (not only up to transformation (12)) which is the main advantage in comparison to CNN-based MD algorithm that does not allow to completely decompose phases [24]. Thus, the proposed techniques require orders of magnitude less multiplications than the best CNN-based MD algorithm presented to date.

Experimental verification
We experimentally tested the FF algorithm for a 5-mode fiber. We assembled an experimental setup shown in Fig. 1(b). We used a piece of SM2000 fiber and a laser with a central wavelength of 1060 nm. This fiber supports 5 modes at the laser operating wavelength. The corresponding V number of the fiber at the specified wavelength is 3.91. A lens with a focal length of 4.5 mm was used to convert NF distribution at the fiber facet to the FF distribution on the camera. We used a polarizer beyond the lens to select only one polarization. The output intensity distribution after the lens was captured by a CCD camera Ximea mq013rg-e2. The digital decomposition speed was much higher than the speed of the camera (60 frames per second) and to test the accuracy and performance of the algorithm without limitations of camera speed, a series of images was recorded and then processed. The complex mode structure was achieved by bending and twisting the FMF.
During the pre-processing each picture was cropped to size of 400x400 and the average background intensity was subtracted. The background intensity was measured as a captured intensity distribution far outside the main beam region where the intensity of the main beam is guaranteed to be less than one thousandth of the maximum intensity of the main beam. The noise intensity was measured as a variance of the background intensity. The captured intensity distribution, the background intensity distribution and the histogram of the background are shown in Fig. 13. The calculated SNR value was of approximately 21 dB. We believe that two main sources of noise are present in our case: the light scattered by optical elements and the temperature noise of the camera. To find the center of the beam profile we captured 100 images and guessed some initial values of coordinates of the center as well as the scaling factor of the image. Then we calculated eigenmodes corresponding to these guessed parameters, then performed MD, then reconstructed intensity patterns from these recovered mode weights and finally calculated the mean correlation between measured and reconstructed intensity distributions. Then we numerically optimized the 3 free parameters (coordinates of the center and scaling factor) to maximize mean correlation across all these 100 images. The numerical optimization was performed by Nelder-Mead (simplex) method. The optimized parameters were then used for processing all other experimentally captured intensity distributions. Then the MD algorithm was applied to each intensity distribution. The acquired weights were used to reconstruct the intensity using the known FF distributions and finally the accuracy of MD was analyzed by calculating correlation between the measured and reconstructed intensity distributions. Figure 14 shows examples of captured far field distributions, reconstructed intensities, discrepancies, and corresponding correlations.  Figure 15 shows the correlation for the whole series of images captured and processed by FF algorithm. The mean decomposition time was 52 µs. The correlation is above 0.95 for the whole series of captured intensity distributions.

Conclusion
We developed an intensity-only MD technique based on FF and NFF analysis which is superior compared to the fastest previously published MD algorithm in terms of both number of modes in noiseless case and noise sustainability for MD in the presence of noise. Compared to the best previous results, the number of modes in noiseless case increased from 27 to 46 for FF analysis only and 49 in case of using both FF and NF (NFF). Importantly, we demonstrate a substantial improvement of performance in the presence of noise: the SNR needed for MD with a full error of 0.1 decreased from 20 dB to 15 dB (FF) or 10 dB (NFF) for 3-mode fiber; from 30 dB to 22.5 dB (FF) or 19 dB (NFF) for 5-mode fiber. The NFF also benefits from the ability to completely decompose phase weights of higher modes relative to the fundamental mode. And another advantage of the FF technique is the use of simplified optical scheme compared to NF and NFF methods, that does not only decreases the number of required optical components but also increases the setup compactness what is crucial for the in field applications. The rigorous analysis of computational complexity is provided, and it has shown that the presented algorithm requires orders of magnitude less multiplications than a comparable machine learning-based MD technique. The developed progress in intensity-only MD is another step towards fully passive phase-sensitive multi-mode receivers.