# Multiplicity of periodic solutions in bistable equations

Gregory Berkolaiko, Michael Grinfeld

## Abstract

We study the number of periodic solutions in two first-order non-autonomous differential equations, both of which have been used to describe, among other things, the mean magnetization of an Ising magnet in a time-varying external magnetic field. When the amplitude of the external field is increased, the set of periodic solutions undergoes a bifurcation in both equations. We prove that despite superficial similarities between the equations, the character of the bifurcation can be very different. This results in a different number of coexisting stable periodic solutions in the vicinity of the bifurcation. As a consequence, in one of the models, the Suzuki–Kubo equation, one can effect a discontinuous change in magnetization by adiabatically varying the amplitude of the magnetic field.

Keywords:

## 1. Introduction

In many lattice models, the mean field approximation leads to ordinary differential equations with periodically time-dependent right-hand side. For example, the Curie–Weiss model with Glauber spin dynamics in the thermodynamic limit leads to the so-called Suzuki–Kubo equation (see the derivation in Berglund (1998, ch. 3) or in Rao et al. (1990), Suzuki & Kubo (1968), and a numerical study in Acharyya & Chakrabarti (1994)).

The Suzuki–Kubo equation is(1.1)where m is the average magnetization of a sample, ϵ is the relaxation time, is the amplitude of the applied magnetic field, and β is 1/kT, with T the absolute temperature.

A seemingly similar first-order non-autonomous equation,(1.2)has been used as a simple generic model for switchable bistable systems. If a>0 and b<0, it describes the overdamped dynamics of a particle in a quadratic potential with periodic forcing. Such equations have been studied in the context of laser optics (longitudinal mode instabilities in a semiconductor laser); see Hohl et al. (1995) and Jung et al. (1991) for analysis and references to the optics literature.

The issue we want to address is the number of 1-periodic solutions for various values of h and other parameters. In particular, we will use h as the bifurcation parameter. This issue is of physical importance, since one wants to determine whether it is possible, say, in the Suzuki–Kubo context, to effect a discontinuous change in magnetization by continuously varying the amplitude of the applied magnetic field h (i.e. a first-order phase transition).

To compare the two equations effectively, we rescale the variables in (1.1) in the following way: set, in sequence,(1.3)Then, in the new variables we have(1.4)

Comparing the two equations, we see that the Taylor series expansion of is(1.5)Hence, superficially, it looks that the behaviour of solutions of (1.4) should be similar to that of (1.2) with a=β−1 and b=−β3. In fact, this is the implicit assumption in much of the physical literature. This assumption is incorrect. It turns out that for certain values of β and h equation (1.4) has five periodic solutions, a situation which is impossible for equation (1.2). The numerical simulations confirming this were reported for example, by Berglund (1998, ch. 7). Here, we aim to establish this fact analytically.

The structure of the paper is as follows. After collecting preliminaries below, in §3 we perform the Lyapunov–Schmidt reduction of the general equation,(1.6)assuming that f is a sufficiently smooth odd function (a property shared by both (1.2) and (1.4)). The reduction is greatly complicated by the fact that we do not know the bifurcating solution explicitly.

In §4, we use the results of the reduction to conclude that for small β behaviour of (1.4) is similar to that of (1.2). Then, we consider the case of large β in (1.4) by explicitly examining the case of β=∞. The limiting equation is in fact a differential inclusion which we analyse explicitly, finding the range of h for which five solutions co-exist. Appealing to a continuation theorem (Aubin & Cellina 1984), we conclude that such behaviour persists for large finite values of β as well.

Whenever we sketch the solutions of the equation, we do it in (x,h cos(2πt)) plane, i.e. we plot the curve (x(t),y(t)), where(1.7)

This makes it easier to identify periodic solutions (they are closed curves) and produces familiar hysteresis loop pictures for large h, figure 1. In our numerics, we take ϵ of the order of 0.05 which ‘slows down’ the time and results in nicer plots.

Figure 1

The picture of the (stable) hysteresis loop for large h. In this case, we take , ϵ=0.05 and h=2.0.

## 2. Background information

