Robust homoclinic orbits in planar systems with Preisach hysteresis operator

We construct examples of robust homoclinic orbits for systems of ordinary differential equations coupled with the Preisach hysteresis operator. Existence of such orbits is demonstrated for the first time. We discuss a generic mechanism that creates robust homoclinic orbits and a method for finding them. An example of a homoclinic orbit in a population dynamics model with hysteretic response of the prey to variations of the predator is studied numerically.


Introduction
The theory of hysteresis operators developed in [1][2][3] as well as other mathematical descriptions of hysteresis (such as differential inclusions [4], differential inequalities [5] and variational descriptions of rate-independent processes [6]) allow one to analyze stability [7,8], bifurcations [9][10][11][12] and complex dynamics [13] of systems exhibiting hysteresis. Coupled systems of differential equations and hysteresis operators have been used to model mechanical [2,14,15], electro-magnetic [16,17], electro-mechanical [18,19] and engineered systems [20][21][22], phase transition processes [3,23], dynamics of fluids in porous materials [24,25], dynamics of populations [26][27][28] and economic systems [29][30][31]; see [32] for further examples of applications. In particular, these models include ordinary differential equations coupled with the Preisach hysteresis operator that contain the time derivative of the output of the Preisach operator. Although more complex situations are possible, a generic solution of such a system can be obtained by matching solutions of several ordinary differential systems (the switching points from one to another system include the turning points of the input of the Preisach operator as well as some values stored by the dynamic memory configuration of this operator) [33][34][35]. While these ordinary differential systems are different, they all have the same set of equilibrium points [36].
This picture suggests a possibility of existence of robust homoclinic trajectories in the system with hysteresis operator provided that (a) the same point is a stable equilibrium of one ordinary differential system and, simultaneously, an unstable equilibrium of the other; and, (b) a trajectory of the first system converging to this equilibrium point in forward time could be matched with a trajectory of the second system converging to the same point in backward time in such a way that the system with hysteresis would switch from the second trajectory to the first trajectory. The concatenation of the two trajectories would then form a homoclinic loop of the differential system with the hysteresis operator.
In this work, we show the existence of such homoclinic trajectories. To the best of our knowledge, homoclinic trajectories of systems with hysteresis operators have not been discussed before.
Existence of homoclinic trajectories in the systems with Preisach operator under the time derivative is closely related to the so-called partially stable equilibria [36], i.e. unstable equilibria, which have an open basin of attraction (in the infinite-dimensional phase space of the system). Such equilibria simultaneously attract and repel many trajectories, and they can be compared to a saddle-node singular point of an ordinary differential system, however they are robust. Systems with Preisach operator available in the literature do not demonstrate abundance of partially stable equilibria, and, moreover, it is even harder to find special kind of partially stable equilibria that would lead to the appearance of robust homoclinic orbits in such systems. We present our results in the context of mathematical ecology. Namely, we are motivated by a predator-prey system [37], where the prey switches between two modes of behaviour, safe and risky, in response to varying abundance of predator. Similar type of non-hysteretic safe-risky adaptive behaviour of the prey in response to the pressure of the predation was recently suggested as a possible mechanism leading to the so-called predator pit phenomena, where a bistability between two steady states with nonzero predator-prey densities appears [38]. When the switch between two modes of behaviour becomes hysteretic, one can observe numerically the existence of a robust homoclinic cycle [37].
We first construct a few examples of simple systems that have homoclinic orbits (Section 2). Then we give a numerical evidence of the existence of such orbits in a more complex predator-prey model with a refuge patch where we assume hysteresis in the reaction of the prey to variations of the predator abundance (Section 3). Unlike the familiar homoclinic orbits of ordinary differential systems, the homoclinic trajectories that we obtain are robust, that is they are preserved under variations of parameters and more general (sufficiently small) perturbations. Moreover, our systems typically have a continuous family of homoclinic orbits.

