Stability properties of periodically driven overdamped pendula and their implications to physics of semiconductor superlattices and Josephson junctions

# Stability properties of periodically driven overdamped pendula and their implications to physics of semiconductor superlattices and Josephson junctions

Jukka Isohätälä Department of Physical Sciences, P.O. Box 3000, University of Oulu FI-90014, Finland    Kirill N. Alekseev Department of Physics, Loughborough University LE11 3TU, United Kingdom Department of Physical Sciences, P.O. Box 3000, University of Oulu FI-90014, Finland
July 7, 2019
###### Abstract

We consider the first order differential equation with a sinusoidal nonlinearity and periodic time dependence, that is, the periodically driven overdamped pendulum. The problem is studied in the case that the explicit time-dependence has symmetries common to pure ac-driven systems. The only bifurcation that exists in the system is a degenerate pitchfork bifurcation, which describes an exchange of stability between two symmetric nonlinear modes. Using a type of Prüfer transform to a pair of linear differential equations, we derive an approximate condition of the bifurcation. This approximation is in very good agreement with our numerical data. In particular, it works well in the limit of large drive amplitudes and low external frequencies. We demonstrate the usefulness of the theory applying it to the models of pure ac-driven semiconductor superlattices and Josephson junctions. We show how the knowledge of bifurcations in the overdamped pendulum model can be utilized to describe effects of rectification and amplification of electric fields in these microstructures.

Pendulum and pendulum-like equations are arguably among the most important classes of equations in modern nonlinear sciencesagdeev (). The most often encountered representatives of this family may well be the driven and damped pendulum, , and its first order counterpart, the overdamped pendulum, , the latter type being the topic of this paper. These equations appear, for instance, in the well-known Stewart-McCumbermccumber (); stewrt () and Aslamazov-Larkinaslamazov69 () models of Josephson junctions. The sinusoidal nonlinearity gives rise to a wide class of nonlinear phenomena that have important practical applications: the ac-Josephson effectlangenberg (); levinsen () and the modern voltage standardhamilton () are prime examples of this. More recently, pendulum equations have been found in the theory of semiconductor superlattices where they frequently occur in the limiting cases of the governing differential equationsalekseev98:spont-dc (); cannon01 (); alekseev02a (); dodin98b (); alekseev05:lsslcr (); zharov06 (); dodin03 (). Moreover, overdamped pendulum equations are often encountered in mathematical models of synchronization of nonlinear oscillatorspikovsky01 (). Further recent interest in the overdamped pendula has come from the field of high- superconductors: It has been demonstrated that stacked array of intrinsic Josephson junctions in magnetic field can be synchronized and described by overdamped pendulum-like dynamicsgaifullin08 (); gaifullin09 (). Properties of ac-driven overdamped pendula are also of importance in theories of the amplification of microwave radiation in Josephson point contactskuzmin79 (); kuzmin80 (); *kuzmin80-orig; velichko99 (); likharev-c10-11 (). And last but not least, in our previous workisohatala05:sbpend () we demonstrated how instabilities occurring in the overdamped pendulum are carried over to higher dimensional systems such as the strongly damped second order pendulum equation. Our present paper has two sides: mathematical and physical. Here we develop a mathematical technique which allows to find bifurcations in a class of overdamped pendula models for a wide range of their parameters, including a difficult but physically interesting case of low frequencies of driving force. We also show how this technique can be applied to symmetric physical systems demonstrating pendulum dynamics in some limiting cases. Our main focus is on the rectification and amplification of microwave radiation in unbiased semiconductor superlattices and Josephson junctions.

## I Introduction

We consider the first order ordinary differential equation with a sinusoidal nonlinearity and arbitrary time dependence

 ˙θ(t)+G(t)sinθ(t)=F(t). (1)

The dynamics of the overdamped pendulum has been studied by several people, motivated by direct physical applications mentioned above. In spite of its apparent simplicity, novel nonlinear dynamics, e.g. strange nonchaotic attractorsromeiras87 () have been found. Here, we restrict ourselves to a specific class of periodic forcing, specifically consider bifurcations occuring in this systems, and its applications in physical systems.

We will implicitly assume everywhere that and are real, continuous, and differentiable sufficiently many times. Our focus will be on functions and that have the following property

 F(t+T/2)=−F(t),G(t+T/2)=G(t). (2)

With the above choice of external time-dependence, Eq. (1) remains invariant under the transformation

 t→t+T/2,θ(t)→−θ(t+T/2)+2kπ, (3)

and is an integer. This type of forcing and the associated symmetry are of interest in many pure ac-driven physical systems, in particular bulk semiconductors and semiconductor superlattices, where breaking of symmetry (3) implies generation of a spontaneous dc biasbumyalene89 (); alekseev98:spont-dc (). In our previous workisohatala05:sbpend () we in passing considered Eq. (1) with and . We observed that the only instability that occurs is an exchange of stability between two periodic solutions following symmetry (3) and having the properties and , where stands for time-average across the period of the solution.

We found that the symmetry breaking bifurcation, () dhumieres82 (); swift84 () of the strongly damped second order pendulum

 ¨θ+γ˙θ+sinθ=fcosωt, (4)

reduced to this instability in the limit of very large damping. We conjectured that other ac-driven systems reducible to Eq. (1) undergo a type of bifurcation similar to the one found in the strongly damped pendulum equation near the points where the exchange of stability occurs. This motivates our present extended study of the stability properties of system (1) and applications to a number of physical systems.

In this paper, our mathematical analysis is based on mapping of equation (1) to a particular second order linear differential equation. Such transformations have proved useful in the study of various linear and nonlinear differential equations. Following Prüfer’s application of the idea to Sturm-Liouville problemsprufer26 (), these changes of variables are sometimes called Prüfer transforms. In a sense the reverse of this approach was taken in by Bondeson et al.bondeson85:quasi-pend () where the authors used a similar transformation to study quasiperiodically driven overdamped equation by relating the problem to a Schrödinger equation with a quasiperiodic potential. We take essentially the same approach, but focus on the more specific problem of periodically driven equation.

