Fiber mode excitation using phase-only spatial light modulation : Guideline on free-space path design and lossless optimization

Phase-only spatial light modulators allow to reshape a Gaussian beam by imposing a given phase distribution along the beam cross section. This technique is widely used in the context of mode-division multiplexing to produce, after propagation through a free-space path, the field designed to excite a given fiber mode. In case of orbital angular momentum modes, the target field is approximated as circularly polarized and several complex algorithms have been developed to increase the purity of the obtained modes. Besides their complexity, those algorithms often exploit higher-order diffraction and spatial filtering, hence entailing power loss. In the theoretical work described here, the mode purity is increased in a simple and efficient way by improving the mode approximation adopted to obtain circularly polarized modes and by optimizing two free parameters in the setup, as demonstrated through pertinent simulations.Phase-only spatial light modulators allow to reshape a Gaussian beam by imposing a given phase distribution along the beam cross section. This technique is widely used in the context of mode-division multiplexing to produce, after propagation through a free-space path, the field designed to excite a given fiber mode. In case of orbital angular momentum modes, the target field is approximated as circularly polarized and several complex algorithms have been developed to increase the purity of the obtained modes. Besides their complexity, those algorithms often exploit higher-order diffraction and spatial filtering, hence entailing power loss. In the theoretical work described here, the mode purity is increased in a simple and efficient way by improving the mode approximation adopted to obtain circularly polarized modes and by optimizing two free parameters in the setup, as demonstrated through pertinent simulations.


I. INTRODUCTION
Mode-division multiplexing (MDM) is ideally the source of a limitless increment of capacity in optical communication channels. Relying on the orthogonality of fiber modes, MDM requires each mode involved in the transmission to be excited (and modulated) separately. This is a very challenging task, and several studies have been conducted to accomplish it. For instance, in optical fibers, long-period gratings have been proposed as mode converters, 1 but their employment in an actual MDM setup has not been reported yet. In free space, widely used methods for fiber mode excitation or conversion are based on the manipulation of light by means of phase plates, 2 q-plates, 3 or spatial light modulators (SLMs). [4][5][6][7][8][9] Although SLMs are usually bulky and expensive, hence unlikely to be used in actual optical networks, they represent an effective tool for prototyping and testing MDM-based systems. For this reason, SLMs are considered in this work. The basic diagram for the excitation of one mode only is presented in Fig. 1.
A collimated, linearly polarized Gaussian beam is manipulated by a phase-only SLM. A quarterwave plate transforms the polarization into circular. The resulting beam is finally injected into a fiber by means of a converging lens. The goal is to produce a transverse field at the air-to-fiber interface that matches the one of a selected mode. Despite the wide use of such a setup as elementary building block for SLM-based MDM systems, an accurate and comprehensive description of how to design and simulate it seems to be missing in the literature. Intuitively, in order to calculate the SLM phase mask required for exciting the desired mode, the main step is solving the inverse problem. Namely, the transverse field generated on the SLM plane by inverse-propagating that of a mode needs to be computed. In several works, 4-10 the transverse field of a monochromatic beam propagating through a lens is calculated at the focal plane as the Fourier transform (FT) of the transverse field at the source plane, and hence the inverse propagation is simply given by the inverse FT. The phase of the obtained field is then recognized to be very similar to that of the mode itself in the region where most of the power is located. Consequently, SLM phase masks are generated by encoding the phase of the mode field, conveniently scaled in size. The resulting phase masks are then optimized by means of optimization algorithms. 5,6 However, this method applies if the distance between the source plane and the lens is the focal length of the lens itself and the lens can be modeled as a parabolic phase shift. 11 This condition cannot always be met in practice, for example when the focal length is too small for a phase plate to be placed between the SLM and the lens. Moreover, imposing distances between optical components reduces the degrees of freedom that can be exploited, as demonstrated here, to achieve higher mode purity in the resulting modes without resorting to lossy or more complicated techniques. 9,10 Therefore, a more general approach for calculating SLM phase masks is needed.
Another assumption widely adopted which might not be met in real scenarios is that of low refractive index contrast (RIC) in the fiber. For step-index fibers (SIFs), any approximate full-vectorial mode can be generated by superposing linearly polarized (LP) modes, if the RIC is negligible. 12 Hence, also hybrid modes carrying orbital angular momentum (OAM) can be obtained in the same way. 10,13 This leads to circularly polarized fields and a consequent important simplification in the mode generation system, simplification adopted in many works. 10,14-16 However, the low-RIC approximation is a source of deterministic errors which increase as the RIC becomes larger, the mode order is sufficiently high, or the fiber does not have a step-index profile.
Here, the mathematical model used to describe the evolution of a monochromatic light beam in a basic setup for mode excitation using phase-only SLMs is presented. In the proposed approach, the constraint of a fixed distance between SLM and lens is absent. It is then numerically demonstrated that such a distance can be chosen to achieve satisfactory results in the excited modes without using complicated algorithms or lossy phase masks. Also, a new approximation to achieve circularly polarized modes is introduced, reducing the error between ideal modes and approximate ones in either case of low or high RIC. Derived for SIFs, it is shown that the same approximation method can be extended to any cylindrical waveguide. Simulation results are presented and evaluated in terms purity of the generated modes. The overall analysis aims at leading the reader in the choice of the Gaussian beam radius and, accordingly, the distance between the SLM and lens by means of preliminary simulations. This paper is organized as follows. In Section II, an introduction to the mode excitation problem and the proposed solution are presented. In Section III, the proposed approximation of cylindrical waveguide modes as circularly polarized is derived. The mathematical model of the free-space setup for exciting the derived modes is detailed in Section IV. Simulation results and the method proposed to optimize them are presented in Section V. Section VI is devoted to the conclusion.

