Full vectorial analysis of polarization effects in optical nanowires

We develop a full theoretical analysis of the nonlinear interactions of the two polarizations of a waveguide by means of a vectorial model of pulse propagation which applies to high index subwavelength waveguides. In such waveguides there is an anisotropy in the nonlinear behavior of the two polarizations that originates entirely from the waveguide structure, and leads to switching properties. We determine the stability properties of the steady state solutions by means of a Lagrangian formulation. We find all static solutions of the nonlinear system, including those that are periodic with respect to the optical fiber length as well as nonperiodic soliton solutions, and analyze these solutions by means of a Hamiltonian formulation. We discuss in particular the switching solutions which lie near the unstable steady states, since they lead to self-polarization flipping which can in principle be employed to construct fast optical switches and optical logic gates.


Introduction
The Kerr nonlinear interaction of the two polarizations of the propagating modes of a waveguide leads to a host of physical effects that are significant from both fundamental and application points of view. Here, we develop a model of nonlinear interactions of the two polarizations using full vectorial nonlinear pulse propagation equations, with which we analyze the nonlinear interactions in the emerging class of subwavelength and high index optical waveguides. Based on this model we predict an anisotropy that originates solely from the waveguide structure, and which leads to switching states that can in principal be used to construct optical devices such as switches or logical gates. We derive the underlying nonlinear Schrödinger equations of the vectorial model with explicit integral expressions for the nonlinear coefficients. We analyze solutions of these nonlinear pulse propagation equations and the associated switching states by means of a Lagrangian formulation, which enables us to determine stability properties of the steady states; this formulation provides a global view of all solutions and their properties by means of the potential function and leads, for example, to the emergence of kink solitons as solutions to the model equations. We also use a Hamiltonian formalism in order to identify periodic and solitonic trajectories, including solutions that allow polarization flipping, and find conditions under which the unstable states and associated switching solutions are experimentally accessible. The nonlinear interactions of the two polarizations of the propagating modes of a waveguide have been studied extensively over the last 30 years [1]- [13]. Different aspects of the inter-actions have been investigated, for example Stolen et al. [1] used the induced nonlinear phase difference between the two polarizations to discriminate between high and low power pulses. In the context of counterpropagating waves, the nonlinear interactions have been shown to lead to polarization domain wall solitons, [8]- [10], [14] which are described as kink solitons representing a polarization switching between different domains with orthogonal polarization states. The nonlinear interactions can also lead to polarization attraction [9], [11]- [13], [15,16] where the state of the polarization of a signal is attracted towards that of a pump beam. For twisted birefringent optical fibers, polarization instability [2,5] and polarization domain wall solitons [17] have been reported. The nonlinear interactions also induce modulation instability which results in dark-soliton-like pulse-train generation [6,7]. Large-signal enhanced frequency conversion [18], cross-polarization modulation for WDM signals [10], and polarization instability [3] have also been reported and attributed to nonlinear polarization interactions. Stability behavior has been studied in anisotropic crystals [19].
The nonlinear interactions of the two polarizations can also be studied in the context of either nonlinear coherent coupling or nonlinear directional coupling in which the amplitudes of two or more electric fields, either the two polarizations of a propagating mode of a waveguide or different modes of different waveguides, couple to each other through linear and nonlinear effects [20]- [22]. Nonlinear directional coupling is relevant to ultrafast all-optical switching, such as soliton switching [23]- [27] and all-optical logic gates [28]- [30]. The interaction of ultrafast beams, with different frequencies and polarizations, in anisotropic media has also been studied and the conditions for polarization stability have been identified [24,31].
In previous work ( [32], Chapter 6), the nonlinear interactions of the two polarizations are described by two coupled Schrödinger equations. These equations employ the weak guidance approximation, which assumes that the propagating modes of the two polarizations of the waveguide are purely transverse and orthogonal to each other within the transverse x, y plane, perpendicular to the direction of propagation z. Based on this, the electric fields are written as where A i (z,t) are the amplitudes of the two polarizations, with e 1 (x, y) e 2 (x, y) = e 1 (x, y)e 2 (x, y)x ŷ = 0, where e 1 (x, y), e 2 (x, y) are the transverse distributions of the two polarizations,x,ŷ are unit vectors along the x and y directions, and it is understood that fast oscillatory terms of the form exp(−iωt ± β i z) are to be included for the polarization fields. The weak guidance approximation also assumes that the Kerr nonlinear coefficients for the self phase modulation of the two polarizations are equal because their corresponding mode effective areas are equal [32]. We refer here to models of nonlinear pulse propagation based on the weak guidance approximation simply as "scalar" models, since these models consider only purely transverse modes for the two polarizations. The weak guidance approximation works well only for waveguides with low index contrast materials, and large dimension structure compared to the operating wavelength. This approximation is, however, no longer appropriate for high index contrast subwavelength scale waveguides (HIS-WGs) [33]- [35]. These waveguides have recently attracted significant interest mainly due to their extreme nonlinearity and possible applications for all optical photonicchip devices. Examples include silicon, chalcogenide, or soft glass optical waveguides, which have formed the base for three active field of studies: silicon photonics [36]- [40] In order to address the limitations of the scalar models in describing nonlinear processes in HIS-WGs, we have developed in [33] a full vectorial nonlinear pulse propagation model. Important features of this model are: (1) the propagating modes of the waveguide are not, in general, transverse and have large z components and, (2) the orthogonality condition of different polarizations over the cross section of the waveguide is given by e 1 (x, y) × h * 2 (x, y) ẑ dA = 0, rather than simply e 1 (x, y) e 2 (x, y) = 0 as in the scalar models. These aspects lead to an improved understanding of many nonlinear effects in HIS-WGs; it was predicted in [33], for example, that within the vectorial model the Kerr effective nonlinear coefficients of HIS-WGs have higher values than those predicted by the scalar models due to the contribution of the zcomponent of the electric field, as later confirmed experimentally [46]. Similarly, it was also predicted that modal Raman gain of HIS-WGs should be higher than expected from the scalar model [49].
Here, we extend the vectorial model to investigate the nonlinear interaction of the two polarizations of a guided mode. The full vectorial model leads to an induced anisotropy on the dynamics of the nonlinear interaction of the two polarizations [50], which we refer to as structurally induced anisotropy, in order to differentiate this anisotropy from others, such as those for which the anisotropy originates from isotropic materials. The origin of the anisotropy is the structure of the waveguide rather than the waveguide material.
The origin of this anisotropy in subwavelength and high index contrast waveguides has also been reported by Daniel and Agrawal [35], who considered nonlinear interactions of the two polarizations in a silicon rectangular nanowire including the effect of free carriers. In their analysis, however, they ignore the coherent coupling of the two polarizations, considering the dynamics of the Stokes parameters only for a specific waveguide and ignore the linear phase.
This anisotropy in turn leads to a new parameter space in which the interaction of the two polarizations shows switching behavior, which is a feature of the vectorial model not accessible through the scalar model with the underlying weak guidance approximation. We also show that the resulting system of nonlinear equations, for the static case, can be solved analytically. Due to the underlying similarity of the nonlinear interaction of the two polarizations and the nonlinear directional coupling of two waveguides, the anisotropy discussed here can be also applied to the case of nonlinear directional coupling, in which the two waveguides have different effective nonlinear coefficients for the propagating modes.
This work develops and expands on results reported for the first time in [50, 51], in particular we derive here (in Section 2) the equations that describe the nonlinear interactions of the two polarizations within the framework of the vectorial model, with explicit integral expressions for the nonlinear coefficients. In Section 3 we determine properties of the static solutions, classify the steady state solutions, and determine their stability using a Lagrangian formalism. We also discuss a Hamiltonian approach and how the phase space portrait provides a complete picture of the trajectories of the system, including the periodic and solitonic solutions (Section 3.5). We derive analytical periodic solutions by direct integration of the system of equations in Section 4, and then discuss switching solutions and their properties. We relegate to the Appendix a mathematical analysis of the exact soliton solutions, which are relevant to the switching solutions, with concluding remarks in Section 5.

