Structurally Invariant Higher-Order Ince-Gaussian Beams and Their Expansions into Hermite-Gaussian or Laguerre-Gaussian Beams

: Paraxial beam modes, which propagate in space and focus without changing their transverse intensity pattern, are of great value for multiplexing transmitted data in optical communications, both in waveguides and in free space. The best-known paraxial modes are the Her-mite-Gaussian and Laguerre-Gaussian beams. Here, we derive explicit analytical expressions for Ince-Gaussian (IG) beams for several first values of the indices p = 3, 4, 5, and 6. In total, we obtain expressions for the amplitudes of 24 IG beams. These formulae are written as superpositions of the Laguerre-Gaussian (LG) or Hermite-Gaussian (HG) beams, with the superposition coefficients explicitly depending on the ellipticity parameter. Due to simultaneous representation of the IG modes via the LG and HG modes, it is easy to obtain the IG modes in the limiting cases wherein the ellipticity parameter is zero or approaches infinity. The explicit dependence of the obtained expressions for the IG modes on the ellipticity parameter makes it possible to change the intensity pattern at the beam cross-section by continuously varying the parameter values. For the first time, the intensity distributions of the IG beams are obtained for negative values of the ellipticity parameter. The obtained expressions could facilitate a theoretical analysis of properties of the IG modes and could find practical applications in the numerical simulation or generation of such beams with a liquid-crystal spatial light modulator.


Introduction
When light is used for steering micromachine components and generating multiple optical traps for microscopic particles, laser beam modes are a convenient tool, as they conserve their structure both upon free-space propagation and at the focus of a spherical lens.Another benefit is that mode symmetry can be adjusted by varying certain parameters and passing from symmetry in the Cartesian coordinates to circular symmetry via elliptical symmetry.Ince-Gaussian (IG) modes can be utilized for this purpose, as they are transformed into either Hermite-Gaussian (HG) beams or Laguerre-Gaussian (LG) beams when the ellipticity parameter is changed.
An initial solution to the Whittaker equation in the form of IG functions was obtained by F.M. Arscott in [1,2].In W. Miller's book [3], this solution was obtained in separated variables in elliptical coordinates.The applications of IG beams in optics research were considered in [4].In [5], expressions for the amplitudes of IG beams were obtained for p = 0, 1, and 2, independent of the ellipticity parameter.In [6], elegant IG beams were investigated.These works have created an impetus for the broad adoption of IG beams, along with other well-known laser modes (LG and HG), in optics.In [7], IG beams were generated with a digital hologram.In [8], vector (classically entangled) IG beams were generated as superpositions of a right-handed, circularly polarized, even IG beam and a left-handed, circularly polarized, odd IG beam.IG beams with quantum entanglement were generated in [9].Nonlinear transformation of IG beams by spontaneous parametric down-conversion was investigated in [10,11].In [12], IG beams were generated also parametrically, but without conversion.In [13], elegant IG beams were studied in a parabolic medium.Propagation of IG beams in uniaxial crystals was studied in [14].IG beams are used for underwater data transmission [15].Propagation of IG beams in a turbulent atmosphere was considered in [16].
As seen from the above brief review of works dealing with IG beams, they have been actively studied in optics.However, no analytical representation of these beams via LG or HG modes with an explicit dependence on the ellipticity parameter has yet been proposed.In this work, in an effort to fill in this gap, we derive a number of particular analytical formulae for even and odd IG modes with the indices p = 3, 4, 5, and 6, which are expressed via LG and HG modes and exhibit explicit dependence on the ellipticity parameter ε.This explicit dependence of the expansion coefficients on the ellipticity parameter makes it possible to control the intensity pattern of IG beams by continuously varying this parameter.In this work, for the first time in the context of optics, we consider a situation wherein the ellipticity parameter ε can be not only positive, but negative as well.We analyze how the IG beam changes with changes in the sign of this parameter.The simultaneous representation of IG modes via several LG modes and several HG modes, as proposed herein, suggests that when the ellipticity parameter ε approaches zero, any IG mode is converted into a specific LG mode (more exactly, to its real or imaginary part), and when the ellipticity parameter ε approaches infinity, the IG mode aligns with a specific HG mode.