Examples of homoclinic orbits 2.1. Preisach operator
In this work we consider a coupled system of differential equationṡ and an operator equation where t ≥ t 0 , and f, h, g : The Preisach operator is defined by where v = v(t), t ≥ t 0 is the input; the function η 0 = η 0 (α R , α S ) which takes values 0 and 1 is the initial state function; µ(α R , α S ) is an integrable density function; and R α R ,α S is the non-ideal relay operator with thresholds α R , α S satisfying α R ≤ α S : We represent states (4) of the relays graphically on the half-plane {(α R , α S ): α R ≤ α S }, which is divided into two parts by a staircase polyline Ω = Ω(t) with the relays in state 1 below (to the left) of this line and in state 0 above (to the right) of this line, see Fig. 1 (left). At any moment t ≥ t 0 the right end point of Ω(t) is the point α R = α S = v(t) [1]. For simplicity, the measure µ(α R , α S ) is often assumed to have triangular support α m ≤ α R ≤ α S ≤ α M , hence we are interested only in the representation of the states of the relays inside this triangle.
The staircase polyline Ω(t) changes over time in accordance with formula (4) in response to variations of continuous non-negative input v(t). We illustrate the evolution of Ω(t) with time t ≥ t 0 using the following example with α m = 0. Let Ω(t 0 ) be the horizontal segment Fig. 1, right). If v(t) monotonically decreases on an interval t ∈ [t 0 , t 1 ], then Ω(t) acquires the vertical link α R = v(t), v(t) ≤ α S ≤ v(t 0 ) as shown in Fig. 2, left. Now, suppose that the input monotonically increases for t ∈ [t 1 , t 3 ]. Then Ω(t) acquires a new Fig. 2, right). If the input reaches its initial value v(t 0 ) at some moment t * 2 , then Ω(t * 2 ) coincides with the initial horizontal segment Ω(t 0 ) ( Fig. 1, right). For t ∈ [t * 2 , t 3 ] the input further increases and Ω(t) remains horizontal until v(t) reaches a maximum v(t 3 ) (Fig. 3, left). If v(t) decreases after the moment t 3 , then Ω(t) acquires a vertical link, see Fig. 3, right. This example demonstrates how the staircase polyline Ω(t) can acquire and loose multiple vertical and horizontal links in response to variations of the input v(t); for more details see [13]. Below Ω(t) is referred to as the polyline of the state η(t) of the Preisach operator or, simply, the state.
In equations (1), (2), the derivative of the output of the Preisach operator is used. For the evaluation of this derivative, the most right link Ω e = Ω e (t) which is attached to the right end point α R = α S = v(t) of the staircase polyline Ω(t) is important, see Fig. 4, left (if Ω has infinitely many links, then Ω e = ∅). Denote by (v m , v), (v, v M ) the end points of the segment Ω e , where v m = v if Ω e is a vertical segment and v M = v if Ω e is horizontal. If v = v(t) increases, then the time derivative of the output of the Preisach operator satisfies [33] If v decreases, then    Figure 2. Evolution of the staircase polyline Ω = Ω(t) from the initial state Ω(t 0 ) shown in Fig. 1 (right) in response to an input v(t). For t 0 < t ≤ t 1 the input v(t) monotonically decreases, and the line Ω consists of two segments (left panel,

Homoclinic solutions
connects to the horizontal link, which is a part of the segment Ω(t 0 ) shown in Fig. 1 , and the state Ω has three links as shown on the right panel for the moment t = t 2 . The leftmost horizontal link and the vertical link of the line Ω(t 2 ) are parts of the staircase Ω(t 1 ) presented on the left panel, and the rightmost horizontal link is . The state Ω(t * 2 ) becomes the horizontal segment (which coincides with Ω(t 0 )) and remains the horizontal segment that is, as long as the input increases), as shown on the left panel. After the moment t 3 the input v(t) decreases again, and the state Ω(t 4 ) presented on the right panel for a moment t 4 > t 3 is similar to the staircase in Fig. 2 (left). .
Assume that system (1), (2) has an equilibrium. Without loss of generality, assume that u * = v * = 0. Since the output of the Preisach operator (2) is defined on a semiaxis t ≥ t 0 , we use the following definition of a solution to system (1), (2) on the whole axis t ∈ R.
In this work, we focus on the piecewise monotone homoclinic solutions {u(t), v(t), η(t)} with monotone tails.
For such trajectories we can define the limit of the state of the Preisach operator at t → ±∞. We will write η(t) → η * as t → −∞ (or, as t → ∞) if the polyline Ω(t) of η(t) converges to the polyline Ω * of η * (as a set on the plane (α R , α S )).