The examples of function f(x) that we are going to consider here, f(x)=tanh βxx, β>1 and f(x)=γxx3, γ>0 share many important properties. Both are odd functions and the equation f(x)=0 has exactly three solutions. If we consider equation (1.6) for h=0, there will correspondingly be three stationary solutions. The central stationary solution x(t)=0 will be unstable and the other two are locally asymptotically stable.

It follows from the implicit function theorem that for small h, there will be three periodic solutions, two stable and one unstable. Below, we show that this situation persists at least until h0=maxx>0f(x). We will also show that for large h, there is only one, stable, solution. The main aim of this paper is to clarify the bifurcation picture as h is increased. The simplest scenario, where three solutions can merge into one, is a subcritical pitchfork as shown in figure 2a. There, for simplicity, we sketch the mean of periodic solutions, , versus h.

Figure 2

Two possible minimal scenarios of the emergence of a single periodic solution. Stable solutions are represented by solid lines, unstable are represented by dashed lines. In scenario (a), three solutions merge in a subcritical pitchfork bifurcation. In scenario (b), the central solution first undergoes a supercritical pitchfork bifurcation, emitting two unstable solutions. These unstable solutions disappear in fold bifurcations upon meeting the stable solutions.

Contrary to the above expectations, numerical investigations of the Suzuki–Kubo equation (1.4) performed by Tomé & de Oliveira (1990) show that for large values of β, there is an interval in h in which the equation has three stable periodic solutions (and therefore five periodic solutions altogether). The bifurcation diagram would then look like figure 2b. For small values of β, Tomé & de Oliveira (1990) reported two stable periodic solutions at most, a situation compatible with figure 2a. The critical value of β at which the bifurcation picture changes is referred to as the tricritical point (TCP).

The existence of the TCP was confirmed by Acharyya & Chakrabarti (1994) (see also Berglund (1998, ch. 3)) but disputed by Zimmer (1993) (see also Korniss et al. (2002)), who studied equation (1.2) numerically and reported that, although the stable non-symmetric solution changed fast in the vicinity of bifurcation, it was not disappearing in a ‘blue-sky catastrophe’ as in figure 2b.

Our investigation concludes that, from a mathematical point of view, both results are correct: there are at most three solutions in (1.2) and there is a TCP in equation (1.4).

When looking for periodic solutions, it is convenient to think in terms of the Poincaré map. The map establishes a correspondence between x0 and the value at t=1 of the solution x(.) of (1.6), which satisfies x(0)=x0. For both functions f(x), the Poincaré map is a continuous function defined on the whole of . The periodic solutions of (1.6) correspond to the fixed points of the Poincaré map and, most importantly, stability properties are preserved by this correspondence.

### (a) Periodic solutions for small h

Let f(x) be a continuously differentiable function satisfying(2.1)Put . Then for any h<h0, there is a unique periodic solution to equation (1.6) satisfying for all t. This solution is stable.

The idea of the proof is to show that the interval (a,b) is a trapping region: any solution that enters it must stay there. Therefore, (a,b) is invariant under the action of the Poincaré map. Then, we can appeal to the Brouwer fixed point theorem to infer that there is a fixed point inside this interval, corresponding to a periodic solution. Analysing stability of the periodic solutions in the trapping region, we find that any periodic solution in the region must be stable. But the Poincaré map, which is continuous, cannot have more than one stable fixed point without having some unstable ones. Therefore, the solution we have found is unique.

For any t, x(t)=a implies , so that any solution to the right of a will remain there. Similarly, solutions to the left of b will never increase past b. Thus, (a,b) is a trapping region and there is at least one periodic solution.

Linearizing equation (1.6) around a periodic solution x0(t) shows that it will be stable if(2.2)is negative. But since f′(x)<0 for all x∈(a,b), every periodic solution is stable. ▪

Let xm be such that , where f(x) is either tanh βx−x or γx−x3. Then for h<f(xm), equation (1.6) has exactly three periodic solutions, one in each of the intervals (−∞,−xm), (−xm,xm) and (xm,∞). The solution in the interval (−xm,xm) is unstable and the other two are stable.

The statement is a direct consequence of theorem 2.1. Indeed, starting with the interval (xm,∞), we apply theorem 2.1 with a=xm and b=y (where y is large enough to ensure f(y)<0) to conclude there is a unique periodic solution in the interval (xm,y). Since y was arbitrarily large, the uniqueness extends to (xm,∞).