On the other hand, applications considered in this paper are based on the connection of an exchange of stability in the pendulum with the physical phenomena of amplification and rectification. Here we consider effects of microwave rectification and amplification in two pure ac-driven systems reducible to the overdamped pendulum: single-band lateral semiconductor superlattice and point-contact Josephson junction.

Formally, by rectification we mean the conversion of pure ac excitation into response at even harmonics of some quantity that is a odd function of ; for instance itself or . As an example of rectification, consider Eq. (4) as a toy model where the drive corresponds to some ac applied field and that current is . Rectification would then imply , i.e. obtaining a direct current response from a pure ac excitation, hence the term “rectification”. For symmetric solutions rectification is impossible, since will only have odd harmonics, excluding possibly zeroth harmonic that is a multiple of . Thus, symmetry breaking is a prerequisite for rectification.

Note that rectification due to spontaneous breaking of symmetry in solutions, Eq. (3), should be distinguished from phase-dependent rectification due to breaking of symmetry in the equationsflach00 (). The latter requires explicitly introducing that does not follow Eq. (2), for example an additional phase-shifted second harmonic .

For large damping we get overdamped first order pendulum for which exchange of stability arises for same parameters as symmetry breaking in second order pendulum. Rectification in pendulum is, however, rather artificial model which does not correspond to any real physical system. Nevertheless, in dynamical systems describing realistic physical situations, symmetry breaking bifurcation is realized near values of parameters that are close to the values necessary for the exchange of stability in the overdamped pendulum.