Existence of robust homoclinic orbits
To demonstrate rigorously the possibility of existence of homoclinic orbits, we consider the following example of system (1):  Figure 4. Evolution of the polyline Ω(t) for an input v(t) that has one maximum point v(t 1 ) = v 1 ; Ω e (t) is the segment of this staircase polyline with the end point α R = α S = v(t) on the bisector. The limit state Ω * (solid line on the left panel) in this example consists of a horizontal segment and the vertical segment Ω * e . During the time interval (−∞, t 1 ] when the input v(t) increases, Ω(t) has three segments, the segment Ω e (t) is horizontal (thick dashed line). The right panel shows the state Ω(t 1 ) (solid line) at the moment when the input achieves its maximum v(t 1 ) = v 1 . After this moment, the input decreases, Ω(t) consists of four segments, the segment Ω e (t) is vertical (thick dashed line). and assume that the density function of the Preisach operator (2) is uniform, µ(α R , α S ) ≡ 1. This system has an equilibrium at zero, u = v = 0.
First, we prove the existence of a simple homoclinic orbit with the limit state η * shown in Fig. 4 (left) such that the polyline Ω * has a sufficiently long vertical segment Ω * e = {(α R , α S ) : α R = 0, α S ∈ [0, v M ]}. We aim to construct a homoclinic trajectory with a component v(t) that increases for t < t 1 , decreases for t > t 1 and has one maximum point v(t 1 ) (cf. Assumption 1), see Fig. 5.
κc > 0, and the segment Ω * e of the polyline Ω * of the state η * is vertical with the end points which is a point of maximum. That is, system (2), (7) has infinitely many robust homoclinic trajectories satisfying η(t) → η * as t → ±∞.
Robustness means that the continual set of homoclinic trajectories persists under perturbations of system parameters.
Proof. A homoclinic solution of system (2), (7) described in this proposition satisfies relation (5) with v m = 0 for t < t 1 and relationship (6) with v M = v 1 for t > t 1 . Hence, for such a solution, x = vv for t < t 1 andẋ = (v 1 − v)v for t > t 1 (where we use the relationship µ(α R , α S ) ≡ 1 when calculating H(v, 0) and V (v, v 1 )). Substituting these expressions forẋ in (7), we see that on the semiaxis t < t 1 ; and, the linear differential systeṁ on the semiaxis t > t 1 ; the point (u(t 1 ), v(t 1 )) = (u 1 , v 1 ) lies on the nullclinev = cu + dv = 0 of these systems, as v 1 is the maximum value of the v-component. Furthermore, to prove the proposition it suffices to show that for every : (a) The solution of system (9) with the initial conditions u(t 1 ) = −dv 1 /c, v(t 1 ) = v 1 satisfieṡ v > 0 for all t < t 1 and v → 0 as t → −∞; (b) The solution of system (10) with the same initial conditions u(t 1 ) = −dv 1 /c, v(t 1 ) = v 1 satisfiesv < 0 for all t > t 1 and v → 0 as t → ∞.
Here (a) follows from conditions (8), which ensure that the zero equilibrium of system (9) is a proper unstable node; (b) follows from the relationship v 1 > v 1m , which implies that the zero equilibrium of system (9) is a proper asymptotically stable node.
The case κc < 0 can be considered similarly. In this case, the following analog of Proposition 1 is valid.
Proposition 2. Assume that conditions (8) and κc < 0 hold and the segment Ω * e of the polyline Ω * of the state η * is horizontal with the end points (α R , α S ) = (0, 0) and there is a homoclinic trajectory with min t∈R v(t) = v 1 . That is, system (2), (7) has infinitely many robust homoclinic trajectories satisfying η(t) → η * as t → ±∞ such that the v-component of the trajectory has a single extremum point, which is a point of minimum.
As the last example of this section, we consider homoclinic trajectories which have more complex limit states η * at t → −∞, see in Fig. 6, and more than two monotonicity intervals of the v-component, see Fig. 7.
Proposition 3. Assume that conditions (8) hold, κc > 0 and the polyline Ω * of state η * includes three links connecting the corner points (0, 0), Then system (2), (7) has a robust homoclinic trajectory with the limit state η(t) → η * as t → −∞. The v-component of this trajectory has a single maximum point and a single minimum point with max t∈R v(t) = v * M and min t∈R v(t) < v * m . We note that the homoclinic solution described in this proposition has different limit states at t → ±∞, that is lim t→∞ η(t) = lim t→−∞ η(t) = η * . As a matter of fact, there are again infinitely many homoclinic trajectories with the maximum value v 1 = max t∈R v(t) of the vcomponent ranging over the interval v 1 ∈ (0, v * M ]. We consider one of these trajectories as an example.
Proof. The proof is similar to that of Proposition 1. Here, it suffices to show that: with v 3 = v min starting from the initial point (u(t 3 ), v(t 3 )) = (−dv min /c, v min ) on the nullclinė v = cu + dv = 0 satisfiesv > 0 for all t > t 3 and v → 0 as t → ∞. (8), which ensure that the zero equilibrium is a proper unstable node of system (9).