The case of the interval (−∞,−xm) is completely analogous. For the middle interval (−xm,xm), we reverse the time , which results in the change and apply theorem 2.1 with a=−xm and b=xm. ▪

### (b) Uniqueness of periodic solution for large h

For large h, we have the following general theorem, which, in particular, is valid for both functions f(x) that are of interest to us.

Let f(x) be a continuous almost everywhere differentiable function. If its derivative satisfies(2.3)for some α, γ and x0, then, for sufficiently large h, equation (1.6) has a unique stable periodic solution.

We shall prove the theorem for ϵ=1, since taking any other value of ϵ results only in a trivial rescaling of the function f and constants α and γ.

Existence. Property (2.3) implies that eventually monotonically as x→±∞. This means, in particular, that the equations f(x)=h and f(x)=−h have exactly one solution each for sufficiently large h. Denoting these solutions by xh and xh+, we notice that the interval [xh,xh+] is trapping for the trajectories of (1.6). Hence, there must be at least one periodic solution.

Uniqueness. We will prove stability of every periodic solution which will imply that there can be at most one (see the proof of theorem 2.1 for more explanations). Let be a periodic solution, which is stable if λ, defined by equation (2.2), is negative. Denoting the time spent by the solution inside [−x0,x0] by t+, we can estimate λ by(2.4)We aim to show that the time t+ satisfies t+→0 as h→∞ and hence, λ<0.

To do that, we define three regions in the plane (x,h  cos(2πt)), which add up to the strip |x|<x0, see figure 3. Region A is defined by ; region B, by ; and region B′, by . For now, we leave the choice of y0 open.

Figure 3

Sketch of an orbit crossing [−x0,x0] and the three regions dividing the strip.

From the definition of region A, the time spent by the solution in A is at most(2.5)Let M be the maximum of |f(x)| on [−x0,x0]. If the solution is in B or B′, then . Therefore, the time spent by the solution in B or B′ during one visit is at most(2.6)

It is easy to see that a periodic solution can visit B at most once during one period. Indeed, by definition of B, a solution can be there only while . There is only one such interval of t during one period. To visit B more than once during one period the solution must leave and then re-enter B, while −h cos(2πt) remains above y0. The solution then must leave through the left boundary of B, see figure 3, and therefore it must re-enter through the left boundary too. But the derivative x′(t) on the left boundary of B is negative, making the re-entry impossible.

The same applies to region B′. Thus, we conclude that the total time spent in the strip |x|<x0 is(2.7)Now, if we take y0 to be equal, for example, to , the right-hand side will tend to zero as h→∞. ▪

### (c) Existence of symmetric solution when f is an odd function

In this section, we assume that function f(y) in equation (1.6) is odd. It immediately follows that if x(t) is a solution of (1.6), then so is . We aim to show that there is always a ‘symmetric solution’ that satisfies (and is therefore periodic).

Let f be a continuously differentiable odd function. Then for all h equation (1.6) has a solution x(t) satisfying .

First of all, if a solution x(t) satisfies for some t0 then and by uniqueness theorem for ODEs for all t. Now assume that retains the same sign for all t. Then is of the opposite sign. By continuous dependence on the initial condition, there exist y0 between x(0) and such that the corresponding solution y(t) satisfies . This is the solution we seek.

In the above derivation, we implicitly assume that all solutions exist for sufficiently large intervals of time. If f is unbounded, we can overcome this by restricting our attention to a trapping region, and reversing time if necessary. Such a region exists for every h precisely because f is unbounded: take a solution xh of f(x)=h and consider the interval [−xh,xh]. ▪

## 3. The Lyapunov–Schmidt reduction

We shall perform the Lyapunov–Schmidt reduction (LSR) for equation (1.6) with odd f aiming to determine the character of the bifurcation of the symmetric solution as we change the parameter h. We shall denote by hcs (the subscript ‘cs’ stand for critical symmetric), the corresponding critical value of h.