Here we consider a model of ac-driven lateral semiconductor superlatticedodin98b (); alekseev05:lsslcr () which is described by two first-order nonlinear balance equations which can be reduced to a sort of overdamped pendulum (see Eq. (30) in the limit of strong nonlinearity. Symmetry breaking in balance equations of lateral superlattice corresponds to rectification of applied ac electric field alekseev05:lsslcr (). We demonstrate how analysis of instabilities in Eq. (1) can provide a quite useful information on the parameter space of rectification in these nanostructures.

Our another application is related to amplification of infinitesimally weak signal in Josephson point contact described by Eq. (1), in which – voltage, – current and itself is difference of phases of wave functions of superconductors in the junction. As a rule, an amplification of a small signal is observed near the onset of a dynamical instabilitywiesenfeld85 (); wiesenfeld86 (). Here we show a small-signal amplification near an exchange of stability and at frequencies of signal close to even harmonics of pump. Therefore, this effect of amplification can be considered also reminiscent of symmetry breaking bifurcation in strongly damped second order pendulum Eq. (4). Despite here even harmonics are forbidden by symmetry, nontrivial amplification of additional weak signal does exist at pump amplitudes and frequencies close to those necessary to realize real symmetry breaking.

The outline of this paper is as follows. In the next section we will introduce the change of variables that yields a second order linear differential equation and briefly recapitulate on some known properties of its solution and their implications on Eq. (1). We will then proceed to the more specific problem of forcing following Eq. (2) and show that an exchange of stability is the only instability occuring in this system. We then consider perturbations of the Eq. (1) and essentially prove our earlier conjecture that the exchange of stability is a limit of pitchfork bifurcations. In the subsequent section, we shift to a more practical approach: An approximate condition for the instability to occur will be derived in nontrivial case of large and . Finally, we go on to apply the results to relevant physical problems. Technical details are presented in three Appendixes.

## Ii Equivalent linear equation

We start by introducing a change of variables from to new variables as

 θ(t)=2arctan(q1(t)q2(t)). (5)

We will denote the vector by and the change of variables by . Since Eq. (5) alone does not fix the functions , , we have some freedom in choosing the differential equations for the new variables. Here we opt for a particularly symmetric form of the equations

 ddt(q1(t)q2(t))=12(−G(t)F(t)−F(t)G(t))(q1(t)q2(t)). (6)

The coefficient matrix on the right-hand side of Eq. (6) will be denoted by . We consider only periodic and , and therefore Floquet theory can be directly applied to the problem. We adopt the following notations for the Floquet solutions , :

 Φi(t)=eBitPi(t)=eRe{Bi}t~Pi(t), (7)

where are the (complex) characteristic exponents and are -periodic functions. Functions contain the oscillating parts of the Floquet solutions. Additional lower indices will label the component of , , and , e.g. . Due to the vanishing trace of , the characteristic exponents have always the form (a) , where is real and are integers, or (b) where is real and not an integer. In the next section, we will show that the latter case is never realized if symmetry (2) applies, and thus case (b) will not be addressed in what follows. However, the special case of () will be covered.

The Floquet solutions give a complete description of dynamics of . Supposing case (a) from above holds, the general solution of the pendulum equation is

 θ(t)=C[cosψ02eB0t~P1+sinψ02e−B0t~P2], (8)

where , is a constant that depends on the initial value of , and whose value over uniquely determines the solution up to modulo . Clearly, when there are exactly two periodic solutions , , that are simply given by the Floquet solutions:

 θi(t)=2arctan(Φi,1(t)Φi,2(t)). (9)

From Eq. (8) it follows that stable solutions of the overdamped pendulum correspond to the unstable solutions of the linear equation and vice-versa. This can be also seen from the relation that appears in Refs. johnson82, ; bondeson85:quasi-pend, and also applies here

 −Λ=1T∫T0G(t′)cosθ(t′)dt′=2|B0|, (10)

where is the average exponential rate of growth for a infinitesimal perturbation of . Negative of its absolute value coincides with the maximal Lyapunov exponent of Eq. (1). Periodicity of the asymptotic solutions in the sense that also follows immediately from the form of the Floquet solutions.

## Iii Symmetric case

We now turn to our findings regarding Eq. (6) with forcing following Eq. (2). Matrix transforms in the -shift as , where

 ~e=(100−1). (11)

This property enables us to write the principal matrix , , , on the latter half of a drive cycle in terms of the former

 U(t+T/2)=~eU(t)~eU(T/2). (12)

The monodromy matrix can now be factored into a square of the matrix , and the eigenvalue equation determining the Floquet solutions, , can be solved using the matrix instead of . Determinant of equals , and thus the solution to the characteristic equation becomes

 B1,2=⎧⎪⎨⎪⎩−2Tarcsinh(12tr~M),−2Tarcsinh(12tr~M)+i2πT. (13)

This shows are always real or real plus an integer multiple of . We will make frequent use of real part of , .

For the periodic parts of the Floquet solutions, the following now holds. Applying Eq. (12) and one finds that . Component-wise this property reads:

 Pi,1(t+T/2) = Pi,1(t), (14a) Pi,2(t+T/2) = −Pi,2(t). (14b)

That is, the first component of is -periodic, and the second -antiperiodic. It immediately follows that the periodic solutions of Eq. (1) are symmetric in the sense of Eq. (3). Note that our definition for , Eq. (13), includes an imaginary component, which contributes to the oscillating part of . Therefore, is not real, and it is then more convenient to use the functions instead. Now , and so has the same periodicity as . On the other hand , and thus is -antiperiodic and is -periodic.

Eq. (14) also determines two properties regarding the rotations and the average value of the periodic solutions . Here we assume that roots of are simple, that is, if then . We aim to connect the number of zeros of over , here denoted , to physically relevant properties of – it will be shown that indeed has significance to dynamics of real physical systems we are considering. Since is symmetric, we can write , where counts the positive direction crossings of the line . From Eq. (9) it can be seen that these crossings occur at simple zeros of , and from Eq. (1) that the direction of the crossing is given by the sign of . Using again the symmetry, the average of over , , is . Now, it is easy to see that the parities of and are the same, and thus . From Eq. (14) it is clear that is odd for and even for , and so and , both modulo . This shows that our previous findingisohatala05:sbpend () regarding the averages of regarding the case and holds in general.

Further, relates to the oscillations of and other quantities that are periodic . For instance, consider is such that () for (). As rotates it passes the upright vertical position times and always in the positive direction, and so gives the minimum and maximum number of oscillations of or in one half drive period. This has implications to the physical systems we are considering, since in these and have relevant physical interpretations.

### iii.1 Exchange of stability

Having established that the characteristic exponents are always real plus integer multiples of , it then follows that the only possible type of instability is an exchange of stability where one Floquet solution loses stability and the other gains it. This in turn occurs when and as consequence . Noting that eigenvalues of are never equal, one finds that vectors , , are linearly independent for all and for any . Thus, the Floquet solutions never map to same solution of Eq. (1) and the corresponding asymptotic solutions never cross each other as a parameter is varied. Consequently, the stability is exchanged without the solutions colliding, in contrast to a transcritical bifurcation. We will later show, however, that the exchange of stability can be seen as a type of pitchfork bifurcation. From Eq. (8) it is clear that when all solutions to the overdamped pendulum are periodic. Since the superposition that gives the general solution, Eq. (8), does not have periodicity analogous to Eq. (14), the corresponding solutions are not symmetric.

Further, crossing the instability has a clear effect on some relevant quantities. As was discussed above, the number of simple roots of relates to the rotations and the average of . Let then be the unstable Floquet solution, the stable periodic solution of the overdamped pendulum, and the number of simple roots of over . Because at the instability switches between being and , changes by one. Consequently, the average value of the stable periodic solution jumps by , and further, since the value of is intimately connected to the oscillations of and , these quantities exhibit a change in in their frequency spectrum. We will later apply this finding in the section on lateral semiconductor superlattices. Clearly, the integer partitions the parameter space into disjoint regions with the instability separating them. Thus, serves as a convenient label for different regions of parameter space.

To help illustrate the bifurcation, we introduce a Poincaré map . Naturally we take this to be the stroboscopic map, defined so that for some fixed which we take to be zero. If , the two fixed points , , of are given in terms of the eigenvectors of , or Floquet solutions at , :

 θ∗i=2arctan(Φi,1(0)Φi,2(0)). (15)

From the discussion above it follows that at , for all . In Fig. 1 we have plotted a representative bifurcation diagram. We take and , where , and plot the fixed points as a functions of the forcing amplitude . From the diagram the bifurcation scenario can be easily visualized. At (outside the plot range) we have two fixed points: , , where the latter is naturally the unstable point, since it corresponds to the upright position of the pendulum. At we find the first bifurcation. The initially unstable state becomes the stable one and vice-versa as the critical is crossed. Exactly at the bifurcation, every point is a fixed point of the Poincaré map. From there on, increasing further, we find the the fixed points exchange their stabilities again at and , with the marginally stable points spanning the whole phase space exactly at the bifurcation.

The bifurcation described above bears resemblance to a pitchfork bifurcation. In fact it can be seen as a degenerate pitchfork bifurcation (PB), since two new branches of fixed points emerge at the critical point. By degeneracy, we mean that these branches exist only exactly at the bifurcation point, span the whole phase space, and and are thus only marginally stable. This is in contrast to the (non-degenerate) PB, where the two additional branches of fixed points exist before or after the bifurcation, that are either stable (supercritical case) or unstable (subcritical case). We note that this is reminiscent of the scenario observed for the second order, strongly damped pendulumisohatala05:sbpend (), where a (non-degenerate) pitchfork bifurcation was found near the criterion for exchange of stability in the overdamped equation.

In terms of normal forms, non-degeneracy is understood as non-vanishing of a number of higher order derivatives with respect to the variable of the flow at an equilibriumkuznetsov98 (). We will see that this is indeed the case when we map Eq. (1) to an autonomous equation that is effectively a normal form on the circle . The form of the general solution, Eq. (8) suggests a natural way of mapping Eq. (1) into an equivalent autonomous form: the periodic parts describe the rotations of , while the exchange of stability is determined by the exponential factors and the superposition phase . It seems natural to use a trial function where the constant phase is replaced by a time-dependent , that also accounts for the exponential factors:

 θ(t)=2arctan⎛⎜⎝cosψ(t)2~P1,1(t)+sinψ(t)2~P2,1(t)cosψ(t)2~P1,2(t)+sinψ(t)2~P2,2(t)⎞⎟⎠. (16)

We substitute this into Eq. (1) and use Eq. (6) to obtain an equation for :

 ˙ψ=−2B0sinψ. (17)

No approximations were needed to derive Eq. (17). This equation describes the approach to limit-cycles for the whole class of systems. Obviously, it is also the simplest non-trivial overdamped pendulum that has symmetry (2). Eq. (17) has equilibria , when , in the case every point is an equilibrium, cf. Fig. 1. To compare bifurcations of Eq. (17) to the pitchfork bifurcation, we recall that normal form of a PB is , where is the bifurcation parameterkuznetsov98 (). From Eq. (17), we see that near the equilibrium (), , where . Thus, we see that whereas in the non-degenerate case, only the leading order term vanishes at the bifurcation, here the right-hand side becomes identically zero at .

To help visualize how the Floquet solutions relate to the solutions periodic solutions of the pendulum, we have plotted in Fig. 2 the periodic solutions , together with the Floquet solutions . We have used the same and as above and take to be near the third bifurcation shown in Fig. 1. The left-hand side subfigures (a), (c), and (e) show the solutions for which is just smaller than the critical value of , while subfigures (b), (d), and (f) are plotted for . From (a) and (b) it can be seen that the change little across the bifurcation point, only the stability is exchanged. Similarly, the Floquet solutions remain roughly unchanged as the as the critical is crossed, excluding the fact that the exponential envelope switches from decaying to diverging or vice-versa.

### iii.2 Effect of perturbations

From Eq. (17) it is evident that near a bifurcation the system is structurally unstable. An important question is then how dynamics change when a perturbation that allows for breaking of the symmetry (3), or explicitly breaks (2), is introduced. An exhaustive study of such perturbations is beyond the scope of this paper. Nonetheless we wish to show that the new type of dynamics appear at the exchange of stability when a perturbation is introduced. This is because to a large extent our motivation has been to show that the exchange of stability is in a sense a limit of bifurcations that occur in realistic physical systems that reduce to the overdamped pendulum. The perturbation is then to be understood as the terms removed from the original nonlinear system describing the real physical system to obtain the overdamped pendulum. Therefore we wish to show that symmetry breaking, or other type of dynamics appear exactly at the exchange of stability.

We use trial function of the form given in Eq. (16) to probe the response of the system to small additional terms. We introduce a perturbed system

 ˙θ+G(t)sinθ=F(t)+εH(θ,t), (18)

where and does not necessarily hold. Substituting the trial function of Eq. (16) into Eq. (18) we obtain the equation

 ˙ψ = −2B0sinψ−εΞ(ψ,t)TΞ(ψ,t)H[θ(ψ,t),t], (19)

where . We have fixed the normalization of the Floquet solutions so that isohatala09_footnote2 (). Unlike in the case of Eq. (16) further approximations are needed. Smallness of and near a bifurcation can be used to simplify the above equation.

First, we prove our earlier conjecture that the exchange of stability is the limit of pitchfork bifurcations appearing in more realistic systems following Eq. (2) or equivalent. We consider that follows the symmetry , but contains even harmonics of . We note that the trial superposition has the property . In this case, it then follows that also Eq. (19) has symmetry (2). Using the averaging methodverhulst05 (), we find that the averaged equation for will have the form

 ˙ψ≃−2B0sinψ+εa1sinψ+εa2sin2ψ+⋯ (20)

where are constants which in general are nonzero. Cosine harmonics of in the averaged equation are forbidden by symmetry. Possible bifurcations of the perturbed system can then be qualitatively sketched by fixing and plotting the roots of . As an example, if and for , we find that the degeneracy of the pitchfork bifurcation has been lifted. A representative bifurcation diagram is shown in Fig. 3, where the equilibria of Eq. (20) are plotted for , for . For comparison, the inset shows the degenerate limit of . It can be seen that the degenerate PB is replaced by a pair of non-degenerate PBs that are supercritical for and subcritical for . Between the bifurcations, an equilibrium exists that is not equal to or , and thus is a symmetry broken solution of Eq. (20). In this case, the solution of the pendulum equation will be described by a superposition of two Floquet solutions, and therefore, it too will not in general have symmetry (3). Thus we see that the exchange of stability is the limit of a symmetry breaking PB for symmetrically perturbed systems.

If the perturbation has only explicit, -antiperiodic time dependence, and , then the system will respond strongly when the perturbation frequency is close to an even multiple of . This can be seen by noting that has -periodic zeroth and first cosine harmonic of . Therefore, the right-hand side of Eq. (19) will contain cosine harmonics of whose coefficients oscillate at a frequency . Using again the averaging method, these terms will not vanish but contribute to slow, large amplitude oscillations of . We will later show that this effect has interesting consequences in the problem of weak signal amplification in Josephson junctions. Note also that this is in effect a dual of the pitchfork bifurcation described above – difference is that here the response follows from a near even harmonic of the drive, not the angle .

Finally, if the perturbation does not follow the symmetry (2) but still depends on , one expects to see the cosine terms appearing in Eq. (20). Depending on the perturbation there are several possible outcomes. The various bifurcation scenarios can then be enumerated by selecting the coefficients on the right-hand side of Eq. (20). As an example, the pitchfork bifurcation may become an imperfect PBkuznetsov98 (), or essentially a saddle-node bifurcation, or the equilibria may be destroyed altogether near exchange of stability.

In summary, we described the scenario for the development of only instability occuring in the system of type Eq. (1) following symmetry Eq. (2). The scenario follows the one found in our previous workisohatala05:sbpend (). At a critical point a symmetric solution of Eq. (1) loses stability and another one gains it. The exchange of stability also marks the point where the average shifts from a minimum of the potential to a maximum or vice-versa. Exactly at the critical point neither of the Floquet solutions is diverging (and both indeed are still linearly independent), and hence the system will remain in some initial superposition indefinitely. In this special case, the symmetry (3) need not be satisfied. This state is, however, only marginally stable, since a perturbation will neither decay nor diverge, and occurs only in a null set of parameter values.

We conclude by reviewing the some of the analytic methods of approaching the overdamped pendulum equation to the exact linearization used here. The point of commonality to these methods is that they in one way or another consider the nonlinearity as a perturbation. For instance, the often used technique of using single harmonic trial function has requires the assumption that the sine term only has the effect of changing the amplitude and phase of the otherwise sinusoidal solution. Although the number of harmonics included in the truncation can be increased, the equations quickly get intractable. Averaging method can also be employed, as it indeed was in our analysis of Eq. (20), however, again the potential term needs to be small compared to the time scale at which varies. In contrast, the mapping to the linear equation, Eq. (6) fully retains the nonlinearity whilst still allowing the use of tools for periodically forced linear equations, such as Floquet theory.

Treating the nonlinearity nonperturbatively also allows us to also approach the limit of low frequency driving. In the following section, we derive approximate analytical formulas for the solutions of the overdamped pendulum. Instead of studying the actual time dependence of , we continue with the focus on finding the critical points where the exchange of stability occurs.

## Iv Asymptotic solution of Eq. (6)

We consider next approximate solutions of Eq. (1) in the non-trivial limit of large and , or equivalently, low-frequency. We introduce a large parameter into the problem by making the change , . In the leading order of , we find that the equation we need to solve comes out as

 ¨y+λ24(F(t)2−G(t)2)y=0, (21)

where or . Eq. (21) is found by taking the derivative of Eq. (6). Keeping only terms of order , one finds the equation (21) for both and separately. We note that Eq. (21) does not share the symmetry of Eq. (6), however, results of the previous section allow us to construct an approximate solution that has the expected properties. This follows from the fact that we need only to solve for the first half of a drive cycle and if need be, use Eq. (12) to obtain the complete solution.

Standard methods of asymptotic analysismurray84:asympanal () can be applied to Eq. (21) to find its piecewise solution in the form (see Appendix A).

 yi(t)=1|R(t)|1/4(aiewiξi(t)+bie−wiξi(t)), (22)

where , and () for (). We denote by the turning-points (points such that ) and by their number, i.e. . Additionally, and . Coefficients are solved from the initial conditions, and standard connection formulas for adjoining subintervals are applied. Using Eq. (22) it is straightforward to construct a solution to any particular .

Although piecewise solutions naturally can be cumbersome, we next show that tractable formulae can be obtained for quantities of interest. Naturally, we apply Eq. (22) to calculating , as its zeros define the critical curves of the system. For simplicity, we limit the discussion to the case of two turning points, again, generalizations are straightforward. With this restriction, the trace of becomes

 tr~M ≃ −2sgnG(0)sinh((κ0+κ1)λ+ln2)cosω1λ (23) − 2sgnG(0)cosh((κ0−κ1)λ)sinω1λ,

where , , and . This equation is one of the central results of this paper, as it allows for calculating the critical curves of overdamped pendulum equations in the non-trivial limit of large and .

In order to keep the treatment more concrete, we now fix and . We consider this restriction reasonable, since it was demonstrated in the previous section that differential equations of the form (1) exhibit the same structure as long as have the appropriate symmetry. Thus, we can choose any such forcing as a representative of the class of equations we are considering. Further, this choice has particular relevance to the physical systems we have in mind.

Interestingly, with this choice, Eq. (21) becomes the Mathieu equation. In addition to constructions (22) and 23), we have the known solutions at our disposal. The trace of comes out as

 tr~M≃−1ωS(f2−28ω2,f216ω2;π). (24)