Nonlinear differential equations of the model
In the vectorial model the nonlinear pulse propagation of different modes of a waveguide is described by the equations: where µ, ν = 1, 2 with µ = ν, and A 1 (z,t), A 2 (z,t) are the amplitudes of the two orthogonal polarizations. These equations follow from the analysis in [33], by combining Eqs. (23,32) of [33], but without the shock term. The linear birefringence is defined by ∆β ν µ = −∆β µν = β ν − β µ and the γ coefficients are given by Here we use the notation (e ν ) 2 = e ν e ν , |e ν | 2 = e ν e * ν and e 2 In these equations e 1 (x, y), e 2 (x, y) are the modal fields of the two orthogonal polarizations, k = 2π/λ is the propagation constant in vacuum, and γ ν , γ µν , µν are the effective nonlinear coefficients representing, respectively, self phase modulation, cross phase modulation, and coherent coupling of the two polarizations, and is the normalization parameter. The coupled equations (2) describe the full vectorial nonlinear interaction of the two polarizations. There are two fundamental differences between these equations and the typical scalar coupled Schrödinger equations (see for example Chapter 6 in [32]). Firstly, the additional terms A * µ A 2 ν , A µ |A ν | 2 , A µ A µ 2 on the right hand side of Eq. (2) represent interactions between the two polarizations. These do not appear in the scalar model since the effective nonlinear coefficients associated with these terms, γ µν as given in Eqs. (6,7,8), contain factors such as e µ e ν which are zero in the scalar model, since the modes are assumed to be purely transverse. All possible third power combinations of the two polarization fields, namely due to the z-component of the modal fields. Secondly, in all effective nonlinear coefficients given by Eqs. (3)(4)(5)(6)(7)(8), the modal fields e and h have both transverse and longitudinal components, unlike the scalar model in which modal fields have only transverse components. The terms containing nonzero e µ e ν provide a mechanism for the interaction of the two polarizations since they allow for exchange of power between the two modes through the z-components of their fields. The last term on the right hand side of Eq. (2), for example, indicates a coupling of power into a polarization, even if initially no power is coupled into that polarization. Although the terms on the right hand side of Eq. (2) that contain e µ e ν are nonzero, they are generally significantly smaller than the remaining terms and are therefore neglected in the following; further investigation of the effects of these terms, and a discussion of their physical significance, will be presented elsewhere. The focus of this paper is to investigate the effect of the z-components of the fields e and h, which influence the values of the effective coefficients, and therefore also the nonlinear interactions of the two polarizations. Hence, from (2), we obtain the equations: These are similar in form to the scalar coupled equations ( [32], Section 6.1.2), however, the coefficients γ ν , γ µν , γ ′ µν , given in Eqs. (3)(4)(5), now contain z-components of the electric field, through both e and h. In the framework of the scalar model, the weak guidance approximation assumes that the effective mode areas of the two polarization modes are equal [32], leading to where we have denoted γ c = γ 12 = γ 21 , γ ′ c = γ ′ 12 = γ ′ 21 . This means that in the scalar model there is an isotropy of the nonlinear interaction of the two polarizations; in order to break this isotropy, one needs to use either anisotropic waveguide materials or twisted fibers, or else couple varying light powers into the two polarizations by using either counter-or co-propagating laser beams. The fact that in the vectorial form (10) of the coupled equations the γ values include the z-component of the fields, as given by Eqs. (3)(4)(5), means that the equalities (11) do not hold in general. As an example, see Fig. 1 in [50] which plots γ 1 , γ 2 , γ c , γ ′ c for a step-index glass-air waveguide with an elliptical cross section; evidently the equalities (11) are not satisfied. One consequence of the vectorial formulation is, as we show in Section 3.4, the existence of unstable states not present in the scalar formulation.