Solution of the Paraxial Equation in Elliptic Coordinates
The paraxial Helmholtz equation is given by [3]: where ( , , ) x y z are the Cartesian coordinates; ( , ) rr  is a 2D vec- tor; z is the coordinate along the optical axis; k is the wavenumber of light; and ( , ) Ez r is the complex amplitude of a monochromatic light field.Below, we use the following dimensionless coordinates: ( , ) 0 Its simplest solution with a finite energy is a conventional Gaussian beam: It is known that separation of variables in Equation (2) into Cartesian and polar coordinates makes it possible to obtain two solution families-HG and LG beams, respectively [17]: with , 0,1, 2, nm= K .Both families are examples of structurally stable light fields, i.e., at any propagation distance z, the transverse intensity pattern of such fields aligns-up to scale-with that in the initial plane 0 z = :  (5) Thus, it will suffice to study both families at 0 z = , omitting the last argument for brevity: ,, ( n m n m  = rr .We also introduce normalized variants of each solutions family: ) When Equation ( 2) is solved, the variables are separated in the elliptic coordinates ( , )  [3]; then, a family of IG beams, which are also structurally stable, is obtained: ). 1 Here, and the parameter  is supposed to be arbitrary.Substitution of Equation (7) into Equation (2) yields equations for the functions () P  and () Q  : Both equations (9) reduce to the Ince equation, with its canonical form being given by [1-3]: where p is an integer nonnegative number;  is a positive ellipticity parameter, as men- tioned in the Introduction; and ( , ) p

  
= is some constant.In works [1][2][3], it is supposed that ε > 0. This assumption can also be seen from comparison of Equations ( 9) and (10), as ε = 2σ 2 > 0. Below, in the current work, we show that ε in the solution (10) can be of any sign.
The constant λ in Equation ( 9) appears on separating the variables  and  in Equation (2).For Equation (10), the constant is a real number, introduced for the existence of solution () Nt in the form of a trigonometric polynomial of degree p (called the Ince polynomial).Once the form of the solution is restricted in this way, we find that  should be a root of a certain polynomial of degree (p + 1), whose coefficients depend on  (i.e., root of a characteristic equation).Thus, in total, (p + 1) roots appear as 01 , , , p    K , and they are all real-valued.In addition, if these roots are sorted in ascending order at 0  = , their ascending order does not change with increasing values of ε (for all positive ε, the characteristic equation has no multiple roots).Thus, by choosing a root ( , ) q p  and by constructing a trigonometric Ince polynomial can be written via the solutions to Equation (10) as follows: ).Thus, even and odd IG modes (11) are usually written as where () ( , ) Ct  and () ( , ) St  are even and odd Ince polynomials, respectively, and in- dices p and q are nonnegative integers ( qp  ).Unfortunately, a tradition of writing the even and odd IG modes differently and studying them separately leads to unjustified complications in the designations.The index q is chosen to be of the same parity as p, i.e.,

 
, 2, 4, q p p p = − − K .If p is odd, then this list stops at 1 q = , but if p is even., it stops at 0 q = for modes IG e and at 2 q = for modes IG o .This separation is probably related to the initial form of the Ince poly- nomials () Nt , as expressed via sin kt or cos kt .If, instead of the sines and cosines, exponentials ikt e were used and if the even and odd IG modes were not separated, then, for any fixed p, the index q would range from 0 to p, and respective numbers ( , ) q p  would be sorted in ascending order.This change would make it possible to investigate the IG modes in a unified way, similarly to the Hermite-Gaussian modes; despite the Hermite polynomials being both even and odd, the family of the HG modes is not split into even and odd subfamilies.However, at this point, we will adhere to the commonly used designations.
As the IG modes are structurally stable, below, we write them at 0 z = , omitting the argument z for brevity, but specifying explicitly the dependence on the parameter  : e o e o e o p q p q p q z  →→ IG r IG r IG r .Such notation is convenient when studying limit- ing cases 0  = and  =, when the IG modes reduce to the LG and HG modes, re- spectively.In addition, similarly to Equation (6), we use the IG modes with and without normalization: , ( , ) preferring the non-normalized variant, ( , )   , ( , ) eo pq  r , in cumbersome cases, as for the IG modes, it is generally much simpler than the normalized one.

Expansions of the IG Modes by the HG and LG Modes at Small Values of p
The simplest case is when 0 p = .Then, 0 q = and the even IG mode reduces to the conventional Gaussian beam: The following three modes are also independent of the parameter ε and were given in [4,5]: ).

Cases with the Characteristic Polynomial of Degree 2
Below, without derivation, we give explicit analytical expressions for those IG modes, for which the characteristic polynomial reduces to a quadratic polynomial and the respective roots ( , )  can be expressed in square radicals.
It is worth noting that if the initial definition is dropped, the parameter  (ellipticity) can be considered to be a negative number (as observed by Edward Ince [18]).Therefore, in addition to the limiting cases of IG modes with 0  = and  =, we also consider a case wherein  = − .
In Table 1, the following designations are used: , () nm HG r and , () LG r are the normalized Hermite-Gaussian and Laguerre-Gaussian modes (6), and Table 1.Normalized IG beams expressed as superpositions of two HG beams or two LG beams and their limiting cases.