A. Problem description
Several techniques can be used to excite fiber modes. The one discussed here is known as end-fire, consisting in illuminating in free space one end face of a waveguide with a monochromatic beam. The fraction of power transferred from the incoming field to the mth guided mode of the waveguide 095111-3 Astorino, Glückstad, and Rottwitt AIP Advances 8, 095111 (2018) is the power coupling coefficient (PCC). For the sake of the present discussion, only the transverse component of the electric field is considered, which implies that the PCC can be calculated in the transverse electromagnetic form as follows. Let E ⊥in identify the electric field in free space and E ⊥m that of the mth guided mode. The subscript "⊥" is used to emphasize that the field is perpendicular to the propagation direction, here denoted by the versorẑ. Assuming that the waveguide end face is also perpendicular toẑ, the PCC is given by where S is the waveguide end face area, n eff,m is the effective refractive index of mode m, the superscript " * " is the complex conjugate operator, and "·" indicates the scalar product. 17 The goal is to selectively excite one mode at a time, which is achieved for mode m if If this condition is satisfied, the orthogonality of modes ensures that T = 0 ∀ m. Hence, with reference to Fig. 1, selectively exciting the mth mode consists in designing the appropriate SLM phase mask and choosing the optimal spacing between optical components, for a given setup, such that (2) is verified.
To accomplish this end, the proposed approach exploits Fourier optics methods. However, such methods apply to scalar (or, equivalently, linearly polarized) electric fields. Thus, each Cartesian component has to be computed separately as scalar. In other words, two SLMs are required, one for each polarization component, and hence the schematic in Fig. 1 should contain two free-space paths that combine before reaching the fiber end. However, this necessity is removed by approximating the modes of interest as circularly polarized (CP), as discussed in detail in Section III B. For them, the transverse electric field xand y-components are represented by the same scalar field, as in the following formula: wherex andŷ denote any pair of directions such thatx ∧ŷ =ẑ ("∧" denoting the cross product), and E m is a scalar electric field. Therefore, the approximate version of E ⊥m , hereafter called E CP ⊥m , is suitable for been generated in free space by using one SLM only and a quarter-wave plate, as detailed in Section IV A 6, and the optimal condition expressed by (2) becomes scalar, namely whereû is a unitary vector oriented along an arbitrarily chosen transverse direction.