Static equations
We find now all solutions of Eq. (10) for the static case, in which the fields A 1 , A 2 are functions of z only. We have therefore the two equations where ∆β = β 1 − β 2 . We express the fields A 1 , A 2 in polar form according to where the powers P 1 , P 2 and the phases φ 1 , φ 2 are real functions of z. It is convenient to define the phase difference ∆φ and an angle θ according to then upon substitution into Eqs. (12,13) we obtain the four real equations: The last equation decouples from the remaining equations, hence we first solve Eqs. (16)(17)(18) for P 1 , P 2 , θ and then determine φ 1 by integrating (19). Eqs. (16,17) show that P 0 = P 1 + P 2 is constant in z. We define the dimensionless variables and the dimensionless parameters In terms of these parameters we obtain the two equations: Since τ takes only positive values, we may regard τ as a time variable which is limited in value only by the length of the optical fiber and by the value of P 0 , and we set the initial values v 0 = v(0), θ 0 = θ (0) at time τ = 0, i.e. at one end of the fiber. The general solution depends on the initial values v 0 , θ 0 and on only two parameters a, b, even though Eqs. (16)(17)(18)(19) depend on the five constants P 0 , γ 1 , γ 2 , γ c , γ ′ c . At the initial time we have P 1 , P 2 > 0 and so we always choose v 0 such that 0 < v 0 < 1. It may be shown from Eqs. (22)(23) that 0 < v(τ) < 1 is then maintained for all τ > 0, i.e. the powers P 1 , P 2 remain strictly positive at all later times. The constraint 0 < v 0 < 1 implies that the initial speedθ 0 is restricted, since it follows from Eq. (23) that |θ | |a| + 2|b| + 1 at all times τ.

