Analytic Methods to Find Beating Transitions of Asymmetric Gaussian Beams in GNLS equations
David Ianetz Jeremy Schiff

The Price of Uncertainty in Present-Biased PlanningWork supported by the European Research Council, Grant Agreement No. 691672.
    [
    
===============================================================================================================================


  In a simple model of propagation of asymmetric Gaussian beams in nonlinear waveguides,
  described by a reduction to ordinary differential eqautions of
  generalized nonlinear Schrödinger equations (GNLSEs)   
  with cubic-quintic (CQ) and saturable (SAT) nonlinearities and a graded-index profile,
  the beam widths exhibit two different types of beating behavior, with transitions between them.
  We present an analytic model to explain these phenomena, which originate in a 1:1 resonance
  in a 2 degree-of-freedom Hamiltonian system. We show how small oscillations near a fixed
  point close to 1:1 resonance in such a system can be approximated using an integrable
  Hamiltonian and, ultimately, by a single first order differential equation. In particular, the
  beating transitions can be located from coincidences of roots of a pair of
  quadratic equations, with coefficients determined (in a highly complex manner) by the
  internal parameters and initial conditions of the original system. The results of the analytic
  model agree with numerics of the original system over large  parameter ranges, and allow
  new predictions that can be verified directly. In the CQ case we identify a band of beam energies
  for which there is only a single beating transition (as opposed to 0 or 2) as the eccentricity is increased.  
  In the SAT case we explain the sudden (dis)appearance of beating transitions for certain
  values of the other parameters as the grade-index is changed.  








§ INTRODUCTION

In the sequence of papers <cit.>
a variational approach was taken to investigate the propagation of
asymmetric (elliptic) Gaussian beams in nonlinear waveguides, with cubic-quintic and saturable nonlinearities and
a parabolic graded-index (GRIN) profile, as described by suitable 
generalized nonlinear Schrödinger equations (GNLSEs). 
The beam widths in the two transverse directions to the direction
of propagation were found to obey a set of ordinary differential equations which can be identified as 
the equations of motion of a point particle in certain rather complicated, but
tractable, 2d potentials. 
Numerical analysis of these equations revealed “beating” phenomena: in addition to fast oscillations, the beam
widths exhibit a (relatively) slow periodic variation. Furthermore, two types of beating were identified: In type I beating
the amplitude of oscillation of the beam width in one direction remains greater than the amplitude of oscillation in the
other direction, whereas in type II, there is an interchange between the widths in the two transverse directions. The type
of beating depends on the parameters of the system and initial eccentricity of the beam. Remarkably, as the initial
eccentricty or other parameters are changed, there can be a transition between types, and this transition
is characterized by a singularity in the ratio of the periods of the beating and of the fast oscillatory motion. 

The intention of the current paper is to provide a theoretical analysis of the beating phenomena and, in particular,
to present an approximate analytic method to find the transitions between types. The relevant tool is the analysis
of small oscillations in
2 degree-of-freedom Hamiltonian systems near a fixed point which is close to 1:1 resonance. The fact that resonance
is the source of “beating” or “energy transfer” phenomena in mechanical systems is well known. A classic example
can be found in the paper of Breitenburger and Mueller <cit.> on the elastic pendulum, which the authors describe as
a “paradigm of a conservative, autoparametric system with an internal resonance”. The paper <cit.> has other features in
common with our work (such as the use of action-angle variables and the fact that the analytic approximation used is 
a single elliptic function equation) but it is in the much simpler context of 1:2 resonance. For other examples
of autoparametric resonance see, for example, <cit.>. The most widely used tool for analysis of
systems near resonance is
the mutliple time scale method, see for example <cit.> for  thorough presentations and many examples. For a typical
modern application see <cit.>. However, averaging techniques present an alternative <cit.>, and in
the context of Hamiltonian systems, working in action-angle coordinates has substantial advantages <cit.>. A typical
study of a system near resonance will involve looking at the bifurcations of special solutions. In this context much attention
has been paid to the definition and identification of nonlinear normal modes — see <cit.> for a review, and
<cit.> for an example in the context of 1:1 resonance. 

The 2 degree-of-freedom Hamiltonian systems we study have a discrete symmetry, and are approximated by a family of systems 
with  1:1 resonance studied  nearly 40 years ago by Verhulst <cit.>. 
Verhulst showed the existence of an approximate second integral and used this to study bifurcations of special solutions
and their stability. Our work differs from that of Verhulst and other works on 1:1 resonance in several regards. The 
bifurcation question we pose depends not only on the internal parameters of the system, but also on the initial conditions.
The question is not only one of identifying different types of solutions of the system, but also seeing how the type of solution changes
as both the initial condition and internal system parameters are varied. We have not seen a similar study in the
highly complex context of 1:1 resonance.  Our methodology uses 
action-angle variables and canonical transformations (though in an appendix we show how to apply standard two time scale
techniques). Unlike in most existing studies, it is necessary
to compute the relevant canonical transformation to second order.
However, this does not affect the result that once the correct canonical transformation has been applied,
the resulting approximating Hamiltonian depends only on a single combination of the angle variables
and is integrable. The equations of motion for the integrable Hamiltonian can be reduced to a single 
first order differential equation, and the rich bifurcation structure of the systems we study can
reduces to understanding the bifurcations of roots of a pair of quadratic polynomials,
with coefficients that depend (in a complex, nonexplicit manner) on the internal parameters
of the systems and the initial conditions. 
Comparison with numerical results shows our method gives high-quality results in a significant region of
parameter space, and allows a variety of interesting new predictions. 

The structure of this paper is as follows.
In the next section we review the relevant models from nonlinear optics and the
collective variable approximation to obtain equations for the propagation of beam widths,
and present the main findings of papers <cit.>  and some
further numerical results.  In section 3, we develop our 
method of integrable approximation for small oscillations in a 2 degree-of-freedom Hamiltonian system 
near a fixed point close to 1:1 resonance.
In section 4 we describe the application of this method to the specific systems relevant to beam propagation,
confirming existing numerical results and presenting new predictions. 
In section 5 we summarize and conclude. Appendix A completes some technical details omitted from the main text, 
and Appendix B describes an alternate
method of approximation of the full equations using a two time expansion. This is a more ad hoc approach than the one
explained in section 3, but we include it as it is more commonly used in the literature, 
and for certain values of parameters it gives better results. 

Before closing this introduction we mention a number of points concerning the relevance of the work
in this paper to optical solitons. We will describe in next section the manner in which we use ordinary differential
equations (ODEs) to study the behavior of solutions of GNLSEs. The use of ODEs to study GNLSEs is widespread, see for example 
   <cit.>
   
   In particular, the last two papers use ODE methods in the study of rotating solitons.
Our work extends the catalog of interesting bifurcations that can be observed in the context of GNLSEs; 
for another example; see the papers <cit.> for a case of a saddle-loop bifurcation. 
Finally, we mention that we neglect dispersive terms in the GNLSEs we study. 
  This is justifiable in the context of  new optical materials
  <cit.>
characterized by Kerr coefficients of the order 10^-11–10^-12 cm^2/ W, making
the critical intensity for self-focusing  small enough that it can be reached
using microsecond pulses and possibly even continuous wave (CW) laser beams. 



§ MODELS, THE COLLECTIVE VARIABLE APPROACH AND NUMERICAL RESULTS


We consider beam propagation in a nonlinear, graded-index fiber, as described by one of the following GNLSEs:

    2iψ_z + ψ_xx + ψ_yy + (  |ψ|^2 - Q|ψ|^4 - g(x^2+y^2) ) ψ   =    0 ,
         
      
     2iψ_z + ψ_xx + ψ_yy + (  |ψ|^2/1+α^2|ψ|^2  - g(x^2+y^2) ) ψ   =    0 .

Here, modulo suitable normalizations <cit.>,
ψ is the strength of the electric field, z is the
longitdinal coordinate, x,y are transverse coordinates, and Q,α,g are parameters. 
The first equation is the case of cubic-quintic nonlinearity (CQ), the second is the case of saturable nonlinearity (SAT).
In the low intensity limit these models are similar, but for higher intensity they display different physical properties. 
In both cases, the higher order nonlinearity prevents beam collapse associated
with the standard Kerr nonlinearity  <cit.>.
The term -g(x^2+y^2)ψ reflects the graded-index nature of the fibre,
that the refractive index n falls
with distance r from the center of the fibre according to the law n^2 = n_0^2 - Gr^2; the physical significance of this 
is explained in <cit.>.   

The collective variable approximation (CVA), introduced for the study of self-focusing beams in
<cit.>, is a variational technique
to approximate solutions of nonlinear Schrödinger-type equations which has been
used and validated in many different situations <cit.>. 




The method replaces partial differential equations such as (<ref>) and (<ref>) by
a system of ordinary differential equations
for the coefficients of an ansatz for the full solution.
The GNLSEs (<ref>) and (<ref>) are variational equations for action principles based on
the Lagrangian densities

    L_ CQ   =    i( ψψ^*_z-ψ^*ψ_z ) + |ψ_x|^2 + |ψ_y|^2 
           - 1/2 |ψ|^4  + Q/3|ψ|^6   
               + g(x^2+y^2) |ψ|^2  ,  
     L_ SAT   =     i( ψψ^*_z-ψ^*ψ_z ) + |ψ_x|^2 + |ψ_y|^2 
      + ln(1 + α^2|ψ|^2)-α^2|ψ|^2/α^4  
      + g(x^2+y^2) |ψ|^2 .

We assume ψ takes the form of the trial function 

    ψ_T(x,y,z) = A(z)
      exp(iϕ(z)
      -x^2/2a_x^2(z) + ib_x(z)x^2
      -y^2/2a_y^2(z) + ib_y(z)y^2
      ) ,

where A,ϕ,a_x,a_y,b_x,b_y are currently undetermined, real functions of only the longitudinal coordinate z.
This trial function describes an elliptic Gaussian
beam with a_x,a_y representing the widths of the beam in the x,y directions. b_x,b_y describe curvatures
of the beam wavefront, A is the normalized amplitude of the electric field, and ϕ is a longitudinal phase factor.
Our choice of a Gaussian shape for the trial function is appropriate because the Gaussian is an exact
solution of the linear Schrödinger equation for GRIN waveguides <cit.>.
Substituting the trial function in the Lagrangian densities
(<ref>),(<ref>) and computing the integrals over the variables x,y we obtain reduced densities
for the functions A,ϕ,a_x,a_y,b_x,b_y. The corresponding Euler-Lagrange equations
in the CQ case are 

    Ȧ   =    -(b_x +b_y)A  , 
    ȧ_x,y   =    2a_x,yb_x,y ,  
    ḃ_x,y   =   1/2a^4_x,y - 2b^2_x,y - g/2
                    - A^2/a^2_x,y( 1/8- QA^2/9) ,  
    ϕ̇   =    - 1/2a^2_x - 1/2a^2_y + ( 3/8 - 5QA^2/18) A^2  
                      .

Here a dot denotes differentiation with respect to z. 
In the SAT case the equations for A,a_x,a_y remain the same, but those for b_x,b_y,ϕ  are replaced by 

    ḃ_x,y   =   1/2a_x,y^4 -2b_x,y^2 - g/2 +
                           ln(1+ α^2 A^2)+Li_2(-α^2 A^2)/2α^4 A^2a_x,y^2 , 
    ϕ̇   =     -1/2a_x^2 - 1/2a_y^2
                           + α^2A^2-2ln(1+α^2A^2)- Li_2(-α^2 A^2)/2α^4A^2 .

Here Li_2(x)=∑_k=1^∞x^k/k^2  is the Spence or dilogarithm function  <cit.>. 
For both CQ and SAT cases we observe that A^2 a_xa_y is conserved <cit.>,
and we write A^2 a_x a_y = 4 E (4E is the beam energy), and use this to
eliminate A(z) . Furthermore, ϕ(z) evidently plays no role in determining the other functions and
can be computed by a simple quadrature once the other functions have been found. 
Furthermore, it is clear that
we can write b_x (b_y) in terms of a_x  (a_y) and its z-derivative. Thus we can reduce the system of 6
equations to a pair of second order equations for a_x,a_y. After some more calculation it emerges that the
equations are simply the equations of motion 

    ä_x = -∂ V/∂ a_x ,     ä_y = -∂ V/∂ a_y
    
for a particle in a potential V(a_x,a_y), where for CQ

    V = V_CQ≡1/2( 1/a_x^2 + 1/a_y^2) 
    - E/a_x a_y
    + 16 QE^2/9a_x^2a_y^2 
    + g/2( a_x^2+a_y^2) ,

and for SAT 

    V = 
    V_SAT≡1/2( 1/a_x^2 + 1/a_y^2) 
    -a_xa_y/4Eα^4 Li_2 ( -4Eα^2/a_xa_y) 
    + g/2( a_x^2+a_y^2) .


Thus integration of equations (<ref>) for potentials (<ref>) and (<ref>) provides a
first approximation to solutions of the GNLSEs (<ref>) and (<ref>). Full numerical solutions
of  GNLSEs have been given in both 
the SAT <cit.> and the CQ <cit.> cases with g=0. 
In  <cit.> it was shown that the breathing frequencies found numerically are
similar to those obtained by the CVA technique.
In <cit.> it was shown that 
the shape of the beam obtained numerically for a saturable medium remains similar to  Gaussian, even
for an asymmetric initial condition. 
However, use of direct numeric methods
to give an overall picture of the behavior of a GNLSE, as a function of all the various parameters, 
remains a computationally overwhelming task, and
having an qualitatively correct analytic or semianalytic model is therefore
useful for developing physical insight
<cit.>. 

Appropriate initial conditions for (<ref>) are

    a_x(0) = a_0 r  ,       a_y(0) = a_0/r ,     ȧ_̇ẋ(0) = ȧ_̇ẏ(0) = 0 .

The latter two conditions are equivalent to taking b_x(0)=b_y(0)=0. 
Note that both the CQ and the SAT system have a scaling symmetry

    [ a_x →λ a_x  ,  a_y →λ a_y , a_0 →λ a_0  ,    r   → r  ,   z →λ^2 z  ,;  Q  →λ^2  Q ,        α→λα , g  →λ^-4  g ,   E  →   E  .               ]

Thus for CQ we do not need to study the dependence of solutions on the 5 parameters Q,g,E,a_0,r, but only on
the 4 scale invariant quantities Qa_0^-2, ga_0^4, E, r. On occasion we will work 
with the scale invariant quantity  K_ CQ = 4QEa_0^-2 instead of the quantity Qa_0^-2.
(For SAT, replace all instances of Q in the previous two sentences with α^2, and K_ SAT=4α^2Ea_0^-2.)
Note that since there is symmetry in both the models between a_x and a_y, there is a r→1/r inversion symmetry,
and thus we need only study r≤ 1 or r≥ 1. 
   
In the papers <cit.> the ODE systems above were studied numerically.
For appropriate choices of the parameters “beating” phenomena were observed: in addition to (relatively) fast
“breathing” oscillations, the beam
widths exhibit a (relatively) slow periodic variation. Two types of  beating were identified: 
In type I beating,
the amplitude of oscillation of the beam width in one direction remains greater than the amplitude of oscillation in the
other direction. In type II beating, there is an interchange between the widths in the two transverse directions. 
This is illustrated in Figure 1, which shows solutions of the CQ system
for E=2.039, K_ CQ=0.71, ga_0^4=0.01, and two choices of r: r=1.14 gives type I beating, whereas 
r=1.16 gives type II beating.





The type of beating depends on the parameters of the system and, as evident from Figure 1, on the  initial eccentricity of the beam.
Remarkably, as the initial eccentricty is increased, or as other parmeters are changed, there can be a transition between  
types. The approach to this transition is characterized by a divergence in the ratio of the periods of the slow beating and of
the fast oscillatory motion. In Figure 2 this ratio (determined from numerical simulations) is plotted as a
function of r^2 for the CQ system,
with E=2.039, K_ CQ=0.71 and  ga_0^4=0,0.01,0.02,0.03. (The reason for the choice of the coordinate r^2 on the x-axis is
simply to make the plot clearer.)
For r just above
1 the beating is type II, then there is a transition to type I, and then a second transition back to type II.
The dependence on the system parameters
of the two critical values of r, which we denote collectively by r_c, 
is explored further in Figure 3. 
In Figure 3a  the values of r_c are plotted as a function of ga_0^4 for three different values of
K_ CQ and a constant value of  E;
in Figure 3b r_c is plotted as a function of ga_0^4 for three different values of
E and a constant value of  K_ CQ.
In general we see that r_c increases as a function of ga_0^4 (for fixed E,K_ CQ). 
From Figure 3b we see that since the (solid) red is above the (dashed) blue is above the (dot-dashed) black,
r_c also increases as a function of E (for fixed ga_0^4,K_ CQ). But in Figure 3a we see there is
difference between the upper  and lower branches of r_c. We deduce that 
the higher value of r_c also increases with K_ CQ (for fixed E,ga_0^4), but the
lower value decreases. 











We shall see later that for other values of ga_0^4,K_ CQ,E there can be just a single transition or
no transitions at all as r is increased from 1.  Transitions between beating types are also
observed in the SAT system, again with a complex dependence on the parameters ga_0^4, K_ SAT and E. 
The aim of this paper is to provide an integrable approximation for 
equations (<ref>) with  potentials (<ref>) and (<ref>) which provides
a theoretical model to predict where the transitions between types take place.   



§ SMALL OSCILLATIONS NEAR 1:1 RESONANCE






In this section we describe a general process of approximation near a 1:1 resonance 
for a 2 degree-of-freedom Hamiltonian system with Hamiltonian  

    H =  1/2(p_x^2 + p_y^2 )  + V(a_x,a_y) .

Here a_x,a_y are the coordinates, p_x,p_y are the conjugate momenta, and the potential V (which typically will depend on
a number of parameters) is symmetric, V(a_x,a_y)=V(a_y,a_x).
We assume that for typical values of the parameters 
the potential has an isolated symmetric minimum (at a_x=a_y=a_ min, say) at which the system is close to 1:1 resonance.
Note that because of the symmetry, ∂^2V/∂ a_x^2 (a_ min, a_ min)
  = ∂^2V/∂ a_y^2 (a_ min,a_ min). 
Thus the Hessian matrix of the potential at (a_ min,a_ min) has
eigenvectors ( [   1; ± 1 ])
with eigenvalues ∂^2V/∂ a_x^2 (a_ min,a_ min)
  ±∂^2V/∂ a_xa_y (a_ min,a_ min).
The condition for being close to 1:1 resonance (i.e. equal eigenvalues) is therefore simply
∂^2V/∂ a_xa_y (a_ min,a_ min) ≈ 0. 
For these systems we study orbits with initial conditions as given in (<ref>).



The process of approximating such a system with an integrable  system has 3 steps. 

The first step is to expand in normal coordinates near the fixed point, retaining only
terms up to order 4 in the potential. Thus we write

    a_x = a_ min + ζ_2+ζ_1/√(2) ,        a_y = a_ min + ζ_2-ζ_1/√(2)

and expand to fourth order to obtain

    H_1 = 1/2( p_1^2 + p_2^2 + ω_1^2 ζ_1^2 + ω_2^2 ζ_2^2 )  + a_1 ζ_1^2 ζ_2 + a_2 ζ_2^3
    + a_3 ζ_1^4 + a_4 ζ_1^2 ζ_2^2 + a_5 ζ_2^4

where p_1,p_2 are the conjugate momenta to the coordinates ζ_1,ζ_2, and ω_1,ω_2,a_1,a_2,a_3,a_4,a_5
are constants that depend on the parameters of the original potential V. The Hamiltonian H_1 has ζ_1 → -ζ_1
symmetry as a consequence of the symmetry of H, and is the general Hamiltonian with this symmetry and a quartic potential. 
Aspects of the behavior of this Hamiltonian at, or close to, 1:1 resonance have been studied previously, for example,  in 
<cit.>. The initial conditions
for this system, corresponding to (<ref>) are 

    ζ_1(0) = a_0/√(2)(r-1/r)
       ,     ζ_2(0) = 1/√(2)(a_0(r+1/r) - 2 a_ min)
       ,     p_1(0) = p_2(0) = 0  .
 
Symmetric solutions correspond to the initial condition ζ_1(0)=0. In regarding H_1 as an approximation for H
we are neglecting terms of fifth order and above.  

The second step is to make the canonical
transformation to action-angle coordinates associated with the quadratic part of the Hamiltonian H_1,
i.e. to substitute

    [      ζ_1 = √(2J_1/ω_1)cosθ_1 ,    p_1 = - √(2J_1ω_1)sinθ_1  ,;      ζ_2 = √(2J_2/ω_2)cosθ_2 , p_2 = - √(2J_2ω_2)sinθ_2     . ]

This gives 

    H_2    =   ω_1 J_1 + ω_2 J_2  
      +   ( 2J_1a_1/ω_1 +3 J_2 a_2/ω_2) √(J_2/2ω_2)cosθ_2  
           
      + (J_2/2ω_2)^3/2  2a_2cos 3θ_2
      +  a_1 J_1/ω_1√(J_2/2ω_2)( cos(2θ_1-θ_2)  + cos(2θ_1+θ_2) )    
           +   3a_3 J_1^2/2ω_1^2
         + a_4 J_1J_2/ω_1ω_2
         + 3a_5 J_2^2/2ω_2^2
       + a_4 J_1J_2/2ω_1ω_2( cos(2θ_1-2θ_2)  + cos(2θ_1+2θ_2) )
         
          ( 2a_3 J_1^2/ω_1^2 +  a_4 J_1J_2/ω_1ω_2)cos(2θ_1)
            + a_3 J_1^2/2ω_1^2cos(4θ_1)
            +(a_4J_1J_2/ω_1ω_2+ 2a_5J_2^2/ω_2^2)cos(2θ_2)   
                   +a_5J_2^2/2ω_2^2cos(4θ_2)    .

Here θ_1,θ_2 are the angle variables, and J_1,J_2 the conjugate actions. 
The initial conditions for the action variables are

    J_1(0) = a_0^2ω_1/4(r-1/r)^2  ,     
    J_2(0) =   ω_2/4(a_0(r+1/r) - 2 a_ min)^2 .

The initial conditions for the angle variables depend on the sign of ζ_1(0) and ζ_2(0). 
If ζ_1(0)>0 (ζ_2(0)>0) then, from (<ref>) we should take θ_1(0)=0 (θ_2(0)=0)
and otherwise θ_1(0)=π (θ_2(0)=π). Due to the ζ_1→ -ζ_1 symmetry of H_1 
the Hamiltonian H_2 has period π (and not 2π) as a function of θ_1 and thus the choice of
the θ_1 initial condition is irrelevant. The choice of the θ_2 initial condition, however, is
important. We are introducing a non-physical discontinuity in the approximation procedure when
the sign of ζ_2(0) changes, i.e. when  r+1/r=2a_ min/a_0. We will see the effects of this
later, in our results for the SAT potential. 

The third step involves a canonical change of coordinates (θ_1,θ_2,J_1,J_2)→ (ϕ_1,ϕ_2,K_1,K_2)
defined by a generating function of the second type G_2(θ_1,θ_2,K_1,K_2) <cit.>, chosen to eliminate the
nonresonant terms from the Hamiltonian (i.e. all the trigonometric terms of order ||J||^3/2 or ||J||^2 except the one involving
cos(2θ_1-2θ_2).) 
The full change of coordinates is given by 

    [ ϕ_1 = ∂ G_2/∂ K_1 ,  J_1 = ∂ G_2/∂θ_1 ,; ϕ_2 = ∂ G_2/∂ K_2 ,  J_2 = ∂ G_2/∂θ_2 . ]

The generating function G_2 should be taken in the form

    G_2    =     K_1θ_1 + K_2θ_2 +  A_1 sinθ_2 + A_2sin 3θ_2 +  A_3 sin(2θ_1-θ_2)
            + A_4 sin(2θ_1+θ_2)  
            +  A_5 sin2θ_1 + A_6 sin 4θ_1  + A_7 sin2θ_2 + A_8 sin 4θ_2
            + A_9 sin6θ_2  + A_10sin (2θ_1 +2θ_2)
            
            + A_11sin (2θ_1 +4θ_2)
            + A_12sin (4θ_1 +2θ_2)
            + A_13sin (2θ_1 -4θ_2) 
            + A_14sin (4θ_1 -2θ_2)

where the coefficients A_1,…,A_14 are functions of K_1,K_2, which  are chosen to eliminate the nonresonant
trigonometric terms in the Hamiltonian to required order.  A_1,A_2,A_3,A_4 are of order ||K||^3/2 and A_5,…,A_14
are of order ||K||^2.  The calculations are long, but straightforward with the help of a symbolic manipulator, and the final
Hamiltonian is found to be simply 

    H_3 =  ω_1 K_1 + ω_2 K_2 
    +  b_1 K_1^2 + b_2 K_1K_2 + b_3 K_2^2 
    + (b_4 K_1^2 + b_5 K_1K_2) cos( 2(ϕ_1-ϕ_2)   )

where 

    b_1    =   3a_3/2ω_1^2 
      -   a_1^2( 8ω_1^2-3ω_2^2)/4ω_1^2ω_2^2( 2ω_1-ω_2 )  ( 2ω_1 + ω_2 )
     
    b_2    =   a_4/ω_1ω_2 - 3 a_1 a_2/ω_1ω_2^3
    - 2a_1^2/ω_1ω_2( 2ω_1-ω_2)  ( 2ω_1+ω_2) 
    
      b_3    =   3a_5/2ω_2^2 -  15 a_2^2/4ω_2^4
    
    b_4    =   ( ω_2-ω_1) a_1^2/2ω_1^2ω_2^2( 2ω_1-ω_2) 
              
    b_5    =   a_4/2ω_1ω_2
    - a_1a_2(4ω_1^2-3ω_1ω_2-4ω_2^2)/2ω_1ω_2^3( 2ω_1-ω_2)( 2ω_1+ω_2)
       - a_1^2/ω_1^2ω_2(2ω_1-ω_2) .


The Hamiltonian H_3 given in (<ref>) is an integrable approximation of the original Hamiltonian H
  given in (<ref>).  H_3 is a normal form for the “natural” Hamiltonian H_1 at or near 1:1 resonance.
Note that in the case of exact resonance  ω_1=ω_2  the coefficient b_4 vanishes. Also close to resonance,  
the corresponding term in H_3 is of lower order than the other terms, and in  <cit.> it is omitted. 
However, we choose to  retain it to avoid any assumption on the relative orders of magnitude of |ω_1 - ω_2|
and ||K||.  The integrability of H_3 is evident, as it only depends on the modified angle variables
ϕ_1,ϕ_2 through the combination ϕ_1-ϕ_2. As a consequence the quantity K_1+K_2 is conserved,
in addition to the Hamiltonian itself. We denote the value of the Hamiltonian by E and the value of K_1+K_2
by P (these should be computed from the system parameters and initial conditions). The full equations of motion are 

    ϕ̇_1 = ∂ H_3/∂ K_1   =   ω_1 + 2 b_1 K_1 + b_2 K_2 +  (2 b_4 K_1 + b_5 K_2) cos(2(ϕ_1-ϕ_2))    , 
    ϕ̇_2 = ∂ H_3/∂ K_2   =   ω_2 +  b_2 K_1 + 2 b_3 K_2 + b_5 K_1 cos(2(ϕ_1-ϕ_2))   ,  
    K̇_1 = -∂ H_3/∂ϕ_1   =    2 K_1 (b_4 K_1 + b_5 K_2) sin (2(ϕ_1-ϕ_2))    ,  
    K̇_2 = -∂ H_3/∂ϕ_2   =    -2 K_1 (b_4 K_1 + b_5 K_2) sin (2(ϕ_1-ϕ_2))  .

Using the two conservation laws it is possible to eliminate K_2 and ϕ_1-ϕ_2 from the K_1 equation of motion
to get a single equation for K_1:

    K̇_1^2     =    - 4 
      ((b_1-b_2+b_3+b_4-b_5)K_1^2+((b_2-2b_3+b_5)P-ω_2+ω_1)K_1+b_3P^2+ω_2P- E)  
          
      ((b_1-b_2+b_3-b_4+b_5)K_1^2+((b_2-2b_3-b_5)P-ω_2+ω_1)K_1+b_3P^2+ω_2P- E) .










Equation (<ref>) is a central result of this paper. 
To solve (<ref>) it is necessary to translate the initial conditions for
J_1,J_2,θ_1,θ_2
into initial conditions for K_1,K_2. This step requires details of the canonical tranformation. 
Due to their length, the full equations determining the initial values of K_1,K_2 are given in Appendix A (equations 
(<ref>)-(<ref>)). Note there are two cases depending on whether θ_2(0) is 0 or π. Note also that 
there  is no guarantee that these equations will have a solution with real, positive K_1,K_2. 
In the case of the SAT system, for a certain range of parameter values we have experienced numerical problems with
the solution of (<ref>)-(<ref>), specifically for initial values of 
J_2 close to zero, close to the jump from θ_2=0 to θ_2=π. However, typically there are values
of K_1(0),K_2(0) close to the given values of J_1(0),J_2(0). 

Once the initial values of K_1,K_2 have been computed, the values of the constants E and P can
be found and equation (<ref>) can be solved. The right hand side of (<ref>) is the product of two
quadratic factors in K_1, with up to 4 real roots, and typical solutions will be  oscillatory between two roots.
When there is a double root then there is the possibility of the period of the oscillation becoming infinite,
marking a bifurcation in the solution. There are two ways that a double root can occur,  by the vanishing
of the discriminant of one of the quadratic factors, or by one of the roots of the first factor coinciding with one
of the roots of the second. The discriminants of the quadratic factors are

    Δ_1    =   ( (b_2+b_5)^2 - 4b_3 (b_1 + b_4) ) P^2 
        + 2 ( (b_2 - 2 b_3 + b_5) ω_1   +  (-2b_1 + b_2 - 2 b_4 +  b_5) ω_2  ) P  
              + 4  (b_1 - b_2 + b_3 + b_4 - b_5) E + (ω_1 - ω_2) ^ 2  ,  
    Δ_2     =   ( (b_2-b_5)^2 -4b_3 (b_1 - b_4) ) P^2 
        + 2 ( (b_2 - 2 b_3 - b_5) ω_1   +  (-2b_1 + b_2 + 2 b_4 -  b_5) ω_2  ) P 
               + 4  (b_1 - b_2 + b_3 - b_4 + b_5) E + (ω_1 - ω_2) ^ 2  .

A simple algebraic manipulation shows that 
the first factor and second factor have coincident roots if either Δ_3=0 or Δ_4=0,
where

    Δ_3    =     b_3 P^2 + ω_2  P -  E ,   
    Δ_4     =    b_1 b_5^2 - b_2 b_4 b_5 + b_3 b_4^2/(b_4-b_5)^2  P^2 +
                   b_4ω_2-b_5ω_1/b_4-b_5 P -  E .

From (<ref>), we see that the first case occurs when the repeated root is at K_1=0. 

It should be emphasized that the occurence of a double root on the RHS of (<ref>) is a necessary 
condition for a bifurcation of the solution (giving rise to a transition between types) 
but not a sufficient condition.  For example, if the solution
is describing an oscillation on the interval between two adjacent roots of the RHS, and the two other
roots outside this interval merge, this will have no effect on the solution. We illustrate, in Figure 4,
with two concrete examples of equation (<ref>) emerging from the CQ system described in Section 2.
In both cases Qa_0^-2=0.077 and ga_0^4=0;
in the first case r=1.01 and in the second case r=1.045. In both cases we plot the roots of the RHS
as a function of the single remaining parameter E.  (The choice to plot the roots for fixed values of Qa_0^-2, ga_0^4
and r and to vary E is just an illustration; we could just as easilly vary any
of the other parameters or a combination thereof.)
In the first case  there are 4 points P_1,P_2,P_3,P_4
at which there are double roots; however, transitions only occur at the two points P_1,P_4
(marked in Figure 4 with large dots). 
In the second case  there are 5 points P_1,P_2,P_3,P_4,P_5
at which there are double roots; however, transitions only occur at the two points P_2,P_4. 
In both cases, the first transition is from type II to type I, and the second transition is from type I to type II,
as indicated by Roman numerals on the plot. 

The theoretical explanation of this is as follows.
In the first case, r=1.01, there are 4 values of E for which there is a double root. The points
labelled P_3 and P_4 on the diagram
are associated with the vanishing of the discriminant Δ_2; the point labelled P_1 is a double root at 0,
associated with the condition Δ_3=0, and the point labelled P_2 is associated with the vanishing of the
discriminant Δ_1. The motion takes place between the root that is at
K_1≈ 0.00015 and an adjacent root: for values of E below P_2 the adjacent root is below,
for values of E above P_2 the adjacent root is above. Thus the double root at P_1 indicates a value of E
for which there is a bifurcation, and the period of oscillation diverges.
The double root at P_2 is a special solution for which
K_1 and K_2 are constant (looking at (<ref>)-(<ref>) it can be seen that there are 3 kinds of solution
of this type, each corresponding to vanishing of one of the three factors on the RHS of this equation; these are related to the
nonlinear normal modes of the system <cit.>). The point P_2  does
not, however, give rise to a transition in behavior of the CQ system;  the beating period diverges there, but the
type does not change. The double root at P_3 also does not mark a transition. This is precisely the case described above,
in which the oscillation is on the interval between 2 roots, and the other two roots outside this interval merge.
The point P_4, however, does mark a second transition, from type I beating back to type II. 



Proceeding to the second example in Figure 4, there are now 5 cases of a double root. P_2 and P_4 are
associated with the vanishing of the discriminant Δ_2, P_3 with the vanishing of the discriminant Δ_1.
P_1 is the case of a double root at zero associated with the condition Δ_3=0, and P_5 is associated with
the final possibility, Δ_4=0. There are however only 2 transitions, associated with the points P_2 and P_4, for
similar reasons to the case described in the previous paragraph. 

In this section we have 
explained how small oscillations of the original Hamiltonian (<ref>) near its fixed point and 
near (symmetric) 1:1 resonance can be approximated using the integrable Hamiltonian (<ref>)
and the single differential equation (<ref>). 
We have arrived at a simple analytic approximation for beating transitions, viz.  
a necessary condition for  a transition between type I and
  type II beating is the vanishing of one of the four quantities Δ_1,Δ_2,Δ_3,Δ_4 given in
  (<ref>),(<ref>),(<ref>),(<ref>). 
It should be emphasized that this far from trivializes the original problem.
There is substantial complexity hidden in the relationship between parameters 
and initial conditions of the original Hamiltonian and those of 
the integrable Hamiltonian. Also determining which of the vanishing conditions gives a
physical transition can be subtle. In  Section  4
we apply the approximation to the CQ and SAT models from Section 2
and validate its predictions against numerical results. 



§ APPLICATION TO THE MODELS




 §.§ The CQ Model


The potential of the CQ model, given   by (<ref>), has an isolated minimum when 

    a_x = a_y = a_ min≡4 E√(2Q)/3√(E-1) C_0

where C_0>0 is a solution of the equation

    1 - C_0^2 =  1024E^4Q^2g/81(E-1)^3 C_0^6   .

The 1:1 resonance condition is E= E_ res where


    E_ res = 4/1 + √(1 + 65536 Q^2g/81)
  
In the case of zero grade index, g=0, we have C_0=1 and the resonance condition is simply E=2. 
The model is valid if the parameters E,Q,g are chosen so that E≈ E_ res and 
the initial conditions (see (<ref>)) satisfy a_0≈ a_ min and r≈ 1.   

The relevant parameters for the quartic Hamiltonian (<ref>) are 

    ω_1^2    =   81/512 (E-1)^2((2-E)C_0^2+E-1)/E^4Q^2C_0^6
    ω_2^2    =   81/512 (E-1)^3(3-2C_0^2)/E^4Q^2C_0^6
    
      a_1    =   243/8192  (E-1) ^5/2( 2(E-3) C_0^2-3(E-1)  )/Q^5/2C_0^7E^5
    
      a_2    =   243/8192  ( E-1) ^7/2 (2C_0^2-5) /Q^5/2C_0^7E^5
    
      a_3    =   729/262144  (E-1) ^3 ( 2(5-E)C_0^2+3(E-1) )  / C_0^8E^6Q^3
    
      a_4    =   729/131072(E-1)^3( 10(3-E)C_0^2+ 21(E-1) )  /C_0^8E^6Q^3
    
      a_5    =   3645/262144(E-1)^4( 7-2 C_0^2 )/C_0^8E^6Q^3
  

The detailed recipe for checking whether a given set of parameters and initial conditions E,Q,g,a_0,r might give
rise to a transition is as follows:


  * Compute the coefficients ω_1^2,ω_2^2,a_1,a_2,a_3,a_4,a_5 using  (<ref>). This is the only
       stage of the recipe that is model dependent. 
       Compute the coefficients b_1,b_2,b_3,b_4,b_5 from  (<ref>). 

  * Compute the initial conditions J_1(0),J_2(0) from (<ref>) and θ_1(0),θ_2(0) from the comments
  following (<ref>). In the case of CQ, all the parameter values which we used gave θ_1(0)=0 (we
  took r>1 throughout) and θ_2(0)=π. 

  * Compute the initial conditions K_1(0),K_2(0) using (<ref>)-(<ref>).   This is
  the only stage of the recipe that is not completely explicit, and involves solving two equations in two
  variables. If no real solution can be found, the method fails. A suitable initial guess for the solution
  is K_1(0)≈ J_1(0) and K_2(0)≈ J_2(0). 

  * Determine the value of E, the constant value of the Hamiltonian H_3 using (<ref>), taking
  cos(2(ϕ_1-ϕ_2))=1. Determine the value of P=K_1+K_2.   

  * Compute Δ_1,Δ_2,Δ_3,Δ_4 from   (<ref>),(<ref>),(<ref>),(<ref>).
  Values of E,Q,g,a_0,r for which any of these quantities vanish are candidates for transitions. 
  







Figure 5 displays results. Figure 5a shows numeric values and candidate analytic approximations of
r_c as a function of E for Qa_0^-2=0.077 and ga_0^4=0.01. The dots denote numeric values of transitions
in the original system.  The solid curves show candidate analytic approximations of 3 distinct types: 
(1) (black) values for which Δ_2=0 (a closed loop with a cusp on the axis at r=1), 
(2) (green) values for which Δ_3=0 (a simple open curve) and 
(3) (red) values for which Δ_4=0 (two crossing open curves). For the values of
Qa_0^-2 and ga_0^4 specified, it seems there are two branches of parameter values for which
there are transitions. We denote
the lower branch (on the plot) by r_c,1(E), which exists for
E greater than a certain value which we denote by E_c,1, 
and 
the upper branch by r_c,2(E), which exists for
E greater than a certain value which we denote by E_c,2,
with E_c,2≈ 1.975 < E_c,1≈ 1.977. 
On the lower branch, as E increases from E_c,1, r_c,1(E) at first follows the approximation
Δ_2=0, until a triple point at which the curves Δ_2=0 and Δ_4=0 intersect.
As E increases further, r_c,1(E) follows the approximation Δ_4=0. Surprisingly, this
approximation stays reasonably accurate for the full range shown on the figure, even though
r_c,1(E) rises to approximately 1.12.
On the upper branch, as E increases from E_c,2, r_c,2(E) at first follows the approximation
Δ_3=0, until a triple point at which the curves Δ_2=0 and Δ_3=0 intersect.
As E increases further, r_c,2(E) follows the approximation Δ_2=0. However, the quality of this
approximation rapidly decreases as E and r_c,2(E) increase further, with the discrepancy already visible
on the plot for r_c≈ 1.06. 




Figure 5b shows numeric values and the correct analytic approximation (made up of pieces of the
curves Δ_2=0, Δ_3=0 and  Δ_4=0) in the cases (1) ga_0^4 = 0, (2) ga_0^4 = 0.01, (3) ga_0^4 = 0.02,
all for Qa_0^-2=0.077. In addition, stars indicate numerical values of transitions obtained for the
quartic system with Hamiltonian (<ref>). For small values of r, the numerical values for transitions
for the exact Hamiltonian and the approximate quartic Hamiltonian (<ref>) are, as we would expect,
very close. However as r increases, we see that the results for the quartic Hamiltonian rapidly
diverge from the results for the exact Hamiltonian, while, remarkably,  the analytic approximation continues
to be a reasonable approximation for the exact Hamiltonian. This may find an explanation in the fact that while the
exact Hamitlonian (<ref>), the quartic approximation (<ref>) and the integrable approximation (<ref>) all
agree close to the fixed point, the global properties of the exact Hamiltonian are expected to be closer to those of
the integrable approximation than the quartic approximation. 

Note also in Figure 5b that the intercepts of the curves on the E axis, that we have denoted above by
E_c,1 and E_c,2, are very close to the values of E
determined by the resonance condition (<ref>), which are (1) E=2, (2) E ≈ 1.977, (3)
E≈ 1.954.
However, even though the intercepts for the two curves obtained for each set of parameter
values are very close, they are not identical. This is something that is difficult to establish a priori by direct
numerics for the original systems (as the beating periods, for values of r close to 1, are very long), but once 
the analytic approximation is available to give accurate candidate values for the transition locations, it is
possible to verify them a posteriori. Thus 
in the small band of values E_c,2<E<E_c,1 there is only a single beating transition as the beam eccentricity
is increased. 
As r is increased from 1 there is immediately type I beating, and as r is increased further there is only a
single transition to type II (as opposed, for example, to the sitution in Figure 2, where as r is changed from 1 type II
beating is seen, and then there are two transitions). Using the analytic approximation it can be shown (see
Appendix A) that the points E_c,1, E_c,2 are determined by the conditions

    ω_1 - ω_2 + P(b_2-2b_3 ∓ b_5) = 0
  
(minus for E_c,1, plus for E_c,2) for a solution with r=1. (The condition r=1  implies J_1(0) = K_1(0)=0,
and then equation (<ref>) gives a single equation from which to determine K_2(0) from J_2(0) = 
ω_2 (a_0 - a_ min)^2.)  
Figure 6 shows the dependence of E_c,1 and E_c,2
on ga_0^4 for two values of Qa_0^-2, as computed by the analytic approximation, along with a few numeric
values (computed a posteriori).
In addition the value of E_ res from (<ref>) is shown, this being the value of E for which there is 
exact 1:1 resonance in the linear approximation.  
We see that the values of E_c,1, E_c,2 and E_ res all decrease monotonically with ga_0^4. 





 §.§ The SAT Model


The potential of the SAT model, given   by (<ref>), has an isolated minimum when 

    a_x = a_y = a_ min≡   a_0 √(K_ SAT/K_0)
   
where K_0, which depends on the parameters E,K_ SAT, ga_0^4,  is a solution of the equation

    K_0^2/4E +  Li_2 ( -K_0)  + ln(1+K_0) -  K_ SAT^2ga_0^4/4E
    = 0   .
  
(Recall that the constant K_ SAT is defined by K_ SAT=4α^2Ea_0^-2.)
The resonance condition can be written K_0 = K_ res where K_ res is the solution of 

    Li_2 ( -K_ res)  + 2ln(1+K_ res) -  K_ res/1+K_ res = 0  .
  
K_ res has numerical value approximately 5.017. 
We recall that for our analytic model to be most effective we need to be near resonance,
and the initial conditions should 
be close to the minimum, i.e. a_0≈ a_ min, or K_ SAT≈ K_0,  and r≈ 1.
These conditions give  K_ SAT≈ K_ res = 5.017  and E ≈ 6.550( 1 - ga_0^4). In practice
we will look at a large range of values of E and K_ SAT, but focus on this
region. We also recall that in our model the sign of ζ_2(0) (as given in (<ref>) plays a critical role. 
From (<ref>) we have ζ_2(0)=0 (or equivalently K_ SAT=K_0) when  <cit.>

    E =   -K_ SAT^2(1-ga_0^4)/ 4(  Li_2(-K_ SAT)+ln(1+K_ SAT))  .


The relevant parameters for the quartic Hamiltonian (<ref>) in the SAT case are 

    ω_1^2    =   2/a_0^4( K_0^2/K_ SAT^2 + ga_0^4 )  
      
    ω_2^2    =   4/a_0^4K_ SAT^2(  -2Eln(1+K_0) + K_0^2  + 2  E K_0/1+K_0)  
    
      a_1    =   √(2K_0)/K_ SAT^5/2a_0^5( 2Eln(1+K_0) - 3K_0^2-2EK_0/1+K_0)  
    
      a_2    =   √(2K_0)/3K_ SAT^5/2a_0^5(-2Eln(1+K_0)  -3K_0^2  + 2EK_0(1+3K_0)/(1+K_0)^2) 
      
    
      a_3    =   K_0/4K_ SAT^3a_0^6( -2Eln(1+K_0) + 5K_0^2 +  2EK_0/1+K_0)  
         
    
      a_4    =   K_0/2K_ SAT^3a_0^6( -2Eln(1+K_0) + 15K_0^2 + 2EK_0(1-K_0)/(1+K_0)^2)  
         
    
      a_5    =   K_0/12K_ SAT^3a_0^6( 2Eln(1+K_0) + 15K_0^2 -  2EK_0(1+10K_0 + K_0^2)/(1+K_0)^3)   .
  



The method is identical to that given for CQ in the previous subsection, so we can immediately present 
results. For fixed values of E and ga_0^4 we look for values of r giving
beating transitions as a function of K_ SAT. Both numeric and analytic results 
suggest there is a qualitative difference in behavior for E above and below a critical threshold, 
and our results are consitent with the value of this threshold being  approximately 6.550( 1 - ga_0^4),
as found above. Figure 7 displays results for ga_0^4=0 and E=6.3 (below the threshold, left) and
E=6.7 (above the threshold, right). The numeric results show that below the threshold, there are two
ranges of K_ SAT for which there is a single beating transition, from type I (for r below r_c)
to type II (for r above r_c). For values of K_ SAT below or above these two ranges, there is only type I
beating, and for values between the two ranges there is only type II beating.  The analytic approximation
reproduces these results well. In this region of parameter space there are values for which Δ_2=0
(indicated in black in the figure) and
Δ_3=0 (indicated in green). It is the latter that are physically relevant, and the
values of r_c predicted by the analytic model are accurate for a good range. 
Moving “above the threshold”, numerics show there is a range of values of K_ SAT for which, as
r is increased from 1, the beating is initially type II, then there is a transition to type I. For some of these
values there is then a further transition back to type II for quite high values of r. It should be mentioned
that these latter transitions were initially discovered using the analytic approximation, and 
confirmed numerically a posteriori. The analytic approximation reproduces the first transition very well, using pieces of
the Δ_2=0 and Δ_4=0  degeneracy curves. The upper transition is not reproduced well, which is
not surprising bearing in mind the values of r involved. Pieces of the
Δ_2=0 and Δ_3=0  degeneracy curves are close to some of the results, but for a small
range of values of K_ SAT and r the model fails
as there is no solution of equations (<ref>)-(<ref>). Two branches of the Δ_3=0 degeneracy
curve come to an abrupt end
(in the plot we have connected the ends with a dashed line, which is not associated
with any degeneracy).
The values of parameters involved
are precisely those for which ζ_2(0)≈ 0.



Figure 8 enlarges upon these results for different values of E and ga_0^4. In the 4 panels here, the upper panels
(a and b) show results for values of E above the threshold, and the lower panels (c and d) show results for values
of E below the threshold. In the left panels (a and c), ga_0^4=0, in the right panels (b and d), ga_0^4=0.02.
For
ga_0^4=0, the values E=6.3,6.4,6.5,6.7,7.0 are shown, the first three of which are below the threshold (in panel c),
and the last two above the threshhold (in panel a). 
For
ga_0^4=0.02, the values E=6.3,6.4,6.5,6.7 are shown, the first two of which are below the threshold (in panel d),
and the last two above the threshold (in panel b). 
Note specifically that for ga_0^4=0 the case E=6.5 is below the threshold (approximately 6.55), while for
ga_0^4=0.02 it is above (as the threshold drops to approximately 6.42). Thus (for example) for E=6.5, K_ SAT=5 and
ga_0^4=0, no beating transitions are observed as the beam eccentricity is increased; but if the grade index is
changed to ga_0^4=0.02, there are two beating transitions.
The analytic theory fully explains this phenomenon. Indeed, for all the cases shown in Figure 8,
the analytic theory  is in excellent quantitative
agreement 
with numerics for lower values of r, and gives reasonable qualitative predictions 
for higher values of r. 

Another conclusion from Figure 8 is that for values of  E below the threshold, we can find two
values of K_ SAT that give rise to a given value of r_c, but for E above the threshold this need not
be the case; furthermore the gap in r_c values increases with the given value of E. In Figure 9
we illustrate this phenomenon more clearly. For the case ga_0^4=0,
we show contours in the K_ SAT,E plane that give rises to the values r_c = 1/0.95≈ 1.053 (black),
r_c = 1/0.92≈ 1.087 (blue) and   r_c = 1/0.895≈ 1.117 (red). It is clear that the “gap” 
between the two branches of each contour increases with r.  Note that in the upper branch of each contour there is a
small section denoted by a dashed line where the analytic method fails
(the dashed line is a straight line between
the last two points on each side for which the method works). As expected,
the regions where the method  fails straddle the curve (<ref>),
incidicated by a dashed turquoise curve.  Note further that in many cases the analytic method works well far beyond
the region in which this is expected, but there are some exceptions. 





§ CONCLUSIONS AND DISCUSSION
 
In this paper we have described the beating phenomena observed in the equations of motion for the beam
widths obtained in a collective variable approximation to solution of the GNLSEs relevant for beams 
in nonlinear waveguides  with cubic-quintic (CQ) and saturable (SAT) nonlinearities and a graded-index profile.
We have described the different types of beating, and the transitions between them.
Arguing that the origin of these phenomena is in a 1:1 Hamiltonian resonance, we have developed an
approximation scheme for small oscillations in 
a class of 2 degree-of-freedom Hamiltonian systems with an isolated fixed point close to 1:1 resonance.
We have shown that such oscillations can be described by an integrable Hamiltonian, or, alternatively, a single
first order differential equation (<ref>). Understanding the bifurcations of the system, which include the
beating transitions, can be reduced to looking at the bifurcations of the roots of a pair of quadratic equations. 
Applying our general methodology to the specific cases of the CQ and SAT models we managed to reproduce
numerical results for beating transitions over a large range of parameter values. The theory
allows us to map out the regions (of parameter space and beam eccentricities) 
where beating transitions do and do not exist. Amongst other things,
in the CQ case we identified a band of beam energies
for which there is only a single beating transition (as opposed to 0 or 2) as the beam eccentricity is increased,
and  in the SAT case we explained the appearance and disappearance of transitions with
changes of the grade-index. 

We expect our methods to have applications to related problems in nonlinear optics, for
nonlinearities other than the ones studied here,  for different beams, such as super-Gaussian
beams <cit.>, and for optical bullets <cit.>. We are encouraged by the fact
that there is some recent experimental evidence
<cit.> of breathing in optical solitons, albeit in a dissipative setting.
We also hope the general theory of 1:1 resonances that we have developed will find application
in the settings of nonlinear mechanics and astronomy, as well as suitable extensions for 
1:1:1 resonances in higher dimensional systems (see for example the recent papers <cit.>). 







§ FURTHER TECHNICAL DETAILS


As explained in section 3, the Hamiltonian (<ref>) is in an integrable approximation
to the Hamiltonian (<ref>), and is obtained from (<ref>) via a canonical transformation
and neglecting higher order terms. The only need for explicit details of the canonical transformation is
to compute the initial conditions of the variables K_1,K_2 from the initial conditions 
of J_1,J_2 given in (<ref>). The equations to be solved are 

    J_1    =     K_1∓4K_1a_1/4ω_1^2- ω_2^2√(2K_2/ω_2)  
        + (   
    a_1^2 ( 48ω_1^4-8ω_1^3ω_2-40ω_2^2ω_1^2+2ω_1ω_2^3+5ω_2^4)/4ω_2^2ω_1^3( 2ω_1+ ω_2)^2( 2ω_1-ω_2) ^2 
        -5a_3/2ω_1^3)  K_1^2  
            + (  ( 40ω_1^3+28ω_1^2ω_2-6ω_1ω_2^2
            -3ω_2^3) a_1^2/ω_1^2ω_2( 2ω_1-ω_2) ^2( 2ω_1+ω_2)^2( ω_1+ω_2) .  
          . 
    +( 12ω_1^3+11ω_1^2ω_2-10ω_1ω_2^2
                 -6ω_2^3) a_ 1a_2/2ω_1^2ω_2^3( ω_1+ω_2)  ( 4ω_1^2-ω_2^2) 
    - ( 3ω_1+2ω_2) a_4/2ω_1^2ω_2( ω_1+ω_2) ) K_1K_2   ,  
    
     J_2    =     K_2∓ 2 (
       ( 2ω_1^2-ω_2^2)K_1a_1/ω_1ω_2( 4ω_ 1^2-ω_2^2) 
    +  K_2a_2/ω_2^2)  √(2K_2/ω_2)
       + (  33a_2^2/4ω_2^5 -5a_5/2ω_2^3) K_2^2  
           
       + ( 16ω_1^4+8ω_1^3ω_2 - 12ω_1^2ω_2^2
                   -2ω_1ω_2^3+3ω_2^4) a_1^2K_1^2/2ω_1^2ω_2^3( 2ω_1-ω_2)^2(2ω_1+ω_2)^2
          
    + (  
    ( 8ω_1^4+16ω_1^3ω_2-10ω_1^2ω_2^2
                   -8ω_1ω_2^3+ω_2^4) a_1^2/ω_1^2ω_2^2( 2ω_1-ω_2)^2( 2ω_1+2ω_2) ^2( ω_1+ ω_2)        . 
          .       
    +( 40ω_1^3+44ω_1^2ω_2-9ω_1ω_2^2-16ω_2^3) a_1a_2/2ω_1ω_2^4( ω_1+ω_2) 
    ( 2ω_1-ω_2)( 2ω_1+ω_2)  
     
    -( 2ω_1+3ω_2) a_4/2ω_1ω_2^2( ω_1+ω_2) ) K_1K_2  .
  
Here the upper signs should be taken in the square roots terms 
in the case θ_2(0)=0 and the lower signs in the case θ_2(0)=π. 

In Section 4.1, in the study of the CQ system, we stated the conditions (<ref>) for the value r_c
giving a beating transition to tend to 1. We briefly describe the origin of these conditions.
The symmetric
solutions with a_x=a_y of (<ref>), arising from the initial condition r=1, 
correspond to solutions with K_1≡ 0 of (<ref>)-(<ref>). From (<ref>), the values of P and
E for such a solution must evidently satisfy b_3P^2 + ω_2 P - E=0, which is just the condition
Δ_3=0, see (<ref>).  As explained in Section 3, a necessary condition for a beating transition is
the vanishing of one of the quantites Δ_1,Δ_2,Δ_3,Δ_4.  To determine E_c,1 in Section 4.1
we want r_c→ 1 for a solution of Δ_2=0. Clearly this requires Δ_2=Δ_3=0, 
and some simple algebra then gives the condition ω_1-ω_2 + P(b_2-2b_3-b_5)=0.  To determine E_c,2, 
however, is not so straightforward, as for this we want
we want r_c→ 1 for a solution of Δ_3=0, and apparently we do not have two equations. The resolution
of this conundrum is as follows: Although we stated above that the symmetric solutions  of (<ref>) 
correspond to solutions with K_1≡ 0 of (<ref>)-(<ref>), the latter in fact provide a blow up of
the former — there is a 3 parameter family of the latter and only a 2 parameter family of the former.
Solving (<ref>)-(<ref>) in the case K_1≡ 0, we obtain K_2=P (constant),
ϕ_2 = ϕ_2(0) + (ω_2 + 2 b_3P)z, and that ϕ_1 must satisfy the ODE

    ϕ̇_1 = ω_1 + b_2 P + b_5 P cos( 2(  ϕ_2(0) + (ω_2 + 2 b_3P)z  - ϕ_1 ))   .

This latter equation can be solved explicitly, and for a general choice of the constant of integration
will give a complicated function ϕ_1(z). However, for a beating transition 
we seek a solution that is characterized by a single frequency, i.e. we need 

    ϕ_1(z) =  ϕ_1(0) +  (ω_2 + 2 b_3 P)  z

Substituting this in the differential equation, we obtain 

    ω_2 + 2 b_3 P  =  ω_1 + b_2 P + b_5 P cos( 2(  ϕ_2(0)  - ϕ_1(0) ))   .

Since the initial conditions 
ϕ_1(0),ϕ_2(0) take the values 0 or π, we deduce that ω_1 - ω_2 + P(b_2 - 2b_3 + b_5) = 0,
as required. 



§ A TWO TIME EXPANSION APPROACH


In this appendix we outline a two time expansion approach <cit.>
which is an alternative to the procedure  based on canonical transformations
given in Section 3.

We wish to look at solutions of the Hamiltonian system with Hamiltonian (<ref>) and initial conditions (<ref>).
We assume the system has an isolated symmetric minimum at which the system is close to 1:1-resonance. To apply a
two time technique we need to introduce a small parameter ϵ explicitly into the equations. Our systems involve
a number of system parameters, for example in the CQ case, the parameters E,Q,g, for which the resonance condition is
(<ref>). We introduce a small parameter by selecting one system parameter and writing this
as its value at resonance plus a small perturbation. However, for reasons described in <cit.>, the “small perturbation” here
  should be quadratic in the small parameter. Thus, for example in CQ, we have to consider two possibilities,
E = E_ res±ϵ^2
  where E_ res (which depends on the other system parameters Q,g) is the value of E at resonance. The two resulting
  expansions will differ just in signs. This is the counterpart in the two time method of the need to choose θ_2(0)
  to be 0 or π in Section 3 and the resulting choice of signs in equations (<ref>)-(<ref>). 
  However, we emphasize that it is not the same, so the resulting method is different, in particular, the “choice” in Section 3
  involves the initial conditions as well as the system parameters.

  Taking, as before,  the minimum of the potential V to be at a_x=a_y=a_ min we now write
  
    a_x = a_ min + ϵã_x  ,      a_y = a_ min + ϵã_y

  and expand to 4th order in ϵ. The order 0 terms are irrelevant and can be discarded. 
  The order 1 terms vanish by definition  of a_ min.
  In the other terms there is dependence on all the system parameters. However by making the
  assignment of the form E = E_ res±ϵ^2, discarding all terms of order higher than 4
  and a suitable rescaling, we obtain an approximate potential of the form 

    Ṽ   =   1/2 C_1  (ã_x^2 +   ã_y^2)       
      + ϵ( C_2( ã_x^3 +   ã_y^3)  + C_3ã_xã_y  (ã_x +   ã_y)   )  
            + ϵ^2 ( C_4( ã_x^4 +   ã_y^4)  + C_5ã_xã_y (ã_x^2 +   ã_y^2)
         + C_6 ã_x^2ã_y^2   +  C_7(ã_x^2 +   ã_y^2)       + C_8 ã_xã_y   )  .

Here C_1,…,C_8 are all functions of the system parameters excluding the parameter replaced by ϵ. 
Note that as a result of the dependence of the system parameters on ϵ there are now quadratic terms
in ã_x,ã_y in the O(ϵ^2) terms. 

Following the usual two time formalism, we seek solutions of the system with potential Ṽ in the form 

    ã_x    =   Λ_1(ϵ^2 z)  cos (√(C_1)z) + Λ_2(ϵ^2 z) sin (√(C_1)z) +
           ϵã_x,1 (z,ϵ^2 z) + ϵ^2 ã_x,2 (z,ϵ^2 z) + …
    ã_y    =   Λ_3(ϵ^2 z) cos (√(C_1)z) + Λ_4(ϵ^2 z) sin (√(C_1)z) +
      ϵã_y,1 (z,ϵ^2 z) + ϵ^2 ã_y,2 (z,ϵ^2 z) + …

Here Λ_1(ϵ^2 z), Λ_2(ϵ^2 z), Λ_3(ϵ^2 z), Λ_4(ϵ^2 z) are
functions of the slow variable ϵ^2 z. Substituting in the equations of motion and equating order-by-order,
the first order terms ã_x,1, ã_y,1 can be determined, and a system of first order equations
is obtained that Λ_1,Λ_2,Λ_3,Λ_4 must satisfy to guarantee the absence of secular terms
in ã_x,2, ã_y,2. Writing

    R_1    =   Λ_1^2 + Λ_2^2 + Λ_3^2 + Λ_4^2  
    
      R_2    =   Λ_1^2 + Λ_2^2 - Λ_3^2 - Λ_4^2  
    
      R_3    =   Λ_1 Λ_3  + Λ_2 Λ_4  
    
      R_4    =   Λ_1 Λ_4  - Λ_2 Λ_3

(c.f. <cit.>) we obtain the system

    R_1'    =    0  
       
      R_2'    =    4 R_4( γ_1 R_1 + γ_2 + (γ_3+γ_4)R_3) 
    
      R_3'    =    -γ_3 R_2 R_4   
    
      R_4'    =      - R_2( γ_1 R_1 + γ_2 + γ_4 R_3)

where the constants γ_1,γ_2,γ_3,γ_4 are certain combinations of the constants C_1,…,C_8. 
(Note R_3^2 + R_4^2 = 1/4 ( R_1^2 - R_2^2 ).)  Thus R_1 is an invariant, as are the quantities 

    Q_2 = R_2^2 + 4(1 + γ_4/γ_3)
                      (  R_3 + γ_1 R_1 + γ_2/γ_3+γ_4)^2 
     ,
           Q_3  = R_4^2  -  γ_4/γ_3( R_3 + γ_1 R_1+ γ_2/γ_4)^2  .
 
Note that R_1,Q_2,Q_3 are related by

    Q_2 + 4Q_3 = R_1^2 - 4(γ_1R_1+γ_2)^2/γ_4(γ_3+γ_4) .
 
Using the invariants it is possible to write a single differential equation for the quantity R_3:

    (R_3')^2 =    -4 γ_4 (γ_3 + γ_4  )
    (   (  R_3 + γ_1 R_1 + γ_2/γ_3+γ_4)^2
         -γ_3/4(γ_3 + γ_4 )    Q_2    ) 
    (  ( R_3 + γ_1 R_1+ γ_2/γ_4)^2
         + γ_3  Q_3/γ_4) .

This has the same form as (<ref>) —  the right hand side is a product of two quadratic factors in R_3 —
and similar techniques can be used to discuss bifurcations of its solutions. Specifically, there can be a double root
if the discriminant of one of the factors vanishes (i.e. if Q_2 or Q_3 vanish), or if the factors have a common root.
The latter happens in the two cases 

    ( (2γ_1±γ_4)R_1+2 γ_2)^2+ 4Q_3 γ_3γ_4  = 0 .
  
As in Section 3, detecting beating transitions requires translating the initial conditions to the constants of
motion R_1,Q_2,Q_3 and checking up to 4 conditions.

We have implemented this method for the CQ and SAT systems and found some satisfactory results which
we do not report here; in certain cases the results were better than those found using the method based on
canonical transformations. However there are numerous reasons to prefer the method based on canonical transformations. 
The two time method requires deciding how to explicitly introduce a small parameter and different ways
of doing this give different results. It also requires advance knowledge of the correct relative order of
magnitude of the oscillations around the fixed point and the deviation of the system parameters from their
resonance values. In general,
the algebraic manipulations required to implement the two time method, most of which
we have omitted in our account here, are substantially more complicated than those required for the method based
on canonical transformations; in particular the reduction of the system (<ref>) 
to a single differential equation (<ref>) is a surprise, that emerges from ad hoc manipulations,
whereas the parallel steps in the canonical formalism are standard, based on the integrability of the Hamiltonian
(<ref>). Finally, from our numerical experiments it emerges that
while the results based on the vanishing of the discriminant
of one of the factors of the right hand side of (<ref>) are good, the results based on conditions (<ref>)
are poor. 












acm