Here denotes the odd solution to the canonical form of the Mathieu equationabramowitz (), , with the (non-standard) initial condition . The quantity vanishes at parameters for which is periodic, whether that period was or . Values of corresponding to a periodic are the Mathieu characteristic values , , and thus the critical curves are described by the equation

 f2−28ω2=bk(f216ω2),k=1,2,… (25)

In addition to Eq. (25), which can be used to calculate the parameters for which vanishes, for practical applications we also need a way of computing for any given . Eq. (24) is not well-suited for this purpose since it requires the use of Mathieu functions with arbitrary parameters. Eq. (23) on the other hand has a more tractable form as it only requires the use of elliptic integrals:

 tr~M = −2sinh(Re{E(f2)}ω+ln2)cosIm{E(f2)}ω (26) − 2sinIm{E(f2)}ω,

where is the complete elliptic integral of the second kindabramowitz ().

In Fig. 4 we have plotted the critical curves as given by Eq. (25) together with the correct numerically obtained ones. For comparison, we have included a high-frequency approximation to the critical lines, where is the Bessel function of order zeroisohatala05:sbpend (). The lines where equals a root of are only in modest agreement for low frequencies, but improves as is increased. Our new result, Eq. (25) is, on the other hand, in very good agreement for low frequencies, and is accurate also in the opposite case of , especially for large . Although not plotted, the condition with given by Eq. (26) also provides very good agreement to the computed critical curves.