Properties of a, b
Of the two dimensionless parameters a, b, evidently b depends only on the optical fiber parameters, whereas a depends also on the total power P 0 , unless ∆β = 0. For the scalar model, when Eqs. (11) are satisfied, we have b = 1 but generally b = 1. In this case a set of steady state solutions appears (the states (24) discussed in Section 3.2 below) which for certain values of a, b are unstable. For fibers with elliptical cross sections we find that b > 1 and the unstable steady states exist provided 1 < a < 2b − 1. We have not, however, been able to eliminate the possibility that b < 1 for other geometries, and so in the following we also analyze the case b < 1. The parameter a can be positive or negative depending on the sign of ∆β and on the value P 0 ; when Eqs. (11) are satisfied we have a = −3∆β /(P 0 γ 1 ) + 1 and hence a can take large positive or negative values for small P 0 .
As an example, we have evaluated b using the definitions Eqs. (3)(4)(5) for step-index, air-clad glass waveguides with elliptical cross sections where the major/minor axes are denoted x, y. The host glass is taken to be chalcogenide with linear and nonlinear refractive indices of n = 2.8 and n 2 = 1.1 × 10 −17 m 2 /W at λ = 1.55µm (as in [52]). Fig. 1(i) shows a contour plot of log 10 b as a function of x, y. We see, as expected, that b approaches 1 as the waveguide dimensions x, y increase towards the operating wavelength. For small core waveguides, however, we find b > 1 with values as large as b ≈ 200. The parameter a, on the other hand, depends on both the structure and the total input power P 0 . For low input powers, specifically for P 0 γ ′ c ≪ |∆β |, a can take large negative values (for ∆β > 0) or positive values (for ∆β < 0) as shown in Fig. 1(ii).
For large values of P 0 , however, a approaches the constant C = (γ 2 − γ c )/γ ′ c , whose contours for elliptical core waveguides are shown in Fig. 1(iii); most such waveguides have positive C values ranging up to 400, but some, those in the region on the left side of the white curve in Fig. 1(iii), have negative or small values of C. The contour plot for ∆β in Fig. 1(iv) shows that ∆β takes a wide range of positive and negative values as x, y vary.

Steady state solutions
There are four classes of steady state solutions of Eqs. (22,23), each of which exist only for values of a, b within certain limits, as follows: provided b = 1 and 0 < a−1 2(b−1) < 1; provided b = −1 and 0 < a+1 2(b+1) < 1; provided |a| 1; and provided |a − 2b| 1. Of these four classes, (26) and (27) lie on the boundary of the physical region 0 < v < 1, but nevertheless influence properties of nearby nontrivial trajectories, and also play a role in soliton solutions. The states (24) lie within the physical region only if the parameters (a, b) belong to either the red or green region of the a, b plane shown in Fig. 2 (i). Similarly the solutions (25) satisfy 0 < v < 1 only in the disjoint regions of the a, b plane defined by either 2b + 1 < a < −1 or −1 < a < 2b + 1. If a, b lie outside these regions, and also outside the strips given by |a| 1 and |a − 2b| 1, there are no steady state solutions.
For special values of a, b these steady states can coincide, for example if a = 1 the solution (26) coincides with the boundary value of (24). Steady states for values of a, b on the boundary of the regions shown in Fig. 2 may need to be considered separately; for example if a = b = 1 then all steady states are given either by (25), or else by cos θ = 1 and any constant v.
In practice, the values of a, b are determined by the waveguide structure, the propagating mode and, in the case of a, the input power P 0 , and hence only restricted regions of the a, b plane are generally accessible. For example, Fig. 1(i) shows that for the fundamental mode of elliptical core fibers we have log 10 b 0, and so the attainable values of b are limited to b 1. We nevertheless include the case b < 1 in our analysis, since this possibility cannot be excluded for other fiber geometries. We discuss the accessible regions for the case of unstable steady states in Section 3.4.  (24), either 1 < a < 2b − 1 (red), or 2b − 1 < a < 1 (green); (ii) the regions of existence for the unstable solutions consisting of (24) (red), and (25) for which 2b + 1 < a < −1 (orange), together with (26,27) for which |a| < 1 or |a − 2b| < 1 (light blue).