B. Proposed solution
Based on the aforementioned assumptions, let E SLM,m be the linearly polarized electric field at the SLM plane that generates exactly E CP ⊥m at the focal plane, hence satisfying (4). In other words, E SLM,m is the scalar field obtained by inverse-propagating E m . In general, E SLM,m does not have a Gaussian amplitude distribution, implying that it is incompatible with the setup in Fig. 1, wherein a Gaussian beam is necessary because the SLM can only manipulate the field phase. However, it is believed that the phase distribution of E SLM,m carries most of the information about E CP ⊥m , in agreement with the argument widely discussed in Oppenheim and Lim's work. 18 Pursuing this idea, the following steps are performed to generate the SLM phase mask: 1) E CP ⊥m is calculated from E ⊥m and projected onto any transverse directionû to compute E m ; 2) E m is inverse-propagated through the lens, up to the position of the SLM (ignoring the quarter-wave plate), and E SLM,m is therefore calculated; 3) The phase distribution of E SLM,m is extracted and converted into the SLM phase mask.
This obviously entails that, when the system is then excited with a Gaussian beam forwardpropagated towards the fiber, the target condition expressed by (4) is never satisfied, as the calculated (optimal) source amplitude is replaced with a Gaussian one. However, as discussed later on, it is possible to reduce such an amplitude mismatch by changing the Gaussian beam radius and the distance between the SLM and the lens. This allows to maximize the PCC for the mth mode only, entailing that the PCC calculated for the undesired modes is minimized.

III. FIBER MODES
The system in Fig. 1 can be employed to excite modes of any cylindrical waveguide, as long as the targeted modes are approximated as CP, as previously mentioned, due to the use one SLM only. In this section, such an approximation is derived starting from modes guided in SIFs. The obtained results are then generalized.

A. Locally elliptically polarized modes
The starting point is the set of modes obtained by analytically solving Maxwell's equations in cylindrical coordinates (r, φ, z) and assuming that the angular dependence of the solution is the complex exponential function exp(iνφ), ν being the radial mode number. For the sake of this work, only hybrid modes are considered, namely those for which ν 0. Such modes are denoted as locally elliptically polarized (LEP), according to the polarization of their transverse field. 19 These modes, written in terms of cylindrical components, are characterized by an overall phase dependence given by exp[i(±ωt ∓ βz + νφ)], where ω is the angular frequency of the electromagnetic wave, β is the mode propagation constant. The phase contribution due to ω and β is not relevant to the following analysis, and hence is omitted for simplicity.
For SIFs, using as a convention the subscripts "0," "co," and "cl" to associate any physical quantity to vacuum, core and cladding, respectively, the set of functions describing the transverse electric field of all LEP modes is where µ 0 is the vacuum magnetic permeability, J ν = J ν (ur) and K ν = K ν ( r) are, respectively, the Bessel function and the modified Bessel function of the νth order, and the prime symbol stands for the total derivative operator. The variables u and are defined as where k co = n co ω √ µ 0 0 and k cl = n cl ω √ µ 0 0 . n co and n cl are refractive indices, 0 is the dielectric permittivity. By indicating with a the fiber core radius, the terms A, B, and D are given by the following equations: whereas C ∈ C is a free, nonzero parameter.

095111-5
Astorino, Glückstad, and Rottwitt AIP Advances 8, 095111 (2018) A and B are expressed in terms of the P-parameter, defined as By definition, in the low-RIC approximation, hybrid modes for which P ≈ −1 are denoted by HE, and those for which P ≈ +1 are called EH. 20 The propagation constant β is the solution of the so-called characteristic equation for SIFs, reported below: 1 uJ where the tilde symbol indicates that the function on which it is placed is calculated at the point r = a.