Statement (a) follows from assumptions
The vector fields (u,v) of systems (9) and (10) coincide on the nullclinev = 0. Asv > 0 for t < t 1 for the solution of (9) with the initial condition (u(t 1 ), v(t 1 )) = (−dv * M /c, v * M ) on the nullcline, it follows that the solution of (10) starting from the same initial condition satisfiesv < 0 for sufficiently small t − t 1 > 0. This solution of (10) remains in the half planev = cu + dv < 0 for all t > t 1 and satisfies v → −∞ as t → ∞, because the assumption v * M < v 1M ensures that the zero equilibrium is a proper unstable node of system (10) with v 1 = v * M . As v * M > 0 > v * m , this solution will eventually hit the line v = v * m at some moment t 2 at a point (u m , v * m ). This proves (b).
The assumption v * M 2 > v 1m implies that the zero equilibrium of system (10) with is a proper asymptotically stable node. As v * m < 0 and the point (u m , v * m ) lies in the half planė v = cu + dv < 0, it follows that the solution of this system starting at (u m , v * m ) will hit the nullclinev = cu + dv = 0 at some point (−dv min /c, v min ), that is (c) holds.  Finally, (d) follows from the relationships v min < v * m < −v 1m < 0, which ensure that the zero equilibrium is a proper asymptotically stable node of system (11) with v 3 = v min .