Lagrangian formulation
We wish to determine the stability properties of each of the four classes of steady state solutions, in particular we look for unstable steady states. These are of interest because polarization states which lie close to these unstable states are very sensitive to small changes in parameters such as the total power P 0 , and so can flip abruptly as a function of the optical fiber length z. Although we may determine stability properties by investigating perturbations about the constant solutions, we find it convenient to reformulate the defining equations (22,23) as the Euler-Lagrange equations of a Lagrangian L which is a function of θ ,θ , and depends otherwise only on the parameters a, b. This also provides insight into the properties and solutions of these equations, and we may then investigate stability by examining the corresponding potential function. From and by substitution into (22) we obtain We consider Lagrangians L of the form where T is the (positive) kinetic energy, V is the potential energy, and the "mass" M is a positive function of θ . The equation of motion is and is identical to Eq. (29) provided We may therefore investigate all possible solutions θ (τ) by analyzing properties of the periodic potential V (θ ); every solution of the system of equations (22,23) corresponds to the trajectory θ (τ) of a particle of variable mass M in the potential V . Steady state solutions are zeroes of V ′ (θ ), and stability is determined by whether these zeroes are local maximums or minimums of V , subject to the constraint that the associated function v should always satisfy 0 < v < 1. Trajectories which begin near a local minimum, with a small initial speedθ (0), oscillate periodically with a small amplitude. On the other hand, trajectories which begin near an unstable point, i.e. near a local maximum of V , can display periodic oscillations of large amplitude with abrupt transitions between adjacent local maximums; we refer to these as switching solutions (previously bistable solutions [50]) since cos∆φ = cos(θ /2) switches periodically between two distinct values. Soliton trajectories also occur in which the particle moves between adjacent local maximums of V , see for example the discussion in [53], Section 2 and [54] for properties of solitons in optical fibers. As mentioned in Section 3.5, soliton trajectories also appear as the separatrix in phase plane plots.
We plot V as a function of θ and either a or b in Fig. 3, showing that V defines a complex surface with valleys and peaks which change suddenly as a or b are varied. Periodic solutions occur for trajectories restricted to a local valley, but there are also unbounded trajectories, in which θ increases or decreases indefinitely, depending on a, b and on whetherθ (0) is sufficiently large. The potential, as a function of θ and a, has saddle points which indicate that a stable solution can become unstable as a is varied; according to the definition (21) we may vary a within certain limits by varying the total power P 0 .
For a = b the potential is essentially that of the nonlinear pendulum under the influence of gravity, namely a simple cosine potential, but with a mass that depends on θ . Provided b > 1 this mass varies between two positive, finite limits. The unstable steady states correspond to a pendulum balanced upright, while the switching states (discussed in Section 4) correspond to trajectories which begin with the pendulum positioned near the top, possibly with a small initial speed, then swinging rapidly through θ = 2π to reach the adjacent unstable steady state. During this motion cos∆φ = cos θ 2 flips rapidly between the values ±1. The soliton discussed in the Appendix is the trajectory in which the pendulum begins at the unstable upright position Fig. 3. The potential V plotted as a function of (i) θ , a for b = 0.8; (ii) θ , b for a = 0. and, over an infinite time, moves through the stable minimum to the adjoining unstable steady state.
Although both M and V are singular when cos θ = b, which occurs only if |b| 1, this singularity is an artifact of the Lagrangian formulation, as is evident from Eqs. (22,23), which have smooth bounded right hand sides for any b. In particular v, which is obtained from Eq. (28) given θ , is a smooth function of τ even if cos θ = b for some τ.
The energy where c is the constant of integration. This constant is determined by first choosing initial values v 0 , θ 0 , where 0 < v 0 < 1, and then findingθ (0) from Eq. (23) which, from (33) evaluated at τ = 0, fixes c. We may integrate (33) to determine θ as an explicit function of τ, expressible in terms of elliptic functions, as discussed further in Section 3.5. A limitation of the Lagrangian formulation is that the constraint 0 < v < 1 is not easily implemented. Whereas every solution of the system (22,23) defines a trajectory θ (τ) in the Lagrangian system (29), the converse is not true, i.e not all trajectories in this system satisfy 0 < v < 1. The initial speedθ (0) must be restricted to only those values allowed by Eq. (23), in which 0 < v(0) < 1, and similarly the constant solutions of Eq. (29) are valid steady states for the original equations (22,23) only in certain regions of the a, b plane. Trajectories which violate 0 < v < 1, while not physical in the context of optical fiber configurations, can nevertheless be viewed as acceptable motions of the mechanical system defined by the Lagrangian (30). We investigate an alternative Hamiltonian formulation in terms of v in Section 3.5.