We start by defining the operator Φ by(3.1)(Sometimes, when the value of h is of no importance, we drop it from the list of arguments of Φ.) The periodic orbits are the zeros of this operator in the space of continuously differentiable 1-periodic functions contained in L2([0,1]). The reduction, based on a particular solution x0(t), leads to the construction of a reduced function such that the solutions of g(x,h)=0 are locally in one-to-one correspondence with the solutions of Φ(x,h)=0, with the solution x0(t) corresponding to the zero solution of g(x,h)=0. It is rarely possible to compute g(x,h) explicitly, but one can examine the bifurcation picture by computing the derivatives of g. For more details on the reduction, we refer the reader to Golubitsky & Schaeffer (1985).

### (a) The bifurcation condition

To start, we describe the symmetric solution x0(t) at the critical point h=hcs. The necessary bifurcation condition is that the differential of Φ has non-zero kernel. This is because if , the solution x0 of the equation Φ(x0,hcs)=0 can be uniquely continued by the implicit function theorem. In our case, the differential of Φ, which we denote by L, is(3.2)Therefore, the necessary condition for bifurcation at is that there is a non-zero 1-periodic solution to(3.3)The general solution of (3.3) is(3.4)This solution is 1-periodic, iff(3.5)which we will refer to as the bifurcation condition.

To proceed with the LSR, we need to find a basis for Ker L and (Range L). The former is spanned by the function v(t) defined by (3.4) with, for example, C=1. The latter satisfies(3.6)for all y. Solving equation(3.7)and again setting the arbitrary constant equal to one, we obtain(3.8)Note that with our choice of arbitrary constants , which will be used often in what follows.

### (b) Symmetries of the reduced function

The function g(x,h) inherits the symmetries of the original operator Φ(x). In our case, the following proposition holds (Golubitsky & Schaeffer 1985, proposition 4.3).

Let and let . If there is a symmetry operator which satisfies(3.9)for all then Rv(t) is equal to either v(t) orv(t). In the latter case, the reduced function is odd in x: .

It is easy to check that the symmetry operator satisfying conditions (3.9) is . Since, v(t) is an exponential (see (3.4)), its sign does not change over [0,1] and must be equal to −v(t). Therefore, our reduced function is indeed odd.

As a consequence, we immediately get(3.10)In addition, gx(0,hcs)=0, since hcs is critical. Thus, to investigate the pitchfork bifurcation of the symmetric solution of (1.6), if the bifurcation is indeed a pitchfork, it suffices to study gxh(0,hcs) and gxxx(0,hcs).

### (c) Bifurcation function

To study the remaining derivatives of g(x,h) at (0,hcs), we compute higher derivatives of the operator Φ (here ):(3.11)

(3.12)

(3.13)

(3.14)

We have already concluded that(3.15)(3.16)These formulae will come in handy below.

The derivative ghx is given (see Golubitsky & Schaeffer 1985) by(3.17)where E is the projection to Range L,(3.18)In our case, we have and, since , the projector E leaves invariant. Taking into account that , we get(3.19)where is the solution of(3.20)Differentiating the identity with respect to h, we notice that . Eventually, we get(3.21)Comparing with (2.2), we conclude that at a stability-gaining bifurcation we should have(3.22)

One consequence of the results of §2 is that the symmetric solution has to undergo at least one bifurcation in which it turns from being unstable to being stable. As of now, we are unable to prove that it undergoes exactly one bifurcation which is equivalent to showing that ghx is strictly greater than zero at all bifurcations.

The character of the stability-gaining bifurcation is determined by the sign of gxxx,

If gxxx(0,hcs) is positive (respectively, negative), the stability-gaining bifurcation from the zero-mean (symmetric) solution is a subcritical (respectively, supercritical) pitchfork.

We expand g(x,h) in Maclaurin series in x, taking into account that the function is odd,(3.23)To understand the bifurcation picture, it suffices to take into account only the first two terms if . Since , it remains so in a small neighbourhood of h=hcs. Therefore, the equation(3.24)has three solutions if gx(0,h)<0 and one solution, otherwise. A result connecting the sign of gx(x,h) with the stability of the corresponding solution of the original equation is given in Golubitsky & Schaeffer (1985, theorem 4.1, ch. 1). In our case, it means that the symmetric solution is stable (unstable) whenever gx(0,h) is positive (corresp. negative). Since the bifurcation is stability-gaining, we get three solutions for h<hcs and one solution for h>hcs. ▪