B. Circularly polarized modes
The transverse field of LEP modes is, by definition, locally elliptically polarized. However, as will be discussed in Section IV C, the complexity of such a field would heavily reflect on an MDM setup, if LEP modes were to be generated using SLMs. A common way to reduce such a complexity is considering approximate fields that lead to CP modes. Under the assumption of low RIC, this can be achieved by approximating hybrid modes as a superposition of degenerate LP modes. 12 Here, a more accurate method is adopted, with the twofold idea of shifting away from the low-RIC assumption and achieving a more accurate approximation in any case.
By defining ν s as the sign of ν and Γ = ω µ 0 C/K ν , the transverse field in (5) can be written, if projected onto Cartesian coordinates, as follows: Equation (10) shows that the elliptical polarization of hybrid modes can be split into two counterrotating, circularly polarized fields, with angular dependences exp[i(ν − 1)φ] and exp[i(ν + 1)φ].

Proposed approximation for step-index fibers
For brevity, the derivation of the proposed CP modes is conducted for ν s = 1 only, meaning that only clockwise rotating (over time and as a whole) modes are considered. Relevant results are extended to the case ν s = −1 by logical considerations. This method is referred to as circularly-polarized-mode approximation (CPMA).
The idea is to isolate only the "predominant" circular polarization in (10) by conveniently choosing, instead of calculating, the value of P. In case of low RIC and low-order modes, it is known that P ≈ ±1, whereas this is not true in general. The CPMA consists in calculating the exact value of propagation constant for a given full-vectorial mode and setting P equal to +1, when positive, or −1, when negative, even in cases when this approximation is evidently too rough. Using this method, hence substituting P = −1 in (10), the equation for the resulting HE ν ,µ modes is Astorino, Glückstad, and Rottwitt AIP Advances 8, 095111 (2018) where the superscript "HE+" has been introduced to emphasize the validity of (11) only under the CPMA and positive values of ν. Similarly, for EH ν ,µ modes, (10) transformes into

Comparison with the standard approximation
Equations (11) and (12) seem equivalent to those proposed by Gloge 12 and a similar derivation of them can be found in Okamoto's book. 21 The difference is that the CPMA does not use the approximation k co ≈ k cl in (9). Therefore, while nearly degenerate modes have same propagation constant and same radial dependence in the weakly guided approximation, in the CPMA they do not. For instance, using the low-RIC approximation, the two polarization components of HE ν +2,µ modes are indistinguishable from those of EH ν ,µ modes and those modes are all classified as OAM ν +1 . This confusion is avoided here by using the superscripts EH± and HE± and by denoting those modes as CP, instead of OAM.
The proposed approximate modes describe the exact ones with higher accuracy, compared to those obtained using the weakly guided approximation, especially when the basic assumptions of low RIC and low mode order do not hold. This has been verified by calculating the overlap integral (OI) between approximate and exact modes with both methods and several fiber parameters. No exception to the hypothesized improvement has been found, although a formal proof is missing.
For a qualitative explanation, since the two CP fields composing (10) have opposite rotation directions (and different angular dependence), they are orthogonal and cannot be simultaneously reproduced in free space through the setup schematized in Fig. 1. Consequently, if one polarization is selected, the highest PCC is obtained when the field with the chosen polarization is matched by E ⊥in , whereas the remaining component of the field can be completely ignored. This is essentially what is achieved by the CPMA when forcing the value of P to be equal to ±1.
As an example of the validity of the CPMA, a waveguide with parameters n co = 1.45, n cl = 1.00 (i.e., air cladding), and a = 5 µm, at vacuum wavelength λ 0 = 1550 nm, is considered. The OI is calculated between approximate modes and actual ones for all the 111 clockwise rotating LEP modes guided by the fiber under test. The quantity ∆ OI is defined as the percentage difference between the OIs calculated when the CPMA and the standard approximation are adopted. The result is shown in Fig. 2 where, for brevity, only the cases for which ∆ OI > 2% are reported.
FIG. 2. Percentage difference between the overlap integrals calculated when the proposed and low-RIC approximations are adopted.