HG r HG r IG r
LG r LG r 2,0 ()

HG r HG r IG r
LG r LG r 3,0 ()

HG r HG r IG r
LG r LG r 2,1 () )

HG r HG r IG r
LG r LG r 0,3 () Im ( ) LG r LG r 3,1 () Im ( ) The coefficients , , , a b c d here are positive, monotonically increasing functions of  , as shown below: We should note a similarity among the IG modes with equal indices p and parity flags ( , ) eo .For instance, to obtain mode , it is suffi- cient to swap the superposition coefficients and change the sign of one of them.It turns out that a similar property is valid in a more general case.Namely, if the whole set of IG modes except one is known for fixed values p and ( , ) eo ,, , this unknown mode can be obtained by manipulating the expansion coefficients of the modes from the set.When the characteristic polynomial is quadratic, there are only two modes and the "manipulation" is just a swapping of the coefficients with a sign change.For more complicated cases, manipulations are also more complicated.However, they are incomparably simpler than obtaining the dependence of the superposition coefficients on the ellipticity parameter.

Cases with the Characteristic Polynomial of Degree 3
For mode 4,2 ( , ) , we obtain the following characteristic equation: whose roots ()  can be found by using the Cardano formula.In ascending order, they read as follows: Asymptotic expansions of all three roots at small and large values of  are given in Table 2. ), We do not present the normalized version in order to avoid cumbersome fractions, and we give the normalized IG modes only in the limiting cases (Table 3).

+HG r
For numerical computation of the IG modes, when the ellipticity parameter is neither very small nor very large, the Cardano formula with expansions (19) should be employed, with the normalizing multiplier added if necessary.If, however, the ellipticity parameter is close to the limiting values, then the asymptotic expressions of the characteristic roots are more convenient.For instance, if 2 m = and 0  → , we ob- tain 24 1 0 144 () and, thus, For many values of  , the roots of the Equation ( 17) can be easily found without using the Cardano formula; for example, this is the case for   .It is enough to choose some integer  as a root and find which  it corresponds to (relative to  , Equation ( 17) is quadratic).In par- ticular, for 2  = , we obtain .
In addition, it can be shown that a simple interrelation exists between the IG modes constructed for the parameters ε and −ε: The first two formulae here are written for the even value of p, whereas the last is written for an odd value of p.This result demonstrates once again that the case of negative ε is no less important in investigating the IG modes than the case of positive ε.In particular, Equation (21) indicates that the IG modes with negative values of the parameter ε are also orthogonal to each other, similarly to the modes with positive values of ε.
The proof of the relationships ( 21) is based upon the following property of the characteristic polynomials: if is e ( ven, is odd. , ) ), if Now, we consider the case wherein 5 p = .For the even IG modes, the characteristic equation is given by We again apply the variant of the Cardano formula wherein the cubic equation has three real-valued roots, and these roots in the ascending order read as follows: ( ) Asymptotic expansions of all three roots at small and large values of  are given in Table 4.For the limiting cases, the normalized IG modes are shown in Table 5.
Similarly to the above considered case for 4 p = , for many values of ε, the roots of the characteristic polynomials can be found without using the Cardano formula.For instance, for the even IG modes, the following values can be chosen: ) .There is one more case with the cubic characteristic equation, the case 6 p = for the odd IG modes:

Even IG Modes at p = 6
Now we consider the case p = 6 for the even IG modes.The characteristic polynomial is quartic: Its roots ()  can be obtained by the Ferrari formula.We use the variant of this for- mula presented in [19].Let us introduce an auxiliary variable Asymptotic expansions of all roots at small and large values of  are written in Table 6.

A A A A A A A A A A A
For the limiting cases, the normalized IG modes are shown in Table 7.
Table 7. Limiting cases of normalized even IG modes at p = 6.

Numerical Methods and IG Modes
It is known that for all characteristic polynomials ( , ) ( , )  by powers of ε (both at 0  → and at  →  ) have expansion coefficients that are rational numbers (see the above formulae as examples of such expansions).Consequently, asymptotic expansions of the IG modes by the LG and HG modes have coefficients that can be represented as series of the powers of ε (if 0  → , then by powers ε, ε 2 , ...; if  →  , then by the neg- ative powers), and the coefficients in these series are also rational numbers.
For instance, let us consider mode In particular, as , then, as a limiting case, we obtain an identity that is already well known: LG r .
If  → + then Here, in the limiting case, we obtain a HG mode: = − IG r HG r .