Stability of steady state solutions
The stability of each of the four classes of steady state solutions in Section 3.2 is determined by the sign of V ′′ for that solution; a positive sign implies that the solution lies at a local minimum of V and is therefore stable, whereas a negative sign implies that the solution is unstable.
For the remaining steady states (26,27), for which v = 0 or v = 1, we have V ′′ = −2 sin 2 θ /|a − b| which in all cases is negative, and so these states are unstable whenever they exist. This is consistent with the observation that v(τ) cannot attain the values 0, 1 at any time τ, provided 0 < v 0 < 1. The regions in the a, b plane where the unstable states exist are shown in Fig. 2 (ii).
Next, we determine conditions under which the unstable steady state solutions (24) are accessible. For elliptical core step index fibers, for which b > 1 as shown in Fig. 1(i), the region of instability is indeed accessible and leads to properties such as nonlinear self-polarization flipping, discussed in Section 4. The region of unstable solutions is given by 1 < a < 2b − 1, These inequalities specify the possible values, if any, of P 0 for which the unstable solutions exist for a fixed fiber. In order to visualize this region we plot a as a function of P 0 in Fig. 4(i), where a is given by (21). The boundaries of the unstable region at a = 1, a = 2b − 1 are shown by the green solid lines. First we consider fibers for which 1 < C < 2b − 1, where C = (γ 2 − γ c )/γ ′ c takes the value shown by the dashed line in Fig. 4(i). Then a has two branches associated with either ∆β < 0 or ∆β > 0; for the branch corresponding to ∆β < 0 (the solid blue line), a is large and positive for small P 0 and asymptotically approaches C for large P 0 . The intersection of this branch with the boundary a = 2b − 1 determines the minimum power P min 1 required in order to access the unstable region. In this case, only part of the unstable region corresponding to C < a < 2b − 1 is accessible, as shown by the blue region. For the ∆β > 0 branch (red solid curve) a is large and negative for small P 0 and asymptotically approaches C for large P 0 . For this branch, P 0 needs to be larger than a value P min 2 . The unstable region is accessible provided 1 < a < C and is a subset (red shaded) of the whole unstable solution region. Fig. 4(i) allows one to determine the minimum and maximum values of a and the minimum power to access the unstable solution region, once ∆β and C are known. For elliptical core fibers these two values are completely determined by the dimensions x, y, see Fig. 1(iii,iv) for plots of C and ∆β .
Besides fibers for which 1 < C < 2b − 1, there are the possibilities C > 2b − 1 or C < 1. From Fig. 1(iii,iv) one can show that these combinations (with ∆β positive or negative) either do not exist, or do not lead to unstable solutions, since the possible values of a do not lie in the unstable region 1 < a < 2b − 1. In summary, the only elliptical core fibers that allow unstable solutions are those with 1 < C < 2b − 1 with either positive or negative ∆β . The case in which ∆β = 0 is discussed separately in [55, 56].
Based on the above discussion, one can find the minimum power P min 0 required to generate unstable solutions for elliptical core fibers. Fig. 4(ii) plots log 10 (P min 0 ) (where P min 0 is measured in watts) as a function of x, y, where the white region corresponds to fibers for which there are no unstable solutions, and the regions below and above the diagonal line correspond to P min 1 and P min 2 , respectively, which have been obtained for the two branches of the function a(P 0 ) shown in Fig. 4(i).