Generalization of the proposed approximation
The question is how to extend the CPMA to modes of any cylindrical waveguide. In fact, when the analytical solution of the wave equation in the structure is unknown, it is not possible to apply an approximation analog to |P| = 1. This problem can be easily solved by noticing that [see (10) and (11)] where j ∈ {co, cl}. Hence, (13) can be used to approximate any clockwise rotating hybrid mode also when the mode field is computed numerically. It can be shown (the proof is omitted for brevity) that counterclockwise rotating modes obey the same approximation rule, provided that i is negated, yielding As a geometrical interpretation, (13) and (14) transform the ellipses described at any point by the transverse electric field of the actual hybrid modes into circles with diameters equal to the average of the axes of the corresponding ellipses.

A. Forward propagation
The setup in Fig. 1 is translated into the mathematical model shown in Fig. 3. Circled crosses are used to indicate the multiplication operator, whereas rectangular blocks are filters whose spatial frequency response is H(·). Following the diagram from left to right, the field evolution is analyzed below.

Gaussian source
Omitting the time dependence exp(iωt) and the field polarization, the propagating electric field is denoted by U(x, y, z), i.e. a scalar function of the transverse coordinates, traveling along the zaxis with a small diffraction angle. As a convention, if the medium has a discontinuity located at z =z, U(x, y,z − ) and U(x, y,z + ) indicate the limits of the field as z approachesz in the propagation direction and in the opposite direction, respectively. At the source z-coordinate, the field is assumed to have a Gaussian intensity profile and a perfectly flat phase front. Since the beam width is of the order of the SLM screen size, normally much larger than the field wavelength, the beam spreading over the distance between the source and the SLM is neglected. The only change along the propagation direction is therefore approximated to a mere phase shift, which can be incorporated in the field complex amplitude, or simply ignored. Hence, the light beam impinging on the SLM can be described by where U s is a complex constant, s defines the beam radius, and z s denotes the z-coordinate of the SLM screen, assumed to be perpendicular to the z-axis. The subscript "s" stands for "source." FIG. 3. Mathematical model of the setup in Fig. 1.

SLM
The beam impinging on the SLM is reflected, and hence it changes direction. However, to simplify the analysis, the SLM is modeled as an ideal, infinitesimally thin film, whose sole effect is imposing a local phase shift, here denoted by φ SLM (x, y). Hence, the propagation direction is unperturbed by the SLM. Despite that, to match the correspondence with real cases, the light manipulated by this device is often referred to as "reflected." After the influence of the SLM, the field is expressed as According to the analysis conducted in Section II, the amplitude of U(x, y, z + s ) is still Gaussian, but its phase distribution, φ SLM (x, y), is that of E SLM,m .

Free-space propagation (1)
The distorted phase front of the still Gaussian (in the intensity) beam entails a diffractive effect due to the propagation in free space from z s to the lens plane, here indicated by z (where the subscript stands for "lens"). This effect is calculated as follows: where F{·} and F −1 {·} denote the 2D FT operator and its inverse, respectively, calculated in the spatial frequency coordinate system ( f x , f y ), and d = |z − z s |. Therefore, it is possible to define the free-space transfer function, corresponding to a propagation over a distance ∆z, as Hence, with reference to the main schematic in Fig. 3, d is used in place of ∆z to indicate the propagation between SLM and lens.

Lensing
In analogy with the SLM model, the mathematical description of the lens consists in a local phase shift on the transverse plane. This phase shift is denoted with φ lens (x, y), to get U(x, y, z + ) = U(x, y, z − )e iφ lens (x,y) .