Applying the Padé Approximants for Approximate Computation of IG Modes
In this section, we demonstrate using a concrete example how the Padé approximants can be employed to compute the IG modes.In the limiting case, IG mode reduces to a HG mode: = − IG r HG r (third row in Table 1).To obtain an approximate expansion of mode

IG r
from the HG modes that is suitable for the whole range   , we use the Padé approximants [20].The Padé approximant of some function fz is given by the expression where the coefficients k a and k b are chosen so that at small values z , it has several first terms of the Taylor expansion by the powers of z , exactly the same as those of the function () fz, whereas at large z , several first terms of the asymptotic expansion by the powers (1 ) z are the same as those of the function ()  fz.Collecting equations (37) and (38) into the one, we get where exact representations of the functions of  are given by Equation (36), whereas their approximations in the form of Equation (39) depend on the choice of the parameters L and M. Since then the Padé approximants of the functions 3,0 () C  and 1,2 () C  should have the fol- lowing form: For each fixed value of the parameter L, the coefficients of the rational functions in Equation (41) can be found by solving systems of equations.These equations are ob-tained if we expand functions (41) into series in powers of  at 0  → and in powers of 1  at  → + , and then make equal these expansions to those already known from Equations (37) and (38).This approach is easy to implement in modern computer algebra systems.We here give only the simplest results when the obtained fractions are not very cumbersome yet: ) . 1 Differences between the exact functions 3,0 () C  and 1,2 () C  and their Padé ap- proximants are depicted in Figure 1.In particular, the maximal difference  −+ is 0.0187 for 2 L = and 0.0063 for 3 L = ; while the maximal difference  − is 0.0631 for 2 L = and 0.0224 for 3 L = .Here, the ex- pansions  − .Increasing the parameter L makes it possible to obtain in- creasingly accurate approximations of the IG mode, although this approach is accompanied by more cumbersome fractional rational expressions for the coefficients.
It is worth noting that using the Padé approximants does not require knowing the exact expressions for functions 3,0 () C  and 1,2 () C  .Instead, it suffices to know the first few terms of the series expansions of these functions (in powers of  at 0  → and in powers of 1  at  → + ).These expansions are obtained from the expansions of the roots of the characteristic equation; to expand the roots by the powers of  or 1  , well-known conventional methods are used.
In this section, the use of the Padé approximants is demonstrated for mode  (dashed) for 2 L = and 3 L = .We intentionally took different signs for both expressions so that the curves fall on opposite sides of the x-axis.

Numerical Simulation
In this section, we obtain intensity distributions for IG beams of different orders and different values of the parameter ε by solving the characteristic equation and by expansion into LG modes.For the computations, expressions ( 19), ( 25), ( 26), ( 30) and (34) are used, as are those in Table 1.ReLG with the amplitude being proportional to cos 3φ, i.e., the in- tensity contains six petals, all located on a single light ring.In addition, Figure 2 confirms that the intensity of the odd mode As is also seen in Figure 2, the intensity pattern of the even IG beam As seen in Figure 3, when the value of the parameter ε is small, IG modes align with their corresponding LG modes: the beam IG , but rotated by π/2.A comparison of Figure 2 for p = 3 and Figure 3 for p = 4 reveals that if p is even, then changing the sign of the parameter ε leads to rotation of the intensity pattern (Figure 3), whereas if p is odd, then changing the sign of ε leads not only to rotation of the intensity pattern, but also to a change in the pattern (Figure 2).This finding from Equation (21) 3).By contrast, for odd values of p, when ε changes sign, not only are the arguments swapped, but the functions are changed (an even to an odd, or vice versa):  As for smaller values of p, Figure 4 demonstrates that if the value of the parameter ε is small, then the IG modes align with the real and imaginary parts of the LG modes, and, on the contrary, if the value of the parameter ε is large by modulus, the IG modes reduce to one of the HG modes.Figure 5 also shows that when the parameter ε is large and positive, IG modes coincide with corresponding HG modes.Thus, the intensity distribution of the mode IG , but rotated by π/2.