To study the sign of gxxx, we use the following formula(3.25)

From the definition of the projector E, see (3.18), and identity (3.16) we conclude that . Further, inverting L we get(3.26)First, we evaluate the second part of (3.25),(3.27)Setting , (3.27) can be written as(3.28)since .

Thus, in our case,(3.29)

The stability-gaining bifurcation of the zero-mean solution of (1.6) with f(x)=xx3 is a subcritical pitchfork.

In this case, and, since v2 is non-zero, gxxx>0. ▪

A different proof of this (local) theorem for cubic f is given by Byatt-Smith (2002, unpublished). That in the optical bistability equation (1.2) one can never have more than three 1-periodic solutions has been proved very elegantly by Pliss (1966).

The above results will be used in §4a to show that locally the situation in the Suzuki–Kubo equation (1.4) is similar to (1.2) if β is sufficiently small.

## 4. The Suzuki–Kubo equation

In this section, we study the case(4.1)in detail.

We aim to demonstrate that for sufficiently large β there are values of h for which there exist at least five periodic solutions of (1.4). We also show that for small β the bifurcation of the zero-mean solution is a subcritical pitchfork (figure 2a) which means that the maximum number of solutions is 3. The boundary between these two regimes is known as the TCP.

### (a) The case of small β

There exists β0 such that for all 1>β>β0 and the zero-mean periodic solution x0(t) satisfying (3.5) we have(4.2)

It is not hard to get the rough estimate β0=4/3, which is sufficient for our purposes.

Put . From (3.3), we have for all n≥1,(4.3)since v(t) is 1-periodic. Differentiating (4.1),(4.4)multiplying it by v2(s) and choosing n=2 in (4.3), we have(4.5)

For the third derivative of f(x), we now get(4.6)where we discarded the integral of a non-negative function and then used (4.5). Thus, the integral on the left-hand side of (4.6) is negative if β<4/3. ▪

It is easy to improve the above estimate to 3/2 if, instead of discarding the integral of ζ4v2 altogether, we deduce from (4.4)(4.7)

### (b) The β=∞ case

When we take β=∞, instead of the ODE (1.4) we are left with the differential inclusion(4.8)where Sgn(x) is the multivalued sign function,(4.9)

For this case, we have the following theorem:

There are values hcs<hca such that for h satisfying hcs>h>hca, the inclusion (4.8) has at least five 1-periodic solutions.

This assertion is verified by explicit computations, which involve constructing solutions of piecewise-linear ODEs. We explain the main steps of the constructions involved.

First of all, note that if a periodic solution x(t) has no internal zeroes, it solves a simple linear ODE and can be computed explicitly. For example, for h sufficiently small, a positive periodic solution is given by(4.10)from which we conclude that this (asymmetric) solution touches zero when(4.11)For example, if ϵ=0.05, hca=1.048187027, to 10 d.p. Periodic solution (4.10), as well as its negative counterpart, are clearly asymptotically stable. To verify theorem 4.3 it is enough to show that there is another asymptotically stable periodic solution for some values of h<hca. Indeed, three stable solutions would mean five periodic solutions altogether. The key is lemma 4.4.

Consider the differential equation(4.12)If(4.13)where(4.14)then equation (4.12) has a solution x0(t) satisfying, for some t0,(4.15)and(4.16)

Gluing two such solutions together, we obtain a symmetric periodic solution similar to the one sketched in figure 4d. Both parts of this solution attract nearby solutions, therefore it is asymptotically stable.

Figure 4

Evolution of the zero-mean solution of inclusion (4.8). (a) Solution is identically zero when h<1. (b) Growth of nontrivial parts when 1<h<hcs. (c) The solution experiences a bifurcation at h=hcs when the trivial parts, where x(t)=0 on an open interval, disappear completely. (d) The stable solution for h>hcs contains no trivial parts.

The proof of the lemma is based on the explicit integration of (4.12) and analysis of the solution satisfying the initial condition(4.17)In appendix A, we explicitly compute the value of h for which this solution satisfies condition (4.15); this value is the critical value hcs. It can then be shown that when h>hcs the solution satisfying (4.17) has another zero at t1<t0−1/2 and remains positive on the interval (t1,t0). Existence of the solution of lemma 4.4 then follows from the continuous dependence on the initial conditions.