Free-space propagation (2)
After the lens, the field propagates again in free space for a distance equal to the focal length, denoted d f (where the subscript "f " denotes "focus"). A different distance might be chosen, although this seldom happens in practice, as it goes against the purpose of lensing. Therefore, selecting the focal length also imposes the distance between the lens and the fiber end face in the model.
At the focal plane z = z f such that d f = |z f − z |, the field U(x, y, z − f ) can be calculated using (17) applied to U(x, y, z + ) and where the distance is d f , yielding In Fig. 3, the corresponding block is labeled H(d f ).

Quarter-wave plate
The quarter-wave plate is physically placed anywhere between the SLM and the lens. For mathematical purposes, instead, its effect is calculated at the focal plane, where it is necessary to leave the scalar representation and express the electric field as 45 • -polarized, namely U(x, y, z − f )(x +ŷ)/ √ 2 or, using Jones's formalism, 22 as Astorino, Glückstad, and Rottwitt AIP Advances 8, 095111 (2018) where the factor 1/ √ 2 has been absorbed by the amplitude of U(x, y, z − f ) for simplicity. Neglecting an overall phase shift, the matrix M ± representing a quarter-wave plate is expressed as where the sign of the imaginary unit, and hence the superscript of M, depends on the orientation of the fast axis of the wave plate: positive, if horizontal, or negative, if vertical. The transverse electric field at the fiber end face becomes The matrices M + and M − generate clockwise and counterclockwise circular polarizations, respectively, and their choice has a fundamental impact on the result. For instance, if U(x, y, z − f ) = E HE+ x,j , a wave plate corresponding to M + would generate which is equivalent to (11a), and hence correctly reproducing the full (approximate) field of any HE ν ,µ mode. On the basis of a similar analysis, if U(x, y, z − f ) = E EH+ x,j , it is necessary to use M − to obtain E EH+ y,j along the y-direction. This is clearly shown by (11) and (12). Notice that it is not possible to generate an EH ν ,µ mode using M + due to the use of the CPMA.
, which has the form of (11) for CP HE ν +2,µ modes, instead of (12). This would cause an error due to the use of a β, and hence a radial field amplitude, belonging to an EH mode to produce a HE mode. Consequently, it must be true that, in case ν < 0, M + and M − generate, respectively, EH and HE modes (a formal proof is omitted for brevity). The allowed combinations are summarized in Table I. The table does not include an exception, represented by the HE ±1,µ modes. Their field, projected on Cartesian axes, do not have angular dependence and therefore both M + and M − are capable to produce the full mode starting from one polarization only. In other words, the two polarization components of HE ±1,µ modes cannot be combined such that the result is an EH mode.

B. Inverse propagation
As discussed in Section II, the proposed method requires the calculation of the SLM phase mask as the phase of the field inverse propagated from the fiber end to the SLM plane. This is achieved by following the diagram in Fig. 4.
The input is the transverse electric field of a CP mode, from which one polarization component is extracted by performing the scalar product with the unit vectorû, arbitrarily oriented on the transverse plane [see (4)]. The result is inverse propagated through a distance d f , the lens, and a distance d . The calculated field carries the spatial phase that can be converted into the SLM image.

C. Two-path vectorial model
For completeness, a model that does not require the CPMA and is able to excite any waveguide mode (including transverse electric and transverse magnetic modes) is briefly discussed. Such a model needs the use of two free-space paths, one for each polarization component of the field, which have to merge into one before reaching the fiber end face. As a mathematical description, the schematic depicted in Fig. 5 is considered.
Notice the presence of the subscripts "x" and "y", necessary to distinguish the two polarization components, generated separately, and the two corresponding phase masks. Also, before combining the two scalar fields at the focal plane, it is necessary to assign a direction to them, achieved with the multiplication by their pertinent versors.
Although this model enables to describe an ideal scenario with higher accuracy, compared to that requiring the use of the CPMA, and allows to generate also TE and TM modes, its implementation in experiments is extremely challenging. In fact, such a model requires the two light beams to have a constant relative phase shift and same intensity when they reach the fiber end face. This means that, due to the finite line width of laser light, the two beams must be generated by the same source (e.g., by using a polarization dependent beam splitter) and experience a controllable, fixed attenuation along their paths. In addition, the distances between each SLM and the lens have to be identical and remain unchanged over time. Any mechanical vibration would alter these distances, hence making the relative phase between the two light paths unstable and, consequently, producing random distortions in the final mode. For these reasons, the two-path vectorial model mostly remains conceptual.