A quantity that will often be needed is the average of , where is the stable solution. Using Eqs. (10, 13) we find the following equation that allows us to express simply in terms of as

 ⟨Gcosθ⟩=2|B0|=2T|arcsinh(tr~M/2)|. (27)

In summary, the main result of this section is Eq. (23) whose roots give the critical curves for that are either large or depend slowly on time, obey symmetry (3), and have exactly two distinct points such that . For the case of constant and sinusoidal our three main results are (i) that interestingly the nonlinear equation reduces to the Mathieu equation. (ii) the Mathieu limit in turn allows us to write the critical curves of the system using the Mathieu characteristic values, with excellent agreement with the numerical data. (iii) Eq. (26) enables us to compute , and consequently also [Eq. (27)], in the limit of slow external drive. This last result will be used in the following section.

## V Applications to physical systems

In this section we touch upon three physical problems to which the theory developed above can be directly applied: Rectification of microwaves in lateral semiconductor superlatticesdodin98b (); alekseev05:lsslcr (), amplification of high-frequency signals in Josephson point contacts, and modeling of Josephson junctions with critical current modulationchesca08 ().

### v.1 Semiconductor superlattices

We present the problem in the form it was introduced in Ref. alekseev05:lsslcr, . A schematic figure showing the geometry of the problem is given in Fig. 5. Plane electromagnetic wave incident on a lateral superlattice is considered. Electric field is polarized along the superlattice axis, that is, parallel to the direction of the current. The electron transport is studied by considering a single miniband with the standard tight-binding energy-quasimomentum dispersion relationwacker02:ssreview (). The electron distribution follows the Boltzmann transport equation. From there, one is able to write ordinary differential equations for ensemble averaged electron velocity and energy. These form the well-known superlattice balance equationsignatov76 (); wacker02:ssreview (). We present the equations here in their scaled form, in which the the maximum (minimum) electron velocity and energy correspond to the value (). Electric field inside the superlattice, denoted , is also appropriately scaled.

 ˙v = −uw−Γv, (28a) ˙w = −uv−Γ(w−weq). (28b)