By studying , one concludes that hcs<hca for all ϵ. Furthermore, expanding this expression in z, we obtain(4.18)from which it is clear that for small ϵ (and hence for very small z) hcs and hca are very close indeed. For example, if ϵ=0.05, hcs=1.048091900 and therefore hcahcs=O(10−4).

Figure 4 explains the bifurcation picture1 of the zero-mean solution. When h<1, the zero-mean solution is identically zero and is unstable. After h crosses 1, the solution gets two symmetric non-trivial parts, but there are still two intervals in t on which x(t) is identically zero. The bifurcation happens when each of these intervals collapses to a point.

When h=hcs, the zero mean solution becomes stable through a super-critical pitchfork bifurcation. It emits two unstable solutions which disappear in a blue-sky catastrophe, i.e. in a collision with the stable solutions (4.10). These unstable solutions are sketched in figure 5 for various values of h, hcs<hhca.

Figure 5

The evolution of the right unstable solution until its annihilation at h=hca.

### (c) The case of large β

It is possible to extend the existence of five solutions from β=∞ case to the case of large β by continuity. In our numerical experiments, we found that β does not have to be very large. Figure 6 contains the plot of three stable solutions of equation (1.4) with β=3, ϵ=0.05, h=0.544 and initial conditions x(0)=−0.413600598, −0.40456237, −1.491606074.

Figure 6

Three coexisting stable solutions of the Suzuki–Kubo equation for β=3, ϵ=0.05 and h=0.544.

For any h satisfying hcs<h<hca and sufficiently large (depending on h) β equation (1.4) has at least five periodic solutions.

The proof of this theorem is based on an upper semi-continuity result for differential inclusions (Aubin & Cellina 1984). To state it, we need to introduce some notation. Let B be the unit ball in n with centre at the origin. Consider a differential inclusion(4.19), where Ω is a bounded subset of n, Λ and F is an upper semi-continuous set-valued map with uniformly bounded compact images. Assume that for all λΛ and any x0 in some set QΩ the solutions of the differential inclusion with x(0)=x0 exist on the interval [0,T] and remain in Ω. The Poincaré map is the set-valued map defined by(4.20)The following theorem is a slight variation of the continuity theorems in Aubin & Cellina (1984, section 2.2).

Let the map satisfy the conditions of the previous paragraph. Then for any λ0Λ and any open set U2n, 0∈U, there is a δ>0 such that |λλ0|<δ implies(4.21)

Clearly, F(t,x,β) defined by and viewed as a differential inclusion satisfies the conditions of theorem 4.6 since all solutions of this inclusion are globally defined and remain bounded.

Fix h∈(hcs,hca). By theorem 4.6, for sufficiently large β the graph of the Poincaré map of the Suzuki–Kubo equation (1.4) is contained in a neighbourhood of the graph of . By theorem 4.3, the graph of has five intersections with the diagonal. Therefore, for sufficiently small neighbourhood, i.e. for sufficiently large β, we get at least five 1-periodic solutions of (1.4). ▪

## 5. Conclusions

In this paper, we have resolved a controversy in the physical literature. In particular, we have proved that for large values of β and certain values of the parameter h there are three stable periodic solutions to the equationThe width of the h-interval in which three stable solutions coexist is proportional to exp(ϵ−1/2), which makes it very hard to detect for small values of ϵ. It is remarkable that the numerical simulations of Tomé & de Oliveira 1990 and Acharyya & Chakrabarti (1994) revealed this interval without a priori knowledge about its existence.

However, some major mathematical questions still remain. One is to find the conditions on f which would guarantee the existence of five periodic solutions for some values of h. We would like to formulate another question as a conjecture.

If the ODE dx/dt=f(x) has precisely two stable rest points, an ODE of the form of (1.6) has at most five 1-periodic solutions.

## Acknowledgments

The authors gratefully acknowledge fruitful discussions with J. Byatt-Smith, B. Duffy and J. Carr.

## Footnotes

• A complete bifurcation picture and animated periodic solutions described in this section are also available at http://www.math.tamu.edu/∼berko/bistable.

• Received June 1, 2005.
• Accepted October 31, 2005.

View Abstract