V. SIMULATION RESULTS AND PARAMETER OPTIMIZATION
A fiber with refractive indices n co = 1.450 and n cl = 1.444, corresponding to a numerical aperture (NA) of about 0.4, is chosen. The wavelength in use is 1550 nm, as an example. In order to have a sufficient number of modes for testing the model, a core radius a = 15 µm is selected, obtaining fifteen LEP modes with positive ν, listed in Table II, along with the assigned values of m and their effective refractive indices. Negative values of ν are not considered as they would give exactly the same results as the positive ones. A spheric lens with focal length of 6 mm is considered.

A. Preliminary test
As a test, d = 400 mm, d f = 6 mm, s = 3 mm, and m = 10 (i.e., the HE 3,2 mode) are selected. The input source is calculated by inverse propagating, from z f to z + s , the x-component of E CP ⊥10 , meaning that, in (4),û is chosen equal tox. The calculated E SLM,10 , i.e. the ideal field at the SLM plane that generates the HE 3,2 mode in the CPMA, is depicted in Fig. 6. As per Section II B, the phase distribution of the obtained E SLM,10 (Fig. 6b) is converted into the SLM phase mask, whereas the amplitude profile (Fig. 6a) is discarded and replaced with the Gaussian one given by (15).
Following the diagram in Fig. 3, the resulting scalar field is then propagated through the system and E ⊥in is calculated.
To evaluate the quality of the obtained field at the fiber end face, T 10 is calculated using (1), where E ⊥in is the obtained field and E ⊥10 is the transverse field of the exact EH 3,2 mode. In this preliminary test, T 10 ≈ 0.444. For a graphical comparison between the obtained field and the desired  one, the real part of E 10 and the x-component of E ⊥in are depicted in Fig. 7. The field imaginary part is not shown for brevity, as it is simply a rotated version of the real part.

B. Parameter optimization
It is clear that the result obtained so far is not satisfactory. The low PCC previously calculated entails that part of the field at the fiber end face is converted into other modes. The remainder, instead, spreads in the cladding and is absorbed by the fiber coating. If the coating is not sufficiently lossy, the non-guided field might be reflected by an external layer, enter the core, and couple to undesired modes. These unwanted effects must be prevented. Thus, as a primary goal, it is necessary to improve the coupling efficiency between the target fields and the produced ones.
Clearly, the CPMA and the Gaussian profile of the source amplitude cause errors that negatively affect the PCCs, even in an ideal scenario. However, the low NA adopted in the considered fiber ensures that the mode approximation has a very small influence on the PCC, even smaller than the power loss due the reflection at the air-to-fiber interface, which is about 3%. Therefore, the low PCC obtained in the previous simulation must be ascribed to the mismatch between the amplitudes of E SLM,10 and U(x, y, z + s ), due to the use of a Gaussian source. In order to improve the obtained results without involving advanced techniques for the generation of the phase mask, 5,6,9 it is possible to exploit the existence of some degrees of freedom. For instance, by plotting T m as a function of s for all the guided modes in the considered fiber, it can be observed that the beam size largely affects the PCC, hence indicating that a convenient choice of s must be made (see Fig. 8).
Similarly, also the distance d between the SLM and the lens affects the results. The proposed approach allows such a distance to be arbitrarily changed, as long as the reflected beam does not exceed the lens size and sufficient space is available to host a quarter-wave plate within d . However, any variation of d implies that the optimal s must be recalculated. Consequently, the optimization problem becomes two-dimensional. As an example, Fig. 9 shows T 10 as a function of d and s . For the HE 3,2 mode, the optimal point is found at d ≈ 0.25 m and ≈ 1.11 mm, yielding a PPC of about 0.852. The new E SLM,10 is displayed in Fig. 10. As before, the amplitude of E SLM,10 is replaced by a Gaussian distribution, whereas the phase of E SLM,10 is converted into the SLM phase mask. The real part of the x-component of the optimized E ⊥in is depicted in Fig. 11, along with the real part of E 10 .