First equation of the set describes the balance between electron acceleration by the electric field and deceleration due to scattering, while the second describes the electron energy gain and dissipation due to scattering processes. The current density is related to the average velocity by , where is the elementary charge and is the areal density of 2D electron gas. The nonlinearity is controlled by the parameter : , where is rate of electron scattering – a high density of the electron gas corresponds to a large nonlinearity, or small . Constant is the scaled equilibrium energy, whose value we set to for convenience.

The interaction of the incident electromagnetic radiation with the conduction electrons is taken into account by employing the Maxwell equations with appropriate boundary conditions. Approximating the lateral superlattice as an infinite conduction sheet, the scaled electric field entering the equations, , becomesdodin98b ()

 u=−u0cosΩt−Γ−1v, (29)

where and are the amplitude and frequency of the external electric field. With the above, Eqs. (28) are rendered nonlinear. Unlike in the case of bulk superlattices where the relation between the total electric field and the average velocity is an additional differential equationalekseev96:dissp-chaos-ssl (); alekseev98:spont-dc (), here the equation is algebraic.

The overdamped pendulum equation is obtained via a formal change of variables, , . In the physically interesting limit of , the dimensionality of the system can be reduced (see Appendix B). The following equation for is obtained:

 ˙θ+(⟨cosθ⟩Γ+Γ⟨cosθ⟩)sinθ=u0cosΩt. (30)

One of the primary interests is the appearance of a spontaneous dc voltage. The term “rectification” was used for the conversion of applied ac irradiation into a dc field and current , via the nonlinearity of these nanostructures. In terms of the variables , a prerequisite for such a current to appear is that a limit-cycle does not follow Eq. (3), i.e. . However, since the governing equations are in fact symmetric in the sense of Eq. (3), and based on the findings of Sec. III.2, the breaking of symmetry implies a pitchfork bifurcation. Consequently, rectification is expected in the real physical system when parameters are such that they correspond to the exchange of stability in the overdamped pendulum.

In the earlier workalekseev05:lsslcr (), we have considered a simplified pendulum equations in which the contribution of was ignored, that is the overdamped pendulum Eq. (1) with and . Our analytic analysis of that equation in alekseev05:lsslcr, was limited to the high driving frequency limit where the instability occurs in the vicinity of (cf. Fig. 4). Comparing with results of numerical solutions of the superlattice balance equations, we observed that the rectification indeed exists nearby the Bessel roots. Here we apply the theory developed in the previous sections to Eq. (30) in order to find the regions of instability in a wider parameter space, including the case of low driving frequencies.

We can solve the functional-differential equation (30) by considering the equation . We require that equals the coefficient of sine in Eq. (30):

 K=⟨cosθ⟩(K)Γ+Γ⟨cosθ⟩(K). (31)

Eq. (27) allows us to write in terms of , and further, Eq. (26) gives using well-known special functions. Thus, roots of Eq. (31) can be easily computed numerically.

Alternatively, parameter space structure of Eq. (30) can be studied by the following way. We introduce a simple change of variables : , where can be found following Eq. (27). Using this transformation, parameter space structure of Eq. (30) in variables can be studied for any by calculating a single dataset of values from the pendulum equation (1) with and .

For the case of this procedure was applied to plot Fig. 6, where branches of the solutions with are displayed. The pendulum limit of the superlattice balance equations [Eq. (30)] becomes invalid as , and thus we have chosen not to plot regions corresponding to . In the case of the simple overdamped pendulum, the position of, say, the th parameter space region of solutions is determined solely by the external drive amplitude and frequency [Fig. 4]. In Eq. (30) the average cosine affects values of that admit a solution corresponding to a particular . Account of the dependence has the effect of shifting the parameter space points corresponding to high values of to the higher values of , which can be seen from the change of variables . For sufficiently low , this shift is enough to make several regions with different co-existing at fixed . The overlap of the different regions of solutions can be seen in Fig. 6.

Comparing to Fig. 2 of Ref. alekseev05:lsslcr, , we find that the branches match the symmetry broken regions found using direct numerical simulation of the superlattice balance equationsisohatala09_footnote1 (). Branch in Fig. 6 falls under the branch (not plotted), which suggests that a symmetry broken region might have been missed in the earlier work. Indeed, our subsequent numerical simulations confirm that there is a region of symmetry broken solutions associated with the instability between regions and . It apparently was not detected because the initial conditions preferred the solution. This immediately demonstrates the usefulness of our analytic results.

The abstract theory developed in Sec. III gives additional insight into the dynamics of semiconductor superlattice electrons. In Sec. III we showed that the integer was connected to the number of rotations of , and consequently oscillations of periodic functions depending on : in short, counts the number of full oscillations of and across half the period of the drive. Since here and correspond to the average electron velocity and energy, respectively, we can now connect the dynamics of the charge carriers in the periodic potential to the stability of the system. At the instability the integer changes by one, and so the instability in fact marks the region where the state with oscillations of energy or current become favourable to the state with oscillations. Further, each of the branches are characterized by the number of current oscillations across half of the drive period.

