{'fraction_non_alphanumeric': 0.03621543870864632, 'fraction_numerical': 0.00923765382771296, 'mean_word_length': 4.879721856793836, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Complexity of Polarized Spatial Patterns in Large Area Square VCSEL

We consider pattern selection process in a wide aperture VCSEL near threshold. We show that for a square geometry of the laser aperture, the patterns formed at lasing threshold can be very complicated because of a possible misalignment between directions of an intrinsic spatial anisotropy of VCSEL and lateral boundaries of its aperture. The analogy with quantum billiard structures is established, and fingerprints of wave chaos are found. Influence of localized inhomogeneous in the pump current is also considered.


I. INTRODUCTION
In the last decades, semiconductor laser devices have played an increasing role in scientific research and applications. Among them, vertical cavity surface emitting lasers (VCSELs) can be distinguished, since these lasers emit normal to the wafer surface in the direction of epitaxial growth. One of the features of VCSELs design is a possibility to mount extremely large (up to hundreds µm in diameter) aperture with high level of spatial homogeneity. Spatial mode structures in such wide-aperture VCSELs have been a subject of many researches [1,2,3,4,5,6,7,8,9,10,11].
In general, the mechanism of pattern formation in VC-SEL is very complicated and involves both a complex structure of a VCSEL cavity (including Bragg reflectors) [10,11,12], as well as peculiarities of light-matter interaction in active quantum well semiconductor layers [15,16,17,18,19,20]. However, near threshold one can invoke the perturbation theory and obtain normal forms governing the evolution of the system. It was shown recently, that due to a slight spatial anisotropy of a VCSEL cavity only a few spatial modes come into play at threshold and the shape of these modes can be analyzed by the linear approximation of the normal form [10,11,12] (in contrast to spatially isotropic systems, where the whole degenerate family of modes have the same critical growth rate at threshold, and the selection process requires consideration of nonlinear competition even at threshold [21,22]). The investigations of the problem using this point of view was started in [12] for VCSELs with circular aperture. In that work, transition from the 'flower-like' modes dictated by circular boundaries to 'stripe-like' ones, which are required by spatial anisotropy of a VCSEL, was described.
In [12], the description was restricted to the light linear polarized in direction coinciding with one of axes of the intrinsic laser anisotropy. In the present article, we extend this approach to the more general "vectorial" case when the laser anisotropy is not so strong and both orthogonally polarized components of the field must be taken into account. This extension allows us to investigate the interaction of two polarization degrees of freedom, that is important for the consideration of competition of different mechanisms affecting pattern formation.
The current research is motivated by the recent experimental investigations [12,13,14], where it is shown that many aspects of spatial structures in VCSEL such as their transverse shape or frequency-versus-lengthscale dependence can be well described in a linear approximation, which validates the 'linearized normal form' approach, used in this article.
We show here that this makes the structures at the laser threshold very complicated even in the case of a simple square aperture and perfectly manufactured VC-SEL cavity. The similar complexity was observed in the resent experimental investigations [4,23,24], where the scarred structures in the square broad-area VCSEL were obtained and treated as coherent states of quantum billiard. So, it is possible to speak in this context about 'quantum chaos' in a wide aperture VCSEL.
The 'quantum chaos' (or 'wave chaos') approach, including the investigation of level (or eigenfrequency) distribution in a domain with 'integrable' or 'nonintegrable' boundary conditions, was applied for investigations of microwave billiards [25,26,27], as well as semiconductor microdisk lasers [28,29]. Very recently, the wave chaos was observed in a VCSEL with deformed circular cav-ity geometry [30]. Though the laser emission far from threshold is strongly affected by nonlinearities, the optical spectra obey the same relations as a quantum billiard, governed by a linear Schrödinger equation.
We develop this approach further by considering the wide-aperture VCSEL near threshold as a kind of quantum billiard where the role of Hamilton operator is played by the operator of linearized order parameter equation for VCSEL. We show, that this operator is considerably more complicated than the transverse Laplace operator which is often used for consideration of quantum billiard systems. This leads to appearance of quantum chaos-like features, such as scarred orbits and nongaussian statistics of the levels even in the simple square geometry without a deformation of boundary conditions, that usually used for demonstrating quantum chaos features.
The structure of the article is following: in the section II we extend the approach used in [12] to take into account both orthogonally polarized components of the light field. In the section III we consider the pattern formation in VCSEL with a square geometry and discuss the influence of the polarization axes rotation against lateral boundaries of the cavity, as well as of slight spatial perturbations of the pumping current. In the section IV we interpret our results in the 'quantum chaos' framework.

II. THE VECTORIAL EIGENVALUE PROBLEM
As has been pointed in the Introduction, the eigenmodes and their decay (or growth) rates can be directly obtained from a linear stability analysis of the nonlasing zero solution of the nonlinear equations governing VCSEL's dynamics [10,11,12]. In comparison with a method when the field is decomposed into transverse eigenfunctions of the empty cavity [9,31,32], this approach takes into account properties of a gain medium and of Bragg mirrors composing the cavity by the strick way.
For the infinitely large transverse-area devices the resulting eigenmodes are plane tilted waves (transverse Fourier modes) with the decay rates depending on the detuning of the longitudinal cavity resonance from the gain maximum and on a polarization of modes [10,11,33]. In contrast, the corresponding modes for the finite devices may be very complicated even in the scalar case. For example, for the circular aperture the modes can range from Bessel-like modes to the modes resembling more tilted waves [12].
The basic nonlinear equations for the vectorial case were obtained in [11] . The short description of these equations is given in Appendix A. The linearization procedure is a generalization of one introduced in [12] to vectorial case, and is described in Appendix B. As in the scalar case, we linearize a complex nonlocal nonlinear operator, obtaining the linear but nonlocal pseudodifferential operator (or speaking more precisely, certain eigen-value problem for such operator) governing the field at threshold. Then, we approximate the operator describing the action of Bragg reflectors in transverse Fourier space and obtain the operator containing only partial derivatives up to fourth order. The total operator obtained after such procedure includes main peculiarities of pattern formation such as the dependence on detuning between longitudinal resonance of the cavity and gain peak frequency, and on the anisotropy of Bragg reflectors. This operator has the following form in the (x, y) space: In contrast to the scalar case, for the vectorial case operatorÔ is a matrix, acting on the vectorial field E = (E x , E y ). Therefore, the coefficients a ij are also 2 × 2 matrices. We also consider inhomogeneities in the pumping current δµ and in the index δn, which play the role of disturbances. In the following, we mainly restrict ourself to only current inhomogeneities, because they are easier to be controlled.
Eigenfunctions E g and eigenvalues λ g of this operator are obtained by solving an eigenvalue problem: The eigenfunction corresponding to the eigenvalue with the largest real part determines a spatial field distribution with maximal growth rate at threshold appearing after onset of generation. Aperture of the device has been simulated by zero field conditions on the boundaries. Besides, because the operatorÔ is of fourth order, one should introduce an additional boundary condition which contains spatial derivatives of the field (see Appendix B) to obtain the well posed eigenvalue problem. This second boundary condition can be chosen by different ways, but all them leads to approximately the same result [12]. Therefore, we select the second boundary condition appropriately, to simplify the solution (see Appendix C).
The eigenvalue problem (2) is a generalization of more conventional one with transverse Laplace operator (6). The important difference however, is that the operator (1) is anisotropic, with a principal directions defined by polarization anisotropy of VCSEL cavity. It should be noted that x-component of the field E x and y-component of the field E y enters to Eq. (2) as just two components of a single vectorial eigenfunction E g . In the following, we suppose that the x component of the field is always strongest one due to the intrinsic polarization anisotropy. However, the orthogonal component can not vanish because of the mixing of both components during propagation in the cavity.

A. Influence of alignment of boundaries and ansotropy direction
To elucidate the mechanism governing the pattern selection, we first consider the case of homogeneous pumping and index profile (δµ = 0, δn = 0). The presence of the aperture is modeled by zero boundary conditions. In this subsection we show, that the resulting spatiotemporal distribution strongly depends on the alignment of the boundary conditions and the main anisitropy axes of the device.
When boundaries and anisotropy directions are completely aligned, the x-component of the first mode at threshold is a standing wave with stripes parallel to the direction of its polarization (see Fig. 1(a), (b)). Though the y-component is not zero due to the polarization mixing effect by Bragg reflectors, the coupling between x and y polarizations is very week since for this case the Fourier transform of homogeneous part ofÔ, which is (k x , k y )-dependent 2x2 matrix, have nearly zero diagonal elements on the line k y = 0 and k x = 0. Thus, y-component is weaker in 100 times approximately for taken parameters.
If the boundaries are rotated with respect to the polarization anisotropy direction, the stripes remain for a small rotation angle α, and their direction coincides with the boundaries rather then with anisotropy direction, i.e. they are rotated with the boundaries (see Fig. 1(e)-(h)). Structures of both polarized components are similar and their intensity become comparable.
However, as α reaches some critical value, the corresponding structure becomes more complicated. Thus, for the parameter of Fig. 1(i)-(l), the spatial distribution of the x-polarized component of the field resembles stripes in the middle of the aperture, whereas near boundaries the pattern is more 'square-like' one. At the same time, the structure of y-polarized component is more complicated and less regular. The critical angle for which stripes are giving place to complicated structures decreases with the size of the device, as well as with a value of anisotropy γ a . However, stripes take place for relatively large size (composing of several tens of oscillation of the intensity).
It should be noted that the structures presented in Fig. 1(i)-(l) are the combination of many transverse Fourier harmonics of different directions (see Fig. 1(j),(l)) in contrast to two travelling waves composing stripes for small α ( Fig. 1(b),(f)).We show in the next section by considering the corresponding eigenvalue statistics, that such structures demonstrate some features which make them close to quantum chaos wavefunctions .

B. Influence of inhomogeneities
In this subsection, we consider an influence of pump inhomogeneities of different shape, amplitude and localization. As it was found earlier for devices with circular aperture [34], any slight inhomogeneities can change the sequence of a few first eigenfunctions in accordance with the order of their decay rates, preserving the shape of modes and their frequencies. In this connection, it is useful to investigate the subsequent eigenmodes of the linear eigenvalue problem (2). The second eigenmode for the parameters of Fig. 1(e)-(h) is shown in Fig.2. The Fourier spectrum of this mode for both polarization components consists of four points: two pairs of two closed points. These pairs are situated approximately at the same places that two spots in Fig.1 (f), (h) and rotated by some angle against polarization direction. So, the main difference between two first modes is the splitting of the Fourier harmonics leading to the strong modulation of stripes in the near field ( Fig.2 (b), (d)). The second eigenmode can be made leading one (one with the highest growth rate) by adding small inhomogeneity to the current δµ = 0.05δ(x, y), where δ(x, y) is a probe function which is nonzero in the vicinity of the right border and zero everywhere else.
If we decrease the detuning δ between the gain maxima and the cavity resonance, the gain-loss dispersion mechanism of Fourier harmonics selection becomes weaker [10,11] and, as result, eigenmodes are more strongly determined by boundary conditions. The first mode for the large angle near π/4 does not have clear stripelike structures in the center of device ( Fig. 3 (a),(b)), and the second mode has the structure even less similar to stripes than previous one (Fig. 3 (c),(d)). However, stripelike structures do not disappear at all, as one can see in Fig. 3 (e),(f) for the next eigenmode. Moreover, the patterns with orthogonal direction of stripes can be easily exited ( Fig. 3 (g),(h)). This could be explained by the fact that these two families of stripes are degenerate when α = 0, and start to deviate from each other with increasing α. Some subsequent eigenmodes with larger decay rate for the same parameters are presented in Fig. 4. It is evident that for small enough (but nonzero) anisotropy the modes can differ from each other in their symmetry properties essentially, as in the case of circular aperture [12]. One can find nearly unordered structures (Fig. 4 (e)-(h)) and the structures which can be considered as 'scared' eigenmodes of quantum billiard (Fig. 4 (a)-(d), (i)-(l)).
It is difficult to change the order of subsequent modes (especially to push one of them to the first position with highest grough rate) by introducing a small inhomogeneity. Since stronger inhomogeneities needed for that mode structures can be changed also.
It should be also noted that whereas patterns of the orthogonally-polarized component (weaker component) of the laser field are usually of the same shape, one can find some exceptions (compare Fig. 4(i) and Fig. 4(k)). However, even in this case the transverse Fourier images of the field components look quite similar. This is connected to the fact, that structures in both polarizations belong to the same eigenmode and coupled to each other in the absence of strong inhomogeneities through boundary conditions (they both and some combination of their derivatives must vanish at the boundary), and through the Bragg reflectors. The exception is only when this connection is very weak (Fig. 1 (a)-(d)). Up to now we have investigated the eigenmodes of homogeneous devise, having in mind that the grough rate of a mode following the first one can become the largest one when a small perturbation of the current δµ is added, which does not sufficiently change the shape of mode.
Increasing the intensity of inhomogeneities leads to changes in the shape of eigenmodes. As an example, we consider here strong inhomogeneities by introducing a peak or a hole into the laser aperture, that is quite natural for experiments.
For the peak in current profile, the first eigenmode is defined by this spot only, as if the rest of aperture is empty (Fig. 5 (a)-(d)). But, the subsequent modes fill all the cavity as there is no peak in the aperture at all (Fig. 5 (e)-(h)). However, they are more disordered then in the absence of inhomogeneity. The corresponding eigenvalues are also more separated as it is described in the next section. On the other hand, if there is a strong hole in the current profile (which is for simplicity modeled by introducing zero boundary conditions on the circle surrounding the hole), the effect is not so visible, i.e. there is not a separation on the "hole" mode and the rest ones. Moveover, the shape of eigenmodes for x polarization is not changed so dramatically far from the hole (Fig. 6), although the structure of weak component is strongly disordered comparing to the nonperturbed case ( Fig. 1(c)). However, the subsequent eigenmodes are more irregular that the first one even for the x-component. The problem becomes close to chaotic one that could be considered in the framework of wave chaos.

IV. SIGNATURES OF QUANTUM CHAOS
Generally speaking, the term 'quantum chaos' or 'wave chaos' is usually attributed to investigation of quantum systems, which posses chaotic features in classical limit [25]. In contrast to classical systems, they obey linear evolution equation (Schrödinger equation). So, the object of investigation is usually the set of eigenvalues and eigenfunctions (stationary states) of some linear operator, describing the quantum system under consideration (usually Hamilton operator). At that, the very important role is played by eigenvalue distribution, which shows fingerprints of chaos in quantum systems when the corresponding classical systems are chaotic.
Among the quantum systems, much attention is devoted to the billiard systems [25,35]. The quantum coun-terpart of classical billiard is described by a wavefunction ψ(x, y, t), which obeys Schrödinger equation: with boundary condition on the boundary of the domain S. Stationary states of this problem are eigenfunctions of the transverse Laplace operator ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 : Depending on the shape of domain S, the corresponding classical problem can be either chaotic or regular. In the limit case of integrable system (for example, when the region has a rectangular shape), the corresponding eigenvalue problem gives set of eigenvalues, obeying Poisson statistics of an eigenvalue spacing s: where p(s) is a probability for corresponding value of eigenfrequency distribution s i = Im λ i − Im λ i+1 . For the opposite limit case of completely chaotic system, the corresponding distribution is the Wigner one: For intermediate situations one can define other statistic families limited by Eq. (7) and Eq. (8) [25,36,37]. One of fingerprints of quantum chaos in these intermediate cases, as well as in the Eq. (8), that the maximum of p(s) is reached for s = 0. This defines so called level repulsion phenomena for chaotic systems.
The eigenfunctions of Hamilton operator also demonstrate fingerprints of chaos. In particular, one can observe so called 'scarred' patterns [23,24,38,39,40], which are localized near (unstable) periodic trajectories of corresponding classical systems. Among the scared patterns, there are other eigenfunctions, which are strongly irregular and fit all the area.
It should be noted, that besides the quantum systems, the above mentioned framework is of the considerable interest for macro-systems, which are described at some level of approximation by Eq. (6) or Eq. (4). Among them are microwave billiards [25,26,27], microdisk lasers [28,29] and VCSELs [23,24,30].
The eigenvalue statistics and other fingerprints of quantum chaos can give a criterion of a complexity of patterns near threshold. Indeed, if a system possess the Poisson statistics, the corresponding eigenfunctions are regular and stable against small perturbations, whereas for the Wigner one the eigenfunctions are much more irregular.
In the present work, we describe the pattern formation in VCSEL by linear operatorÔ (1) obeying Eq. (2) This operator is considerably more complicated then the transverse Laplace operator from Eq. (6). First of all, the eigenvalues of the operatorÔ are complex. Their real parts describe the decay rate of the corresponding eigenfunctions, whereas the complex parts are their oscillation frequencies determining an energy level distribution. Examples of the eigenvalues distribution are presented in Fig. 7. For homogeneous device, all eigenvalues are concentrated along one line and in average their decay rates are increased with the frequencies (Fig. 7 (a)). The closeness of decay rates of the adjacent modes allows to change their order by a small perturbation. It is worth noting that for the taken parameters the eigenvalues are grouped in clusters and become deviate with decreasing of the angle α. In the case of a local peak in the pump current the eigenvalues distribution confirms the separation of the first and the rest modes ( Fig. 7 (b)). On the other hand, the operator (1) is anisotropic, with a preferential directions defined by anisotropy of VCSEL cavity. As one can see from (1), (2), if the boundaries are defined parallel to x and y axes, one can separate variables in the eigenvalue problem and therefore it becomes integrable. Hence, as one can see in Fig. 8(a), when the boundaries and anisotropy direction coincides, the eigenvalues obey more or less Poisson statistics.
When these two directions (anisotropy and boundaries for square shape of the aperture) are not aligned, the variable separation is not possible anymore (as would be for the second order differential operator). As a result, for α = 0 ( Fig. 8(b)) one can clearly see that the maximum probability of the level spacing distribution is shifted to nonzero s, possessing a fingerprint of wave chaos. However, the system can not be considered as purely 'quantum chaotic' because the statistic is of intermediate type, neither purely Poisson nor Wigner one. It should be noted, that for isotropic system obeying Eq. (6) the eigenvalue spacing statistics remain Poisson for any values of α. In our case, the eigenvalue spacing has the maximum value again at zero for the angle α = π/4, that can be explained by possible degeneracy of the operator O (Fig. 8(c)).
As a comparison the eigenlevel statistics is shown for boundary conditions of Fig. 6. It is clearly seen that for this case the statistics for spacing of imaginary parts of eigenvalues can be considered as Wigner one, which is quite expectable, because the same happens with the eigenvalues of the transverse Laplace operator. On the other hand, for the large peak in the aperture the eigenvalues do not obey the Wigner statistics anymore. For this case, the leading eigenvalues (eigenfunctions for which are shown in Fig. 5(a)-(d) and can be described as eigenvalues of a disturbance only) are strongly separated from others (which are resembling Fig. 5(e)-(h)), which obviously breaks the Wigner distribution.
The other signatures of complexity of spatial structures can be found in fingerprints of quantum chaos demonstrated by the shape of eigenfunctions. As one can see, the neighboring eigenmodes for misaligned boundaries and anisotropy directions (α = 0) can differ sufficiently (Fig. 3, Fig. 4), including very irregular ones (Fig. 4(e)-(g)), which is a common property of systems possessing quantum chaos. On the other hand, beside complicated unordered structures, quantum chaos is characterized by 'scared' eigenfunctions, which are located near unstable periodic trajectories of corresponding classical system. The example of eigenfunctions of operatorÔ, resembling that kind of structures (Fig. 4)(a)-(d),(i)-(l). Experimentally these structures were observed recently in VCSEL [23,24]. Bragg reflectors, the presence of polarization anisotropy creates a spatial anisotropy. This allows to investigate the structure of pattern near threshold by linearized order parameter equations, since even in a linear approximation the mode with the smallest decay rate is nondegenerate (in contrast to spatially isotropic systems, where the whole family of modes has the same growth rate at threshold, and nonlinear competition is of strong importance for the selection). Since our studies concentrated on a VCSEL with square aperture, the competition of boundary direction (which can be defined as a vector, parallel to one of the sides of the square), and spatial anisotropy is important mechanism for pattern selection.
To measure the complexness of a pattern structure we use criteria from the topic of quantum chaos, where the main object of research is also an eigenvalue problem of some operator. One of fingerprints of complicated behavior is an eigenvalue spacing statistics.
In this work we have considered a statistics of imaginary parts of eigenvalues, describing frequencies of corresponding eigenfunctions. We have shown that in the case of aligned directions of boundaries and spatial (or polarization) anisotropy, the above mentioned statistics is close to Poisson one. However, when the angle of rotation of spatial anisotropy against boundaries is not zero, the statistics is not Poisson anymore. In the later case, the statistics shows a signature of level spacing repulsion, i.e. maximum of eigenspacing distribution is not zero.
The other fingenprints of quantum chaos can be found in the spatial shape of eigenfunctions, which can be very irregular, as well as 'scared-like' patterns.
To compare these results with experiments, we also have investigated the influence of the inhomogeneities. The considered system is stable in the sense, that small perturbations only change the modes order preserving their shape, and not lead to dramatic change of the whole picture. In this sense, the consideration of subsequent modes (having higher decay rate) are also important, because they can be considered as leading ones of the device perturbed by some (may be unknown) inhomogeneity.
With increasing the inhomogeneity, the shape of eigenmodes also changes, as well as all the characteristics like the level spacing distribution. One can, for example, model the hole in the current by imposing zero boundary conditions somewhere on a circle located inside the squire aperture. In this case, the statistics evidently becomes Wigner, because the corresponding billiard is scattering one. On the other hand, large positive peak in the current leads to separation of eigenvalues and eigenfunctions into two classes, first of them belonging to 'peak family', with eigenfunctions localized near the peak, whereas the second family is distributed across the whole cavity aperture.

Acknowledgments
This work was financially supported by the Deutsche Forschungsgemeinschaft for equipment and by travel grants. We are grateful to Malte Schutz-Ruhtenberg for useful discussions.

APPENDIX A: BASIC SYSTEM
The system under consideration was initially obtained in [11]. The equations for the field E = (E x , E y ) take into account propagation inside the complex VCSEL cavity. The material equations are based on spin-flip model [41], describing four level system with populations differences of relevant transitions determined by their sum N and difference n: and the operatorsM andĜ describes the modal losses and gain, correspondingly,Ω describes the influence of the diffraction in the complex resonator of VCSEL. These operators take into account the influence of reflection from Bragg reflectors and are quite complicated. More details can be found in [11,12]. The operator describes the Lorentz shape of gain contour with a detuning δ between the gain maximum and cavity resonance frequency. κ is the field decay rate in the VCSEL's cavity. The operatorΓ describes the inner anisotropy of VC-SEL's cavity: where γ a and γ p is an amplitude and phase anisotropy, respectively.

APPENDIX B: VECTORIAL EIGENVALUE PROBLEM DERIVATION
To derive the linear evolution equation in the form, (which is then reduced to Eq. 2) using the basic equations (A1-A3), we take into account that at the laser threshold two branches of the steady state solutions cross each other: the zero solution (E = 0), and the nonzero lasing one. Because of the lasing solution at the cross-section point is characterized by E = 0 and N = µ the further analysis is drastically simplified giving the following diagonal operator of the linearized problem: where eachÔ ii is in turn the matrix operator: The diagonal form of (B2) shows, that the field E and the variables N ,n are independent of each other at the laser threshold, and the spectrum forÔ 22 lies entirely in the half-plane Re λ ≤ 0, hence we can consider only the equation for the field Eq. B1 withÔ ≡Ô 11 .Ô includes the current density at threshold N = µ. From (A1),(B2), the explicit value for the operatorÔ is defined by the following action on the field E: However, the expression (B4) is still nonlocal, because all the operators in (B4) are integrodifferential. The equation can be made local in the manner of [12]. The procedure used for a scalar operator in [12], must be applied to each of the operators O (11) ij , which are the scalar components ofÔ. Because the diagonal part of (B3) plays the main role, it is approximated up to the fourth order of k ⊥ , whereas the operators O (11) ij for i = j are approximated with terms of the second order of k ⊥ .
In the presence of index inhomogeneities,Ω = Ω homogen + ilδn (where l = τ ω/n 0 is expressed via the optical frequency ω, round trip time of the cavity τ and the mean index n 0 . The pump current inhomogeneities have to be included in the coefficient of the last term µ = µ homogen + δµ.

APPENDIX C: NUMERICAL IMPLEMENTATION OF EIGENVALUE PROBLEM
To solve the problem, the Matlab PDE toolbox implementation [42] of Arnoldi method for the system of partial differential equations of the second order has been used. For that, the operators of the fourth orderÔ (11) ii have been represented via multiplication of two operatorŝ P i1 ,P i2 of the second order with constant coefficients: where P i1 = ∇ · (c Here c (k) ij are matrices defined by formulas given in [12] using the coefficients a ij (to define c (1) ij and c (2) ij the coefficients (a ij ) (11) and (a ij ) (22) are used, correspondingly). Then, defining a new vector variable one can present the action of the operator of fourth order O (x,y) as a operator of the second order of the typê P = ∇ · (c ⊗ ∇) + a.
acting on a four-component function E ex = (E, E 1 ).
Here a is a 4 × 4 matrix (instead of 2 × 2 for scalar case [12]), c is a rank four tensor which can be described by eight 2 × 2 matrices c (k) ij , mentioned above after equation (C1). The boundary conditions have to be also generalized for the function E ex as in the work [12]: where E 1 is defined by (C2) or (C3).
In the framework of the original system it means that we introduce a second boundary conditions on the vector field E in accordance with that the operatorÔ (x,y) is of the fourth order.