C. Optimized results for all guided modes
The function plotted in Fig. 9 is concave, and hence there exists a point such that the T 10 is maximized. In addition, the magnitude of the gradient of T 10 around its maximum point is sufficiently small to allow fine adjustments and a reasonable margin of error in real setups.
The low steepness and the concavity of the PPC found for the 10th mode are properties valid for all the guided modes. Such properties facilitates the search for the optimal PCC as they enable the use of optimization algorithms such as gradient descent, adopted in this work to find the optimal T m for all modes in Table II. A comparison between the PCCs obtained by using the initial values of d and s and those obtained by using the optimal values of these two parameters is provided in Fig. 12. The achieved improvement is evident. However, for a complete evaluation of the final result, it is necessary to consider also the amount of power coupled into unwanted modes. To this end, let X ,m be the normalized cross-PCC (X-PCC) given by where the th mode is undesired and the mth mode is the target one, for which the field E ⊥in is generated. The result obtained after applying the proposed optimization method is summarized in Table III, where only the values of X ,m < 40 dB are reported for brevity. Out of 210 combinations, only in 10 cases the difference between the power coupled into the desired modes and that coupled into an unwanted mode has been lower than 40 dB. As a worst case, a difference of 17.6 dB has been calculated.
It is important to highlight that those reported here are purely theoretical results. Achieving them in actual setups requires s and d to be accurately measurable and adjustable. Moreover, the beam size of the source should be independently set for each mode. In general, this is not a realistic scenario. However, optimizing the simulated setup using the proposed approach can guide experimentalists through the calculation of the beam radius that allows a satisfactory purity in a subset of modes. This subset can therefore be chosen as the target one in the MDM system. In addition, unsatisfactory PPCs can be further improved by using more advanced techniques, such as simulated annealing. 23 Nevertheless, as in any checkerboard-based algorithm, such an improvement requires the undesired field to be scattered into higher diffraction orders, hence introducing losses. 24 Since the proposed  approach is almost lossless, its usage as a first optimization method finds validity especially when power light utilization efficiency is to be maximized.

VI. CONCLUSION
In this paper, the basic setup for exciting modes in any cylindrical waveguide using spatial light modulation has been studied. In particular, a method to generate the local phase shift to apply to a Gaussian beam, used as a light source, has been provided. A mathematical model of such a setup, based on Fourier optics methods, has been illustrated.
Theory has been developed around the assumption that circularly polarized modes are to be excited. Such modes do not exist in cylindrical waveguides, as they are only a well-established approximation adopted to simplify the setup, approximation based on the assumption of low refractive index contrast in the waveguide. With the primary aim of removing this assumption, a method to achieve a more accurate approximation for circularly polarized modes, compared to the conventional one, has also been proposed. Its effectiveness has been discussed and evaluated.
A software implementation of the developed model setup has been used to simulate the excitation of a set of guided modes in a multimode step-index fiber, chosen as an example. The results have been reported in terms of power coupling coefficients between desired modes and corresponding obtained fields, calculated at the air-to-fiber interface. Moreover, the dependence of the power coupling coefficient on two free parameters in the setup has been analyzed. It has been shown that such a dependence is a concave function, hence ensuring the existence of an optimal point, calculated for each mode guided by the considered fiber. By exploiting this finding, a theoretical high purity in the excited modes has been achieved, hence suggesting that the proposed optimization algorithm and phase mask calculation can be used as a guideline when designing a real setups.

ACKNOWLEDGMENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 642355.