In summary, the theory developed in the previous sections devoted to the mathematical analysis of Eq. (1 allowed us to analytically probe the dynamics of Eqs. (28) in the physically interesting limit of , . This limit is equivalent to very high nonlinearity, and has been largely inaccessibly analytically. Importantly, the analytical results revealed the significant multistability in this system. We found that several branches of solutions can coexist at same parameters . The analytical results also directly suggested, that a large region of rectification was missed. The region of rectification was found to be much larger than expected in the case of . Our new analytical findings are also consistent with the previous result that is required for observing rectification.

### v.2 Josephson junctions

As a second application, we consider the problem of a high-frequency gain (negative absorption) in microwave irradiated Josephson point contacts. The corresponding motion equation for the Josephson phase difference islikharev-c10-11 (); velichko99 ()

 ˙θ+sinθ=fsinωt+εcosΩt. (32)

In addition to the driving current , a probe current has been added. The probe amplitude is assumed small and the frequency is incommensurate to the drive .

Having introduced a weak probe current , we wish to find the power absorbed by the junction, and its dependence on the frequency of . The absorbed power is described by , where is the voltage across the junction. In terms of pendulum variables the absorption takes the form . For we have gain, i.e. a weak signal will be amplified. The system placed in a cavity will radiate at frequencies for which . Linearizing Eq. (32) and using the approach similar to kuzmin80, ; *kuzmin80-orig, we find

 AJJ=ε22∞∑k=−∞b−kdkΩ2Ω2+(⟨cosθ⟩+2ikω)2, (33)

where

 e±∫t0[cosθ(t′)−⟨cosθ⟩]dt′=∞∑k=−∞e2ikωt{bkdk. (34)

Derivation of Eq. (33) is presented in Appendix C.

Near exchange of stability, , the expression for absorption consists of terms of the form . These diverge when ( is an integer), and thus we expect to see strong gain for probe frequencies near an even multiple of the pump frequency when .

Again, results of the previous section can be applied to find that result in such strong response. The range of parameters where gain occurs is rather narrow, and hence an analytic first approximation is useful. In Fig. 7 we plot the regions of negative absorption using Eq. (33) together with the critical curves as given by Eq. (26). Coefficients , were computed numerically from the Floquet solutions using the fact that the expansions in Eq. (34) are equal to [cf. App. C], where is oscillating part of the unstable Floquet solution and . We also made several runs computing directly from Eq. (32). While this direct procedure is in general more slow, its results are in a good agreement with the results obtained employing Eqs. (33) and (34).

To illustrate the strong dependence of the absorption on the parameters we have also plotted typical absorption profiles, dependence of on , in the inset of Fig. 7: the profile forms strong peaks just as the critical curve for exchange of stability is approached. Gain peaks are located in the vicinities of even harmonics of the driving current. Therefore the system can radiate at frequencies close to . This effect is very different from ordinary generation of harmonics. Really, following symmetry (3) overdamped pendulum can generate only odd harmonics of the strong pump while even harmonics are forbidden. From physical point of view, generation at (exactly) odd harmonics is a spontaneous process, whereas radiation at frequencies nearby even harmonics is stimulated emission. From the viewpoint of bifurcations, the effect of gain in overdamped pendulum represents a reminiscent of symmetry breaking bifurcation in strongly damped second order pendulum Eq. (4), where the bifurcation results in appearance of even harmonics. Despite here even harmonics are forbidden by symmetry, nontrivial amplification of additional weak signal does exist at pump amplitudes and frequencies close to those necessary to realize real symmetry breaking.

Interestingly, similar gain profiles, centered near a characteristic Bloch or cyclotron frequency and its harmonics, attract much attention in physics of semiconductor superlattices, where them termed dispersive gain profiles sekine05 (); hyart09_magnetic (). Moreover, the dispersive gain profiles were found in the models of ac-driven semiconductor superlattices as wellhyart08 (); hyart09_quasistatic ().

Recently, relatively strong coherent electromagnetic radiation of technologically important terahertz frequency band has been observed from dc-driven high- superconductorsminami09 () and arrays of niobium point contactssong09 (). Physical mechanisms responsible for the observed radiation are still subject of intensive debates. Nevertheless, we notice that our simple model of a pure ac-driven junction gives a natural framework to generalization of dc-drive configurations considered experimentally, and in this respect a further development of this theory may have a beneficial applied aspect.

Finally, as a separate remark, we note that the general form of Eq. (1) with and the associated symmetry do have relevance to more complex models of Josephson junctions. Recently, a pendulum equation was used to model observed current-voltage characteristics of a type of Josephson junction with a strong applied microwave fieldchesca08 (). The magnetic field induced changes to the critical current were considered significant enough to be incorporated into the model. The resulting equation has, in the low-frequency drive case, a form corresponding to Eq. (1) and, in the absence of a dc current component, also the symmetry considered in current paper.

## Vi Conclusions

We have studied the periodically driven overdamped equation (1) with periodic coefficients by using a type of Prüfer transformation. The linear form Eq. (6) allowed for easy analysis of the nonlinear system. We showed that if the system is driven by external forcing of the form Eq. (2) there exists only one type of instability in the system. The instability was identified as an exchange of stability between two periodic solutions, and it was found to be essentially a degenerate form of a pitchfork bifurcation. We showed that the degeneracy of the bifurcation is lifted when additional terms that also follow the symmetry are added, thus proving that the exchange of stability is in a way a precursor to symmetry breaking bifurcations in realistic physical systems. Also, we found that near the exchange of stability, the overdamped pendulum responds strongly to perturbations whose frequency is close to an even multiple of the drive frequency. We further used the linear form to find explicit solutions to the problem, and used them to construct a condition for the appearance of the instability. For the simple choice of , the instabilities of Eq. (1) can be well described by the Mathieu equation. We wish to emphasize, that up to the knowledge of the authors, no such low-frequency solutions have been constructed so far.

From our purely mathematical findings there is then a natural cross-over into the field of physics. As was shown, near exchange of stability even weak perturbations can generate symmetry breaking bifurcations, or other types of strong response. Thus, regions close to this bifurcation in real physical systems are expected to show novel phenomena. We applied our methods directly the real physical problems, namely rectification of pure ac irradiation in lateral semiconductor superlattices and amplification of a weak signal in Josephson point contacts.

For lateral semiconductor superlattices, making use of our analytical methods, we were able to construct the parameter and phase space structure of this system. Excellent agreement was found with our previous numerical resultsalekseev05:lsslcr (), with the additional finding that the system exhibits strong multistability: parameter space consists of overlapping regions each carrying a solution that is characterized by the extent of oscillations of the electron gas average energy and velocity. Each of these regions was found to have an associated region of instability which in the low frequency drive and strong nonlinearity limit corresponds to rectification of the incident ac electric field – that is, the effect where the nonlinear interaction between the 2D electron gas and the external ac field generates a directed current and voltage across the superlattice. The analytical results revealed the complex way in which the multistability and the related regions of instability appear, and suggested that large regions of rectification were missed in earlier simulations. In fact, rectification persists for higher frequencies than was previously thought. The new results supported the earlier finding that fairly high electron mobilities are required for rectification to appear.

Second physical system that we considered was amplification of a weak signal in Josephson point contacts. We found (i) a direct and striking correspondence between gain, that is, negative absorption and the exchange of stability bifurcation – strong gain was shown to appear in a narrow region just around the point where loss of stability occurs, when also the weak signal frequency was close to an even multiple of the strong ac driving current. Secondly, (ii) our analytical results make it possible to find the regions of gain very accurately, or allow for computing them with little cost. Further, (iii) our findings on the dynamics of overdamped pendulum suggest insight into the physical process of amplification – the overdamped pendulum exhibits only odd harmonics which follows from the constraints set by the symmetry, yet, we found that the system is expected to radiate when the weak signal to be amplified is tuned to even harmonics.

Finally, we would like to point out intriguing similarities between gain profiles found here in the model of ac-driven Josephson junction (ac-driven overdamped pendulum) and dispersive gain profiles recently described in the models of ac-driven bulk [not lateral] semiconductor superlattices in hyart08 (); hyart09_quasistatic (). In these workshyart08 (); hyart09_quasistatic (), a superlattice is in essence modeled by two standard balance equations for electron velocity and miniband energyignatov76 (), that is, Eqs. (28) without dependence in . Taking into account some earlier findingsalekseev02a (), we speculate that bifurcations of overdamped pendulum also determine conditions for high-frequency gain in ac-driven bulk superlattices, described by the standard balance equations, in the limits of high frequencies and rare collisions. A detailed comparison of properties of amplification in ac-driven Josephson junctions with the properties of dispersive gain in ac-driven semiconductor superlattices hyart08 (); hyart09_quasistatic () in the limit of weak dissipation, however, goes beyond the scope of the present paper.

###### Acknowledgements.
We are thankful to Alex Zharov and Andrei Malkin for discussions on nonlinear dynamics in lateral superlattices, Boris Cheska and Marat Gaifullin on experimental aspects of pendulum-like dynamics in Josephson junctions, Leonid Kuzmin, Alexander Klushin, and Kazuo Kadowake – on amplification and generation of high-frequency radiation in Josephson junctions and their arrays, Timo Hyart – on correspondence between semiconductor superlattices and Josephson junctions, and Sasha Balanov for bifurcations in pendulum. We thank Feo Kusmartsev and Erkki Thuneberg for a constant encouragement of this activity. This research was partially supported by AQDJJ Programme of European Science Foundation.

## Appendix A Construction of the asymptotic solution

Here we briefly outline how the approximate asymptotic solution of Eq. (6) is constructed using asymptotic solutions of Eq. (21). Without loss of generality we can assume that . The interval is divided into subintervals where , is a turning point and is the total number of turning points. We additionally set and . Solutions on each of the intervals are approximately given by the usual WKB solutions

 y2k(t) = A2keλξ2k(t)+12B2ke−λξ2k(t)(−R(t))1/4, (35a) y2k+1(t) = A2k+1R(t)1/4sin(λξ2k+1(t)+π4) (35b) + B2k+1R(t)1/4cos(λξ2k+1(t)+π4), ξk(t) = ∫ttk√|R(t′)|dt′. (35c)

Each of the above solutions is only valid in its corresponding interval – when , solutions have the exponential form and when , the oscillatory form is the appropriate solution [cf. harmonic oscillator , real, solutions oscillate for and converge or diverge exponentially when ].

Connection formulas for the coefficients can be derived by solving the problem at the turning points. For first order zeros of , i.e. roots such that , the approximate solutions around are given in terms of the Airy functions as

 xk=A∗k[˙ϕk]−1Ai(λ2/3ϕk)+B∗k[˙ϕk]−1Bi(λ2/3ϕk) (36)

with

 ϕk=(32∫ttk[−R(t′)]1/2dt′)2/3, (37)

where , are constants. Note that above we need to raise a complex number to power , where is either real or pure imaginary. Here, the argument of is chosen to be or , i.e. so that is real. Doing so it follows that if is increasing (decreasing) around , then , , is continuous and decreasing (increasing) around . Away from the turning points the functions asymptote into the WKB solutions given in Eq. (35). This allows one to write a linear relationship between and : , where

 W2k = (2exp(κ2kλ)0012exp(−κ2kλ)), (38a) W2k−1 = (cosω2k−1λ−sinω2k−1λsinω2k−1λcosω2k−1λ). (38b)

Here, and . Denoting , we wish to construct a matrix so that . Above we have derived connection formulas for the superposition coefficients across the whole interval , so all we need are matrices that map the initial values of , , to the coefficients , and the last coefficients to the end values . In other words, we need and so that