Conclusions
In this work, we derived analytical expressions for 24 IG modes with the indices p = 3, 4, 5, and 6 via 2, 3, or 4 LG modes (more exactly, via their real or imaginary parts) and HG modes (see Equations ( 19), ( 25), ( 26), ( 30) and (34) and Table 1).The expansion coefficients of the IG modes in terms of LG and HG modes are expressed via the ellipticity parameter ε.Using the derived formulas, one can immediately derive expressions for the IG modes when the parameter ε is either zero or approaches plus or minus infinity.This explicit dependence of the IG modes on the ellipticity parameter enables control of their intensity through adjustment of this parameter.We have derived the symmetry properties of both even and odd IG modes.In particular, we have demonstrated how these modes are interrelated when the ellipticity parameter undergoes a sign change.The obtained representations of IG modes via HG and LG modes are convenient not only from the theoretical point of view, as they reveal the properties of these beams without simulation, but also from the practical point of view.This practical application makes it relatively easy to program complex amplitudes of the IG modes for their numerical simulation or generation using a spatial light modulator.For large-distance atmospheric probing using IG beams, high-power laser beams are needed.For this purpose, it is convenient to use a method of coherent beam combining, as proposed in [21].Ince-Gaussian beams can be effectively used for encoding (or protecting) data in wireless information transmission, where the same beam waist radius of the Gaussian beam.Then, in dimensionless variables, Equation (1) reads as follows: These formulae are very useful for investigating the limiting cases 0  → and  →.For 4 p = , the even, non-normalized IG modes are given by even IG modes (without normalization) are given by

Table 5 .
Limiting cases of normalized even IG modes at p = 5.Equation (22) establishes relationships between the roots of the characteristic polynomials: Therefore, as per Equation (21), odd IG modes are obtained from the even IG modes by swapping the indices of the HG modes, replacing  with  − , and adding a multiplier ( 1) m − : regardless of the in- dex p and the parity flag ( , ) eo , the coefficient for any term of the form   is an in- teger.Therefore, the series expansions of the roots(


IG r , which is expanded into the series of HG modes.It is easy to see that this ap- proach is suitable for any other IG mode.Its value for numerical computations of modes large values of the index p because, in this case, roots of the characteristic equation cannot be explicitly expressed in radicals.

Figure 2
Figure 2 illustrates intensity distributions of the IG modes at p = 3 for several values of the parameter ε.

Figure 2 .
Figure 2. Intensity distributions (here and in all figures below, black color means 0 and yellow-white colors mean maximum) of the IG modes

Figure 2
Figure 2 confirms that if the value of the parameter ε is small ( 0   ), the mode

3 , 3 o
IG also contains six petals, but they are ro- tated by an angle of π/6.

2 e 4 oIG has an intensity with 8 petals 4 e
with one central light spot and two surrounding light rings, and the beam 4,IG has an intensity with four petals on each of these two rings.The intensity of the beam petals, but they are divided into 4 petals on each of the two rings.The beam 4,located on a single ring.Additionally, Figure 3 confirms that for large values of the parameter ε, IG beams align with the HG beams.Therefore, intensity patterns of the IG beams at 1  ?align with respective HG beams.The intensity of the beam 4,0 e IG has 5 petals that are located on the vertical axis and separated by four horizontal zero-intensity lines; the intensity of the beam separated by two horizontal and two vertical zero-intensity lines; the intensity of the beam 4,IG has 5 petals, which are located on the horizontal axis and separated by four vertical zero-intensity lines; the beam 4,2 o IG has an intensity with 8 petals separated by one vertical and three horizontal zero-intensity lines, and the beam 4,4 o IG has the same intensity as the beam 4,2 o

3 . 5 Figure 4
Figure 4 depicts the intensity distributions of the IG beams with p = 5 for several values of the parameter ε.

Figure 5
Figure 5 illustrates the intensity distributions of the IG beams with p = 6 for several values of the parameter ε.

Figure 5 .
Figure 5. Intensity distributions of the IG modes

Figure 5 2 eee
Figure 5 confirms that if the ellipticity parameter ε is small, then the mode p, q values is transmitted, but the parameter ε is randomly changed.

Table 3 .
Limiting cases of the normalized even IG modes at p = 4.
Then, the non-normalized IG modes are expressed via the LG and HG modes in the following way: e m It is obvious that when index p, which determines the degree of the characteristic polynomial, increases, finding all of its roots in a simple form and without numerical methods becomes increasingly impractical.Nevertheless, it is easy to observe that for  → and  → + , we retain only terms up to ε 3 and ε −4 ,respectively, then we obtain the following asymptotics.
. This equation indicates that for even values p, the IG functions are not exchange of arguments, the intensity rotates by π/2 (Figure IGhas an intensity with three horizontal and three vertical zero-intensity lines, and, finally, the intensity of the mode eIGhas, vice versa, four vertical and two horizontal zero-intensity lines.The intensity of the mode