Population dynamics model
In this section, we give a numerical evidence of the existence of homoclinic trajectories in a predator-prey type model [37], which has nontrivial dynamical properties such as multiple stable equilibria:u where u is the number of prey; v is the number of predator; the term a(u) = ρu − λu 2 describes the logistic growth of the prey with the birth rate ρ and the competition rate λ; is the Holling type II functional response; is the predator interference; σ is the efficiency of conversion of food to growth; and, the term describes death of the predator with the death rate γ (all the parameters are positive). The term h(t) describes the flow of the prey to some refuge patch from which it never returns. System (12)-(13) is analagous to the two-patch predator-prey system proposed in [37], where the safe patch is completely separated from the environment; alternatively, the term h(t) can be considered as an additional predator-induced death rate. We consider the flow rate in response to the change of predator abundance in the following form: where we ensure that the flow satisfies h(t) ≥ 0 using the function x + = max{x, 0}, andv can be replaced with the right hand side of equation (13). The main ingredients here are the constant flow rate k 0 and the hysteretic reaction of the prey to variations of the predator abundance, k d dt (P[η 0 ]v)(t), where hysteresis appears due to a delay of the response of the prey to a change of the trend in predator dynamics (i.e., the change of the sign ofv). A specific choice of the model for this hysteretic response in the form of the Preisach operator is related to a number of phenomenological assumptions. Namely, the environment is divided into many small interconnected patches, and it assumed that the change of the rate of flow of the prey to the refuge from any given patch is delayed until the abundance of the predator drops/increases from its extremum value by a certain positive amount, which is specific to the patch. More detailed derivation of the model is presented in [37]. We assume that the measure µ(α R , α S ) of the Preisach operator P[η 0 ] is zero outside the triangle 0 ≤ α R ≤ α S ≤ 1, and µ(α R , α S ) = 2 inside the triangle. The additional term 2k † |v|(v − v † ) in (14) strengthens the reaction of the prey to the change of the predator abundance if v is below the threshold v † , and weakens this reaction if v is above the threshold. In other words, the prey is attentive to the change of the predator numbers if the predator is not abundant and its attention to the predator becomes saturated for v ≫ v † . This artificial term is similar to the term we introduced in (7) and, indeed, allows us to demonstrate numerically dynamics of (12)- (14), which are similar to dynamics of equations (2), (7). We will consider the system with (k † > 0) and without (k † = 0) this term.
In the case of increasing v, substituting formula (5) in equations (14), we obtain the following expression for the rate of flow to the refuge: wherev can be replaced with the right hand side of equation (13). Similarly, when v decreases, Equilibria of system (12)- (14) can be found from the algebraic system, which is obtained by setting the derivatives of all the variables, including d/dt(P[η 0 ]v), to zero in (12)- (14).
Numerical examples of homoclinic orbits of system (12)- (14) presented in the next section were obtained by an algorithm, which is similar to the method used in the previous section.
If k = k † = 0 (the rate to the refuge (14) is constant, there is no hysteresis), then these equilibrium points of the ordinary differential system (12)- (14) have the eigenvalues (0.136, 0.614), (−0.089, 0.942), and (−1.33, −0.25), respectively. That is, the first equilibrium is an unstable node, the second equilibrium is a saddle and the third equilibrium is a stable node. When the hysteresis term is present (k > 0), we give a numerical evidence that equilibrium (17) can have a homoclinic orbit attached to it. We consider a constant density function of the Preisach operator (3), in the triangle 0 ≤ α R ≤ α S ≤ 1 and set µ = 0 outside this triangle. The integral of µ over the whole half plane α S ≥ α R is normalized to 1.
As the limit state of the Preisach operator, lim t→−∞ η(t) = η * , we choose the polyline Ω(t 0 ) which has two links: a vertical link Ω e (t 0 ) = {(α R , α S ) : where v * is the second component of equilibrium (17) (see Fig. 4, left, where the origin is shifted to the point α R = α S = v * ).
Next, we construct a homoclinic trajectory similar to the one shown on Fig. 7 for system (12)- (14) with k † = 0, k = 700. The parameter of the limit state is set to v * M = v 1 = 0.2596. As in the previous example, we solve system (12)- (13), (15) with v m = v * numerically backward in time starting from the point (u 1 , v 1 ) on the nullclinev = 0 with v(t 0 ) = v 1 , see Fig. 9. This solution satisfiesv > 0 for t < t 0 . Then, we solve system (12)- (13), (16) with v M = v 1 for t > t 0 starting from the same point (u 1 , v 1 ). We observe that v(t) decreases, while the trajectory spirals around the equilibrium (17) for t 0 < t < t 2a , where t 2a ≈ t 0 + 1.612, because (17) is a stable focus for this ordinary differential system. However, for t = t 2a , v(t 2a ) ≈ 0.25741 the flow h(t) in (12) becomes 0; at this point, the term kvV (v, v 1 ) is negative, and |v| becomes large enough to compensate for the term k 0 > 0 in relation (16). We observe that kvV (v, v 1 ) + k 0 < 0 and h(t) ≡ 0 for t 2a < t < t 2b , where t 2b ≈ t 0 + 11.043, v(t 2b ) ≈ 0.0384, and kvV (v, v 1 ) + k 0 = 0 at t = t 2b . After this moment, h(t) becomes positive again. Shortly after this point, the vcomponent of the solution reaches its (global) minimum v(t 3 ) = v 3 ≈ 0.0384 at t 3 ≈ 11.046 as the solution arrives at the nullclinev = 0. From the nullcline, we continue the trajectory for t > t 3 as the solution of ordinary differential system (12)- (13), (15) with v m = v 3 . For this system, the equilibrium (17) (17) is denoted by A.

Conclusion
We have demonstrated the existence of robust families of homoclinic trajectories for operatordifferential equations (1), (2) with the Preisach hysteresis operator. These systems can be considered as switching systems where switching from one planar ordinary differential equation to another occurs when a variable either passes a turning point or achieves a value stored in the memory state of the hysteresis operator. Homoclinic trajectories are possible when the basin of attraction of an equilibrium of one planar system in forward time overlaps with the basin of attraction of the equilibrium of another planar system in reversed time.