Hamiltonian function
Although the Lagrangian formulation in terms of θ is convenient for an analysis of the steady states and their stability, and also for a qualitative understanding of all solutions including solitons, the constraint 0 < v < 1 is more easily implemented by means of a direct formulation in terms of v. This automatically eliminates unphysical trajectories for which one of the input powers P 1 , P 2 is negative. Such a formulation follows by construction of a Hamiltonian function which, being conserved, allows us to firstly integrate the nonlinear equations and obtain analyt- Hence as a function of τ, H is conserved and takes the constant value H 0 = H(v 0 , θ 0 ) on any trajectory. We may investigate all possible solutions, therefore, by analyzing the curves of constant H 0 in the v, θ plane. We have and from Eq.
where Q is the polynomial of 4th degree (provided b 2 = 1) given by Since the left hand side of (37) is positive, solutions exist only if sin 2 θ 0 0 (as follows from Eq. (22)) Q has at least two real zeroes, possibly repeated, and so there is an interval within 0 < v < 1 in which Q(v) > 0, and so solutions always exist. If the initial values v 0 , θ 0 are such that the trajectory begins in a stable steady state, v remains constant for all τ > 0, otherwise the trajectory is nontrivial. There are two types of nontrivial solutions, periodic and soliton solutions.
We can gain insight into possible solutions by plotting contours of constant H(v, θ ) in the v, θ plane, which supplies essentially a phase portrait of the system. Solutions for which both v, θ are periodic in τ form closed loops, and lie close to a stable steady state, whereas nonperiodic trajectories lie outside the separatrix which defines soliton solutions, as we discuss in the Appendix. Fig. 5 shows two examples in which stable steady states are marked in green, and unstable steady states are shown in red or orange. Periodic solutions are evident as closed loops surrounding stable steady states, whereas the separatrix marks soliton trajectories which connect unstable steady states. Apart from these solitons, all other solutions v, cosθ (but not necessarily θ ) are periodic in τ. The switching solutions of particular interest, in which the state of polarization inside the waveguide flips between two well-defined states, are those close to the separatrix.

Periodic solutions
Periodic solutions v of (37) attain both minimum and maximum values, denoted v min , v max respectively, with 0 < v min v max < 1. Sincev = 0 at a maximum or minimum of v, both v min , v max are roots of Q. We can factorize Q as a product of quadratic polynomials, and hence explicitly find all roots, and so identify v max and v min . We integratev = Q(v) over the half-period in which v increases, in order to find τ as a function of v, and also the period T : where τ 0 is the time at which v achieves its minimum, i.e. v min = v(τ 0 ). These integrals may be evaluated in terms of elliptic integrals of the first kind, see for example the explicit formulas in [57] (Sections 3.145, 3.147). In particular, T is expressible in terms of the complete elliptic integral K, and so can be written as an explicit function of a, b, v 0 , θ 0 , i.e. as a function of the waveguide parameters and the initial power and phase of the input fields. The precise formulas depend on the relative location of the roots of Q.
Having found v, cos θ is obtained from Eq. (36) and is also periodic in τ, as isθ which is obtained from Eq. (23), however θ itself need not be periodic. Although it is straightforward to find v, θ numerically as functions of τ, for specified numerical values of a, b and initial values v 0 , θ 0 , the exact solutions are useful because they display the exact dependence of the solution on all parameters, such as the total power P 0 ; it is not necessary therefore to solve the equations numerically for every choice of P 0 , rather the exact solution gives the explicit periodic solution and the period as known functions of P 0 .
For switching solutions, the phase difference between the two polarization vectors experiences abrupt phase shifts through π as the light propagates within the waveguide. As a result, the state of polarization flips between two well-defined polarization states, where the flipping angle depends on a, b and on θ 0 , v 0 . The following are two examples of switching solutions.
As the first example we choose a = 1, b = 4 with the initial values v 0 = ε, θ 0 = 0, where ε = 10 −4 , in which case the input laser beam is linearly polarized and the polarization state is close to one of the principle axes of the waveguide. Hence, the trajectory starts near the unstable steady states (24) or (26), which lie on the boundary of the red region shown in Fig. 2 (i). We plot v and cos θ 2 = cos ∆φ as a function of τ in Fig. 6 (i), showing switching behavior for cos θ 2 , which is periodic and flips abruptly between the values ±1; θ , however, is an increasing function of τ, with jumps through 2π at periodic intervals. The polarization vector experiences an angular flipping associated with the abrupt flipping of cos ∆φ , however, since v 0 = ε and θ 0 = 0, the flipping angle is very small, as depicted in the inset of Fig. 6 (i). Regarded as the trajectory of a particle of mass M in the potential V in Eq. (32) this motion corresponds to a particle moving slowly over the peaks of the potential, which are the unstable steady states, then sliding quickly down the valleys through the minimum values of V and back to the peaks. For a = 1 the potential is flat at its maximum values, since in this case V ′ = 0 = V ′′ = V ′′′ , hence v,θ are each close to zero except when θ moves to an adjoining maximum of V . In terms of the contour plots shown in Fig. 5(i) this trajectory corresponds to the contour which begins just above the unstable steady state (orange dot) and closely follows the separatrix shown in red (which is the soliton solution discussed in the Appendix) with a maximum value ∼ 0.4 for v. As a second example of switching behavior we choose a = b = 2 with v 0 = 1/2, θ 0 = ε, where ε = 10 −4 , which corresponds to a linearly polarized input laser beam in which the polarization vector makes an angle of 45 • to either of the principle axes of the waveguide. Again, the initial value lies close to an unstable steady state (24) and a, b lie within the red region of instability in Fig. 2. We plot v and cos(θ /2) as functions of τ in Fig. 6 (ii), showing the periodicity of these functions and the switching behavior of cos(θ /2). Since v 0 = 1/2, the angular flipping of the polarization vector is π/2, because cos(θ /2) flips between values ±1, as shown in the inset of Fig. 6 (ii). Unlike the previous example, θ is also periodic in τ with a trajectory that corresponds to the motion of a particle in the potential V , starting slowly near the unstable steady state (24) but sliding rapidly through the potential minimum to approach an adjoining unstable steady state. This motion is similar to the periodic oscillations of a nonlinear pendulum (since a = b, see the definition of V in Eq. (32)) with a large amplitude of almost 2π, and v attains nearly all values between 0, 1. In terms of the phase space contours shown in Fig. 5(ii), the motion corresponds to a periodic trajectory which begins near the red dot (unstable steady state) and again closely follows the separatrix which marks the soliton trajectory.

Discussion and Conclusion
Switching states, as defined and demonstrated here through simulation by means of a full vectorial model, are attractive for practical applications, since they allow nonlinear self-flipping of the polarization states of light propagating in an optical waveguide. This flipping is due to the nonlinear interactions of the two polarizations, and has properties that depend on the total optical power and on the specific fiber parameters. These properties can in principle be employed to construct devices such as optical logic gates [58], fast optical switches and optical limiters [55,56], in which small controlled changes in the input parameters lead to sudden changes in the polarization states.
The minimum power necessary to generate such switching states is determined for any waveguide by the inequalities (34) and, for chalcogenide optical nanowires with elliptical core cross sections, is summarized in Fig. 4 (ii). The minimum power required in such nanowires is in the range 1 − 10kW which, although not practicable for CW lasers, can be achieved in pulsed lasers. Although we have limited our analysis to the static case, ignoring the temporal variation of laser light, it is still applicable to slow pulses with pulse widths in the order of nanoseconds depending on the dispersion of the waveguide. A more practical minimum power requirement that achieves switching behavior is by means of asymmetric waveguides, such as rib waveguides, for which ∆β can be reduced to very small values while still having different field distributions for the two polarizations, as discussed in [55,56].
The nonlinear interactions of the two polarizations can be impacted by two factors that have not yet been investigated: (1) interactions with higher order modes in few-moded waveguides and, (2) contributions from nonlinear terms containing different forms of e 1 e 2 , i.e., nonzero values for the coefficients γ (1) µν , γ (2) µν , γ ν in Eqs. (6)(7)(8). (This applies only when e 1 e 2 is no longer approximately zero, as assumed in this paper). In few-moded waveguides, higher order modes contribute to the nonlinear phase of each polarization of the fundamental mode through cross phase and coherent mixing terms. Inspection of Eq. (2) reveals that nonzero γ (1) µν , γ (2) µν , γ (3) µν coefficients significantly change the dynamics of nonlinear interactions of the two polarizations and most likely lead to different parameter regimes for the existence of periodic and solitonic solutions. These factors will be the subject of further studies.
In summary, we have developed the theory of nonlinear interactions of the two polarizations using a full vectorial model of pulse propagation in high index subwavelength waveguides. This theory indicates that there is an anisotropy in the nonlinear interactions of the two polarizations that originates solely from the waveguide structure. We have found all static solutions of the nonlinear system of equations by finding exact constants of integration, which leads to expressions for the general solution in terms of elliptic functions. We have analyzed the stability of the steady state solutions by means of a Lagrangian formalism, and have shown that there exist periodic switching solutions, related to a class of unstable steady states, for which there is an abrupt flipping of the polarization states through an angle determined by the structural parameters of the waveguide and the parameters of the input laser. By means of a Hamiltonian formalism we have analyzed all solutions, including solitons which we have shown are close to the switching solutions of interest.