The dynamics of a driven harmonic oscillator coupled to independent Ising spins in random fields
TThe dynamics of a driven harmonic oscillatorcoupled to independent Ising spins in random fields
Paul Zech, Andreas Otto, and G¨unter Radons
Institute of Physics, Chemnitz University of Technology, Germany. (Dated: June 25, 2020)We aim at an understanding of the dynamical properties of a periodically driven damped harmonicoscillator coupled to a Random Field Ising Model (RFIM) at zero temperature, which is capableto show complex hysteresis. The system is a combination of a continuous (harmonic oscillator)and a discrete (RFIM) subsystem, which classifies it as a hybrid system. In this paper we focuson the hybrid nature of the system and consider only independent spins in quenched random localfields, which can already lead to complex dynamics such as chaos and multistability. We study thedynamic behavior of this system by using the theory of piecewise-smooth dynamical systems anddiscontinuity mappings. Specifically, we present bifurcation diagrams, Lyapunov exponents as wellas results for the shape and the dimensions of the attractors and the self-averaging behavior of theattractor dimensions and the magnetization. Furthermore we investigate the dynamical behavior ofthe system for an increasing number of spins and the transition to the thermodynamic limit, wherethe system behaves like a driven harmonic oscillator with an additional nonlinear smooth externalforce.DOI: 10.1103/PhysRevE.101.042217
I. INTRODUCTION
This work is motivated by problems, which arise if adynamical system contains a hysteretic subsystem. Hys-teresis phenomena in general can be found in many dif-ferent research fields, such as magnetic effects at oxideinterfaces [1], shape memory alloys [2], ultrathin single-layer films [3], organic ferroelectrics [4], soft porous crys-tals [5], atomtronics [6] and metallic nanoparticles [7].An overview over the field of hysteresis can be found in[8]. In general, hysteresis means that the instantaneousoutput depends not only on the current input value butalso on its history. Thus, systems with hysteresis are sys-tems with memory. For example, in the case of magneticmaterials this means that the magnetization and the ori-entation of the magnetic domains depend not only on thecurrent external magnetic field but also on its past behav-ior. In contrast to simple bi-stability, complex hysteresisis characterized by multistability, i.e. multiple internalstates are possible for a single input value, and non-localmemory, i.e. different internal states are connected to agiven output value. As a consequence, not only one ma-jor loop but also various small subloops may appear insystems with complex hysteresis.One of the most prominent model for complex hystere-sis is the Preisach model [9]. It is a purely phenomeno-logical model and a superposition of rectangular hystere-sis loops, which are the elementary building blocks, alsocalled Preisach units. In contrast to the Preisach modela more physical way to model hysteresis is the zero-temperature RFIM [10, 11]. In the Ising model hysteresisappears because of the interaction between spins, whichrepresent, for instance, magnetic or dielectric dipole mo-ments of atoms. This paper serves as preparatory workfor studies of the zero-temperature RFIM initially estab-lished to study phase transitions with a renormalization group approach [12]. In addition to the usual Ising modeleach spin in the RFIM has its own local quenched disor-der field, which in general leads to ”smooth” hysteresisloops instead of ”hard” jump-like loops appearing in thenormal Ising model. The RFIM shows many properties,which can be also found in the Preisach model [13], but incontrast to the Preisach model, the RFIM is a spatiallyextended model.Typically, the dynamical interaction of a hystereticsubsystem (hysteretic transducer) with some environ-ment can be considered in two different ways. In thefirst scenario, the hysteretic transducer handles the in-put from the environment and produces the output ofthe system without feeding back to the environment. Incontrast, in the second scenario the feedback of the hys-teretic transducer to the environment is not negligible. Inthis case, a dynamical model, e. g. in form of an OrdinaryDifferential Equation (ODE), is necessary for describingthe environmental behavior. Many results in the litera-ture on hysteretic systems can be attributed to the firstscenario [13–20]. Other studies without feedback focuson thermal relaxation processes [21, 22], or the powerspectral density of stochastically driven Preisach models[23–28]. On the other hand, not much is known aboutthe second scenario, i.e. hysteretic systems coupled toits environment via a feedback mechanism. The gen-eral difficulty for such problems lies in the model for thehysteresis and the resulting dynamical systems. For ex-ample, for the Preisach model or the RFIM one obtainscoupled ODE-Preisach-operator equations or piecewisesmooth hybrid dynamical system, respectively. Somework has been done on ODE’s coupled to a Preisach op-erator. The ferroresonance phenomenon in LCR circuitswith an inductance modeled by a Preisach operator isstudied in Refs. [29, 30], and the mechanical equivalent,an iron pendulum in a magnetic field has been studied a r X i v : . [ n li n . C D ] J un in Ref. [31], where the hysteresis appears because of theinteraction between the ferromagnetic iron mass and themagnetic field.In a general manner, we are interested in such dynam-ical systems with hysteresis. As a prototypical examplewe consider a driven harmonic oscillator similar to [31],but in contrast to [31], the dynamics of the magnetiza-tion of the iron pendulum is modeled by a bulk of Isingspins.As a first step, especially in this paper, we neglectspin-spin interactions and we will focus on systems withnearest neighbor interactions in following papers. Theabsence of spin-spin interactions means, that the systemdoes not have non-local memory and no hysteresis be-tween the intensity of the magnetic field and the magne-tization of the pendulum is possible. However, alreadythis simplified system of a pendulum coupled to a RFIMwith independent spins shows very complex behavior. Onthe one hand side the dynamics of the system is quite in-teresting, because of the hybrid character of the systemwith discrete internal states of the Ising spins and a con-tinuous nature of the pendulums motion. Such kind ofsystem are called piecewise-smooth hybrid systems andcan be found in many fields and in every system, where asudden change of the dynamics appear. Typical examplesare relay feedback systems [32–34] or mechanical systemswith strong impacts, such as impact moling, ultrasonicassisted machining [35], gear dynamics including back-lash [36], or systems exhibiting dry friction [37, 38]. Anactual overview over this topic can be found in [39]. Onthe other hand we are indeed neglecting the memory andtherefore the hysteresis in the system, but the disorder ofthe random fields of the Ising spins can cause some inter-esting issues when dealing with physical properties of thesystem because of the dependency on the actual disorderrealization. Therefore questions of self-averaging arise.The paper is organized as follows. In Sec. II we give abrief introduction to our model which basically consistsof two parts: an oscillator model and an Ising model. Be-cause of the hybrid character of the system, in Sec. III webriefly introduce specific methods for piecewise-smoothsystems, derive the thermodynamic limit in the case ofan infinite number of spins and make some remarks onthe numerical calculation of the trajectories and the Lya-punov exponents. In Sec. IV we present characteristic re-sults on the dynamics of the system with one spin as wellas the system with many spins, the transition from thepiecewise-smooth hybrid system to the smooth system inthe thermodynamic limit and the self-averaging behaviorof the attractor dimensions and the magnetization. Weend with a conclusion and an outlook on future work inSec. V. II. MODEL
Our dynamical system basically consists of two sub-systems. One part is the continuous subsystem given by c k BA cos ωtx F M FIG. 1. Illustration of the prototypical example of a dynam-ical system with external force F M (red curve), because ofthe interaction between the iron mass (black disk) and theexternal magnetic field (grey arrows). Here the external forceshows non hysteretic behavior because of the neglected near-est neighbor interaction. a periodically driven damped harmonic oscillator as de-scribed in Sec. II A. The second part is the discrete sub-system eventually given by the full RFIM and describedin Sec. II B. A mechanical example for such a system,which can be realized in experiments, is illustrated inFig. 1. A. Oscillator model
We consider a periodically driven damped harmonicoscillator with an iron mass subject to an external mag-netic field (see Fig. 1). The position of the iron mass x ( t ) can be determined by [31]: m x (cid:48)(cid:48) ( ϑ ) + c x (cid:48) ( ϑ ) + k x ( ϑ ) = A cos ωϑ + F M ( ϑ ) , (1)where m , c and k are the mass, damping and stiffnessof the oscillator and A and ω are the amplitude and theangular frequency of the periodic excitation. F M denotesthe additional force, which comes from the interaction ofthe iron mass with the external magnetic field and isdescribed in detail below.The oscillator model Eq. (1) is put into a dimension-less form by using the transformations t = (cid:112) k/mϑ and q ( t ) = kA x ( ϑ ). Thus, we obtain¨ q ( t ) + 2 ζ ˙ q ( t ) + q ( t ) = cos Ωt + F ( t ) , (2)where ζ denotes the damping ratio, Ω is the dimension-less excitation frequency and F ( t ) = F M ( ϑ ) A . Note thatthe angular eigenfrequency of the oscillator in Eq. (2) isequal to one, which means that the resonant forcing isgiven by Ω = 1 and the period of the resonant oscilla-tion is equal to 2 π . Later we will see, that the additionalforce will be piecewise constant F ( t ) = const = CM ,which allows us to give an explicit solution for Eq. (2)with q ( t = 0) = q and ˙ q ( t = 0) = v in case of moderate OscillatorRFIM Output: x Input: B ∼ x Output: M Input: F M ∼ M FIG. 2. The feedback mechanism of the RFIM coupled to aharmonic oscillator: The actual position x of the oscillatoraffects the magnetic field B (input) of the RFIM. Thereforethe spin configuration and also the magnetization M (out-put) depend on x . Then the magnetization again acts on theoscillator as an additional external force. damping ζ > q ( t ) = CM + 1 κ (cid:2) Ωζ sin Ωt − ( Ω −
1) cos Ωt (cid:3) + 12 κη e − ζt (cid:2) η cos ηt ( Ω − κ ( q − CM )) − sin ηt (2 ζ ( Ω + 1) − κ ( v + q ζ − ζCM )) (cid:3) , (3)where a and b are given by: κ = ( Ω − + 4 Ω ζ (4) η = (cid:112) − ζ . (5)We will use this solution later to avoid numerical inte-gration, when simulating the system (see III D) and alsoto analytically calculate bifurcation points for single-spindynamics (see IV A). B. Random field Ising model
For completeness and later reference we introduce herethe general RFIM with nearest neighbor couplings andits zero-temperature dynamics. It simplifies considerablyfor independent spins as detailed in section II C.The RFIM is used to determine the magnetic force F ( t ) that acts on the iron mass. The input and outputof the RFIM is the external magnetic field B and themagnetization M of the iron mass, respectively, whichin general is depending on the position of the mass. (cf.Fig. 2). The total magnetization results from a super-position of the magnetization of N discrete spins, whosestates are given by σ i ∈ {− , +1 } . Consequently, theRFIM is a discrete subsystem because the total magne-tization can only take a discrete set of values, which isalso known as quantization (see e.g. [40]). Moreover, thespin flips are assumed to occur instantaneously, whichmeans that also the change of the magnetization occursinstantaneously equivalent to an impact. The energy ofa specific configuration of the RFIM can be given by itsHamilton function H = − ˜ J (cid:88) (cid:104) ij (cid:105) σ i σ j − µ B (cid:88) i ( B + ˜ b i ) σ i , (6) where ˜ J is the coupling constant between two spins, µ B is the magnetic moment and ˜ b i is the local field of the i th spin. The (cid:104) ij (cid:105) indicates a sum over nearest neighborpairs, where each pair of spins is only counted once. Sincewe are interested in a coupling of the mechanical oscilla-tor with the RFIM, the external magnetic field B = B ( q )is assumed to be a function of the oscillator position q . Inthe following, we focus on the case of a linear dependence B ( q ) = β + β q, (7)which means that the oscillator displacements are lin-early coupled to the variations of the magnetic field. Sim-ilar to the continuous subsystem we will use a dimension-less Hamiltonian for the RFIM H = H µ B β = − J (cid:88) (cid:104) ij (cid:105) σ i σ j − (cid:88) i ( q + b i ) σ i , (8)where J = ˜ Jµ B β and b i = ˜ b i + β β are the dimensionlesscoupling constant and the dimensionless local disorderfield, respectively. The values b i , i = 1 , . . . , N are pa-rameters of a specific realization of the RFIM, which arechosen randomly and kept fixed during the time evolutionof the system (local quenched disorder). In particularwe choose the random fields to be Gaussian distributedand uncorrelated variables with b i b j = R δ ij and b i = 0.Here X denotes a quenched average of X , i. e. an averageover all disorder realizations of the system.In this paper we consider the RFIM at zero tempera-ture, which means that there are no stochastic spin flipsand the system dynamics is fully deterministic. The so-called single-spin-flip dynamics is used to describe theinternal dynamics of the RFIM subsystem [41, 42]. Inthis case the RFIM changes its spin configuration only ifa single spin flip would lower the energy of the subsys-tem, i.e., the Hamiltonian function H in Eq. (8) after thespin flip is smaller than the initial value before the spinflip. The energy difference ∆H i for a single flip of the i thspin can be given by ∆H i = 2 σ i J (cid:88) j ∈ nn( i ) σ j + q + b i , (9)where (cid:80) j nn( i ) indicates the sum of j over the near-est neighbors of the i th spin. Hence, the i th spin flipsif ∆H i <
0. Eq. (9) can be used to define so-called metastable states , which are spin configurations whereno single spin flip would lower the energy of the RFIMsubsystem, that is ∆H i ≥ i = 1 , . . . , N . Hence,these meta stable states full fill the so called metastabil-ity condition : σ i = sgn( F i ) (10)where F i = J (cid:80) j ∈ nn( i ) σ j + q + b i is the local field of the i th spin.Now, the dynamics of the RFIM subsystem can bedescribed as follows. The position q of the oscillator isupdated according to Eq. (2) until the energy differencefor any single spin flip is lower than zero ( ∆H i < ∆H i fora possible following spin flip is calculated for the newspin configuration. This is necessary because one spinflip may cause an avalanche of spin flips. If there isanother spin with ∆H i <
0, this spin is also reversedand the procedure is repeated until the system reachesthe next metastable state. If the RFIM has reached thenext metastable state and Eq. (10) is full filled for everyspin, the position q of the oscillator is updated againand the spin configuration does not change until themetastable state becomes again unstable. The presentupdate mechanism for the RFIM is known as sequentialupdate because the next metastable state is achieved bya sequence of single spin flips. It can be shown that dur-ing avalanches also other update mechanism, as for ex-ample parallel or synchronous update, lead to the samemetastable state. Moreover, also a different order of thesingle spin flips would not change the next metastablestate. This is due to the so-called no passing rule [43],which means that no spin flips more than once (eitherfrom down to up or vice versa) during the transition tothe next metastable state [44]. It might be worth to notethat the analysis of numerical algorithms for updatingthe internal states of the RFIM and its connection tograph theory is an enduring field of research [45–52].The output of the RFIM is the normalized dimension-less magnetization M defined by M = 1 N (cid:88) i σ i . (11)The connection between the dimensionless magnetization M and the original magnetization M of the iron mass isgiven by M = ρ m M , where ρ m is the magnetic dipoledensity. The magnetic force F M on the oscillator can bedetermined by F M = − ∂ x H , which leads to the dimen-sionless force F ( t ) = CM ( t ) , C = µ b β kNA . (12)At each time t the oscillator position q ( t ) determines thespin configuration of the RFIM, and therefore the mag-netization M = M ( t ). In comparison to the oscillatordynamics the update of the RFIM can be characterizedas adiabatic limit because at each time t the discrete sub-system is always in a metastable state. C. Independent spins
In this paper, we consider the case of independentspins, i.e., J = 0. In this case, there are no nearestneighbor interactions and no spin avalanches. Also thephenomena of first-order phase transition in dependency of the randomness R and of the dimension of the systemvanish. Since there are no nearest neighbor interactionsthe spatial arrangement of the spins is irrelevant. FromEq. (10) we find, that the condition for metastable statestakes the simple form σ i ( t ) = sgn( q ( t ) + b i ) ∀ i = 1 , . . . , N. (13)Eq. (13) can be understood as an definition of the spindynamics of the system. For a given position of the os-cillator at time t , each spin of the RFIM points in thedirection of its local field. Thus the magnetization canbe calculated by M ( q ( t )) = 1 N (cid:88) i sgn( q ( t ) + b i ) . (14)For our case of independent spins Eq. (14) directly deter-mines the magnetization M ( q ( t )) in dependence of q ( t ).Therefore the metastable state is always equivalent to theground state of the RFIM with the lowest possible energy.As a result for J = 0 the hysteresis feature vanishes andno memory develops in our system.Nevertheless, the case of independent spins is not triv-ial because there are still discrete changes of the force F ( t ) and the hybrid character of the system does notvanish. In fact, the limitation J = 0 gives us the pos-sibility to calculate exact results for the smooth systemin the limit of N → ∞ and to study the transition fromthe hybrid system for large but finite N and the smoothsystem with infinite N . III. METHOD
In this section we want to explain some details of thedynamics of the hybrid system and the calculation of Lya-punov exponents for such systems. In addition, we derivea smooth representation in the thermodynamic limit withinfinitely many spins ( N → ∞ ) and make some remarkson the implementation of the numerical methods. A. Dynamics of the piecewise-smooth system
If the spins are ordered according to their local disorderfields b i , it can be seen that the function M ( q ) is a stepfunction with N + 1 different levels of the magnetizationat which the force F acting on the oscillator is constant M i = 2 iN − , i = 0 , . . . , N, (15)implying, that M ( q ) is piecewise constant. Such systemsare called piecewise-smooth systems. The regions { S i } of constant F are separated by N boundaries at whichthe spin flips occur. Hence, one can argue that there isone ODE in each of these regions of the phase space andthe system state is completely determined by knowing t , S i − S i Φ i − x ∗ i − ,i Σ i Φ i FIG. 3. An illustration of our piecewise-smooth system. Thetwo smooth regions S i − and S i are separated by the bound-ary Σ i . The intersection point of the flow Φ i with the bound-ary is called x ∗ . q ( t ), and ˙ q ( t ). In other words, our system behaves likea piecewise-harmonic oscillator with same stiffness andsame damping ratio but with a forcing, which dependson the regions S i .Thus, we can write our system as a combination of aset of N + 1 ODEs˙ x ( t ) = F i ( x ( t ) , t ) = v ( t ) − ζv ( t ) − q ( t ) + cos φ ( t ) + CM i Ω , (16)with i = 0 , . . . , N . The state variable x = ( q, v, φ ) T is anelement of the smooth regions x ∈ S i , with (cid:83) i S i = D ⊂ R . The two-dimensional manifold Σ i = { x : H i ( x ) = 0 } with the indicator function H i ( x ) = q + b i separates twoneighboring regions S i − and S i . The intersection pointof the trajectory with the boundaries and the flow inthe region S i is denoted by x ∗ i − ,i and Φ i ( x ( t ) , t ), re-spectively (see Fig. 3). An exemplary trajectory of thesystem with N = 3 spins in the state space is illustratedin a projection in Fig. 4. We consider a point on thetrajectory with q < − b (black bullet ( • ) in left figure).This means that all spins are in the down-state. The sys-tem is evolved with magnetization M = −
1. After sometime we have − b < q < − b and the first spin has beenflipped to the up-state. In this region of the phase spacethe system further evolves with magnetization M = − .After passing the next boundary with − b < q < − b thenext spin will flip and the magnetization is M = . For q > − b all spins are in the up-state with M = 1 (rightfigure). When the oscillator moves in the other direction,the spins flips occur in reverse order. B. Thermodynamic limit N → ∞ For N spins there are N boundaries and at each bound-ary the magnetization increases by the value N . Forincreasing N , on one hand, the number of boundaries in-creases and, on the other hand, the changes of the mag-netization decrease. Thus, in the limit N → ∞ (ther-modynamic limit) the hybrid character of the piecewisesmooth system should vanish. In the following we de-rive a smooth function M ( q ) for the magnetization in dependence on the oscillator position q , representing thebehavior of the system in the thermodynamic limit.For q → −∞ each spin is in the down state and themagnetization is given by M = −
1. For increasing q themagnetization increases monotonically and in the limit q → + ∞ we have M = +1. The specific shape ofthe function M ( q ) is determined by the positions of theboundaries. Since the location of the boundaries is deter-mined by the local disorder b i , the probability density ofthe boundaries is a Gaussian distribution p ( q ) with zeromean and standard deviation R . The associated cumu-lative distribution can be written as P ( q ) = q (cid:90) −∞ p ( q (cid:48) ) d q (cid:48) = 12 (cid:18) (cid:18) qR √ (cid:19)(cid:19) . (17)It determines the ratio between the number of spins in thedown-state and the number of all spins in dependence ofthe oscillator position q . Therefore, in the limit N → ∞ we can substitute the ratio iN in Eq. (14) by P ( q ) andobtain a smooth function for the magnetization M ( q ( t )) = 2 P ( q ( t )) − (cid:18) q ( t ) √ R (cid:19) . (18)This means that in the limit N → ∞ our smooth dy-namical system consisting of a driven damped harmonicoscillator coupled to a RFIM with infinitely many inde-pendent spins can be described by˙ x ( t ) = v ( t ) − ζv ( t ) − q ( t ) + cos φ ( t ) + C erf (cid:16) q ( t ) R √ (cid:17) Ω , (19)with a smooth nonlinearity given by the error function.Thus, in Eq. (19) the feedback from the RFIM is char-acterized by an additional nonlinear external force thatdepends on the oscillator position q ( t ). In Sec. IV wecompare the dynamics of the piecewise smooth systemEq. (16) with a large but finite number of spins N andthe dynamics of the smooth system Eq. (19) with in-finitely many spins. C. Lyapunov exponents
Lyapunov exponents are defined as the average rateof divergence or convergence between a reference trajec-tory and a perturbed trajectory, where the perturbationsare infinitely small. For the smooth dynamical systemEq. (19) we use the standard method [53] for calculat-ing Lyapunov exponents. We use the same method inthe smooth regions S i of the piecewise smooth systemEq. (16). However, if the reference trajectory crossesa discontinuity boundary, the determination of the dy-namic behavior of the infinitesimal perturbations is notstraightforward because the perturbed trajectory mayhave crossed or may cross the boundary at an earlier − − − − V e l o c i t y v − − − b − − − b − b − − − b − b − b Position q − − . . M ag n e t i z a t i o n M FIG. 4. Example of a projection of a state space trajectory of the harmonic oscillator coupled to N = 3 spins. There are foursmooth regions S i with piecewise constant magnetization ( ) corresponding to four different spin configurations ( ↑↓↑ ). Ateach of the three boundaries a jump discontinuity appears in the acceleration ¨ q . The values of the local disorder are b = − . b = 1 . b = − . or a later time instant, respectively. In general, in theneighborhood of the discontinuities a careful treatmentof the determination of the perturbations is necessary be-cause the switching behavior of the perturbed trajectoryis typically different from the switching behavior of thereference solution. The compensation of such deviationscan be done by using the concept of the so-called Discon-tinuity Map (DM) [54–57]. In the following, we describeat first the basic concept of the DM for transversal inter-sections of the boundary, where the reference trajectorycrosses the boundary. Later we recall the concept of grac-ing intersections, where the reference trajectory hits theboundary tangentially. a. Transversal intersection – discontinuity map-ping During the calculation of the Lyapunov exponents weonly know the time instant t ∗ , where the reference tra-jectory reaches a discontinuity boundary and where weswitch between two different ODEs. At t ∗ the perturbedstate is in the neighborhood of the boundary but mayhave crossed it in the past or will cross it in the future.If the perturbed trajectory crossed (will cross) the bound-ary at an earlier (a later) time t ∗ + δt with δt < δt > t ∗ + δt by using knowl-edge of the state x ∗ , approximates the perturbed stateat the crossing time, applies the effects from the disconti-nuity crossing, and approximates the perturbed state attime t ∗ by evolving the perturbations before intersectingthe boundary to the to perturbations after the dynamicshas been switched. In other words, the DM immediatelyincorporates the effects of a discontinuity crossing even ifthe state is only in the neighborhood of a boundary andthe crossing appeared in the past or will appear in thefuture. For our system, the DM from a region S i to aregion S j can be given by the map Q ij : x → Q ij ( x ) = qC ( M i − M j ) δt + vφ , (20)where δt = q ∗ − qv ∗ . Since there is only a force jump atthe boundary the DM only changes the velocity of the state variable. Note that only changes with | i − j | = 1are possible and that the change in the magnetization is ∆M = M j − M i = ± δ x the Jacobian X = ∂ x Q ( x ) of the DM evaluated at the intersection point x ∗ can be used to calculate the perturbation δ x + ( t ∗ ) afterthe crossing by δ x + ( t ∗ ) = X δ x − ( t ∗ ) . (21)where δ x − ( t ∗ ) denotes the perturbation before the in-tersection with the boundary. The matrix X is called saltation matrix . For our system it has the form X = v ∗ C∆M . (22)Then, for a trajectory with only one crossing at time t ∗ the largest Lyapunov exponent would be defined as λ = lim t →∞ t ln | Y ( t, t ∗ ) XY ( t ∗ , δ x || δ x | , (23)where Y ( t, t (cid:48) ) is the fundamental solution of the varia-tional equation δ ˙ x ( t ) = D δ x ( t ) of the ODE (16) fromtime t (cid:48) to t and D denotes the Jacobian. The calculationof the Lyapunov exponent via Eq. (23) is illustrated inFig. 5 and can be explained as follows. We start withan initial perturbation δ x at time t = 0, and evolvethe perturbations up to the intersection time t ∗ , whichis calculated from the reference trajectory. At this pointwe have δ x − ( t ∗ ) = Y ( t ∗ , δ x . Then, the effect of theDM is captured by applying the saltation matrix X to δ x − ( t ∗ ). After the intersection we can use Y ( t, t ∗ ) again,to compute the perturbation up to time t . Note that theJacobian D and consequently the fundamental solution Y ( t, t (cid:48) ) is independent of the associated phase space re-gion S i because only a constant term, i.e. the magneti-zation M i , changes at the boundaries in the ODE (16). δ x (0) Y ( t ∗ , x Φ δ x − − b x ∗ δ x + Φ δ x ( t ) Y ( t, t ∗ ) X ( t ∗ ) v q FIG. 5. Illustration of the calculation of the Lyapunov expo-nent for piecewise-smooth systems. To evolve a small pertur-bation δ x in the smooth region, we can use the fundamentalsolution Y ( t , t ) of the variational equation of the ODE. Theeffect of the boundaries on δ x is described by the saltationmatrix X . b. Grazing intersection – Poincar´e-section and zerotime discontinuity mapping Note, that the saltation matrix has a singularity for v ∗ →
0, when the trajectory hits the boundary tangen-tial. This case is called grazing intersection . In generalgrazing occurs if the trajectory hits the boundary tan-gentially with the velocity normal to the boundary beingequal to zero. The point x ∗ , which belongs to the graz-ing intersection is called grazing point . Similar to thecase of transversal intersections, we have to make a cor-rection, when calculating Poincar´e maps for trajectoriesstarting near a orbit, which grazes the boundary. Thesecorrection arise, because some trajectories will not inter-sect the boundary, whereas others do cross. There aretwo common corrections to this, called Poincar´e SectionDiscontinuity Mapping (PDM) Q P and Zero Time Dis-continuity Mapping (ZDM) Q Z [58]. Both correctionsare constructed in the same manner as the transversalDM and they are describing the same grazing scenario.But whereas the PDM is defined with respect to a givenPoincar´e section, the ZDM is defined, such that zero timehas been elapsed between perturbation before and afterthe intersection with the boundary.It has been shown, that in case of the degree of smooth-ness of one, one has to add square-root terms of x in themapping of points near grazing [57]. For the artificialexample illustrated in Fig. 5 the grazing point is givenby x ∗ = ( − b, , φ ∗ ) T , hence the PDM takes the followingform: x → Q P ( x ) = (cid:40) x for H = q + b < , ν z √ q + b for H = q + b > , (24) ν z = 2 √ CΩ ( M − M )( b + cos φ ∗ + CM )( b + cos φ ∗ + CM ) , (25) and the corresponding ZDM can be written as: x → Q Z ( x ) = (cid:40) x for H min < , ν √ H min for H min > , (26) ν = 2 √ C ( M − M )( b + cos φ ∗ + CM )( b + cos φ ∗ + CM )( b + cos φ ∗ + CM ) , (27)where H min ( x ) is the minimum value of H ( x ) obtainedalong the flow Φ ( x , t ), where x is a arbitrary point nearthe grazing point x ∗ .For piecewise-smooth systems with a degree of smooth-ness of one it has been shown, that the dynamics – es-pecially the different scenarios, which can occur in bifur-cation diagrams – can be explained by piecewise-smoothdiscontinuous square-root maps like (25) [58]. The anal-ysis of piecewise-smooth square-root maps reveals, thatthose maps can describe various bifurcation scenarios in-cluding period-adding and robust chaos [59–61]. D. Numerics
In general for J (cid:54) = 0 the generation of trajectories ofthe hybrid system is done in the following way. The ex-ternal force of the initial metastable state of the RFIMis determined and used to solve the continuous subsys-tem until the next boundary is reached, i. e., until themetastable state of the RFIM becomes unstable. Then,the new metastable state of the RFIM is calculated byusing the single-spin-flip update, and the external forcecorresponding to the new metastable state of the RFIMis used to solve the continuous subsystem up to the nextboundary crossing.In our case of J = 0 the simulation of the RFIM be-comes obsolete, because of the absence of memory in thesystem. Hence, from the indicator function H i ( x ) = 0 = q + b i , we can determine the boundaries in the phase spacedirectly before starting the actual propagation of the tra-jectories (see III A). Also, since the continuous subsystemis linear in between the boundaries, it is possible to an-alytically calculate the trajectory of the oscillator for afixed external magnetic force (see Eq. (3)). Neverthelesswe are not able to calculate the time the oscillator needspropagating from one boundary to the next. In doing sowe would have to solve Eq. (3) for t and this can notbe done in a analytic way. Therefore for a low numberof spins, i. e. a “large” distance between two boundaries,we then use the analytical solution with a fixed stepsize ∆t to propagate the trajectory until a boundary givenby 0 = q + b i is crossed. Then, a root finder is used tocalculate the exact time t ∗ and state x ∗ at the intersec-tion point. The stepsize ∆t is chosen according to theresults in Fig. 6 such that no boundary crossings will beskipped.For an increasing number of spins in the system, thenumber of boundaries becomes large and the distances − − − − − . ∼ Number of spins N M e a n o f m i n . b o und a r y d i s t . Z m i n R = 2 . R = 1 . R = 1 . R = 0 . R = 0 . Z min on the number ofspins for different values of randomness R . The mean Z min was calculated for 500 different realizations of b i . between two boundaries decrease. Thus if the numberof spins is large, we use an adaptive scheme to solve thesystem, because the performance of the above-mentionednumerical solution with a root finder is low for large N . However in this case, there are still regions in thephase space, where the distance between two boundariesis large. This holds, for example, at very low and highvalues of q , where the probability for the occurrence ofa boundary is small. In these regions, we again use theanalytical solution of the continuous subsystem in com-bination with the root finder. On the other hand, if thedistance between two boundaries becomes “small”, wecalculate the next intersection point ( t ∗ , x ∗ ) directly byassuming a constant velocity of the oscillator betweentwo boundary crossings. This is similar to a linearizationof Eq. (3). By using this adaptive method, we are ableto generate bifurcation diagrams with a finite but veryhigh number of discontinuities (see Sec. IV).To get quantitative information on the size of “small”and “large” distances between two boundaries and howto chose the corresponding stepsize ∆t , such that noboundary crossings will be skipped, we calculated the ex-pectation value for the minimum distance between twoboundaries from the probability distribution of the localdisorder values b i in dependency on the number of spinsand the randomness R . Specifically, the local disordervalues b i can be written as a multivariate random vari-able b = ( b , b , . . . , b N ), where the b i are independentand identically distributed random variables with b i ∼N (0 , R ). Since we are interested in the distances be-tween two boundaries, we consider one sorted realizationof b , called s = ( s , s , . . . , s N ) with s ≤ s ≤ · · · ≤ s N .Hence, we can define the distance between consecutiveboundaries for one realization as z = ( z , z , . . . , z N − )and z i = s i +1 − s i , and denote z min the minimum of z . Because z min depends on the realization of b , thereis a related random variable Z min . The dependency of the expectation value Z min on the number of spins N for different values of R is presented in Fig. 6, where X indicates the expectation value for the random variable X over different realizations of b . It can be seen, thatfor an increasing number of spins, the mean minimal dis-tance decreases algebraically to zero with an exponent ofroughly N − . . IV. RESULTS
In this section we present the main results, startingwith the single-spin and many-spin dynamics. Next westudy the transition to the thermodynamic limit, i.e., thetransition from the piecewise-smooth system with an in-creasing number of spin-flip boundaries to the smoothsystem. Moreover, we present some results on the fractaldimensions of the chaotic attractors and take a look atthe behavior of the magnetization for an increasing num-ber of spins. All numerical simulations have been doneby fixing the normalized eigenfrequency as Ω = 1 . ζ = 0 . A. Single-spin dynamics
We start by investigating the dynamic behavior of thesystem with one spin N = 1. We have calculated thelargest Lyapunov exponent and bifurcation diagrams fordifferent initial values q and v . We have found mul-tistablility in a various parameter ranges. According toEq. (13) and (16) with N = 1, an exemplary basin ofattraction for the parameters b = 0 . C = 1 .
65 ispresented in Fig. 7, where chaotic solutions with a largestLyapunov exponent greater than zero and periodic solu-tions with vanishing Lyapunov exponent are indicated byblack and white boxes respectively.One can see the effects from the discontinuity at theposition of the spin flip q = − b = − . v = 0,where the saltation matrix has a singularity. Multistabil-ity can be observed for the asymmetric case with b (cid:54) = 0.In contrast, for b = 0 the spin flip position would beequivalent to the equilibrium position of the oscillatorand in this case no multistable behavior can be observed.A typical bifurcation diagram for the symmetric case( b = 0) is presented in Fig. 8. It shows the displacementof the oscillator q ( φ = 0) at the Poincar´e section φ = 0and the corresponding Lyapunov exponent of the asymp-totic solution. The bifurcation diagram is generated withthe fixed initial conditions q = − . v = 0 . C . The system shows the typicalscenarios, which are known for piecewise-smooth square-root maps [59–61]. This is in accordance with the ac-tual square-root dependency of the PDM and ZDM fromEq. (25) and Eq. (27) found for our system. The threecorresponding bifurcation scenarios, which appear dueto the discontinuity with degree of smoothness one, areoutlined by the three colored boxes in Fig. 8. The green − − − − − − V e l o c i t y v Position q FIG. 7. The system with one boundary at x = − . (cid:4) ) indicate chaos with a largestLyapunov exponent greater than zero while the white boxes( (cid:3) ) correspond to periodic behavior. The parameters are b =0 . C = 1 .
65 in Eq. (13) and (16) for N = 1. box demonstrates an overlapping period-adding cascade,which in the case of decreasing values of C starts at C ∗ = 10. This is in agreement with the prediction wecan make by using the solution of the harmonic oscilla-tor with constant magnetization from Eq. (3). Whenstarting at the left side of the boundary q ∗ = 0 we cancalculate for large t the maximum q -values of the peri-odic orbit of the system. By assuming, that this orbittouches the boundary if q max = q ∗ , we find a formula for C ∗ : C ∗ = 1 M (cid:18) q ∗ − √ κ (cid:19) . (28)For M = − q ∗ = 0 and ζ = 0 .
05 we find C ∗ = 10. Theblue box illustrates period-adding with chaotic segmentsin between and the red box shows an immediate jumpfrom chaos again to a periodic solution.For a non-zero disorder parameter b (cid:54) = 0, in gen-eral, the qualitative behavior of the bifurcations is similarto the bifurcations in the symmetric case. However, inthe asymmetric case b (cid:54) = 0 the location of the periodicwindows and the chaotic regions can slightly change de-pending on the specific initial condition. Moreover, it isworth to emphasize that the system without discontinu-ity ( N = 0) does not show any chaos because it reducesto the dynamics of a damped harmonic oscillator withperiodic excitation. This means that the origin of chaosin the system with one spin is the piecewise constantmagnetization that jumps at the spin flip position. B. Many-spin dynamics
For systems with a few number of spins the dynamicbehavior and the bifurcation diagrams look similar to − . .
05 0 2 4 6 8 10 λ m a x C − − P o i n c a r ´ e s ec t i o n q ( φ = ) FIG. 8. Bifurcation diagram for the totally symmetric sys-tem with N = 1 and b = 0. The system shows the bifur-cation scenarios expected from piecewise-smooth square-rootmaps, illustrated by the different colored boxes (from left toright): immediate jump to robust chaos ( ), period-addingwith chaos ( ) and overlapping period-adding cascade ( ). the one spin case and only the number of discontinu-ities may be different. However, if the number of spinsis much higher than one, the characteristic properties ofthe system change. A bifurcation diagram and the corre-sponding Lyapunov exponent for a high number of spins N = 20 000 and a fixed realization of the b i is presentedin Fig. 9. In this case and for all following numericalcalculations the degree of randomness of the disorder ischosen as R = 1 .
7. First we report, that the largest Lya-punov exponent λ max for the many-spin system evaluatedfor chaotic regions is roughly two times larger than λ max of the system with only one spin. This is not obvious,because for an increasing number of spins the height ofthe magnetization jumps at the boundaries goes to zeroand the saltation matrix converges to the identity. On theother hand, the typical bifurcation scenarios from grazing(period-adding, immediate jump to chaos) vanish, whichis also not obvious because the number of boundariesand discontinuities is much higher than for the one spinsystem. This indicates that the chaotic behavior onlyarises due to transversal intersections with the bound-aries. The dynamic properties of the system with manyspins are, in general, very similar to the dynamic prop-erties of the smooth system in the thermodynamic limit( N → ∞ ). This can be seen, for example, by comparingthe bifurcation diagram and the maximum Lyapunov ex-ponent of the piecewise-smooth system with N = 20 000spins and their counterparts calculated from the smoothsystem, which are presented by the red curves in Fig. 9.Overall both the black curves for the piecewise-smoothsystem and the red curves for the continuous system lookvery similar. But a more detailed view (see Fig. 10) ofthe bifurcations for C ∈ [2 ,
3] shows, that the two dia-grams are slightly different. This is due to the fact, that0 − . . . λ m a x C − − P o i n c a r ´ e s ec t i o n q ( φ = ) (cid:4) N = 20 000 (cid:4) N → ∞ FIG. 9. Comparison of the bifurcation diagrams for thepiecewise-smooth system with N = 20 000 ( (cid:4) ) spins and thesystem in its thermodynamic limit ( (cid:4) ). It can be seen, thatthe typical grazing scenarios (immediate jump to chaos andperiod-adding cascades) vanish, whereas the main behavior ispretty similar. − . . . . . . . C − (cid:4) N = 20 000 (cid:4) N → ∞ λ m a x P o i n c a r ´ e s ec t i o n q ( φ = ) FIG. 10. Zoom of the bifurcation diagram from Fig. 9. Itcan be seen, that besides the general similarities between thepiecewise-smooth system and the system in its thermody-namic limit there are some differences in the actual behaviorof the bifurcations. Mainly this is because of the dependenceof the dynamical properties of the system for N = 20 000 onthe actual realization of the local disorder { b i } . even for N = 20 000 the behavior of the system dependsnoticeable on the actual realization of the disorder { b i } .The same observation can be made in the comparison ofthe chaotic attractors of the smooth and the piecewise-smooth system, which are presented in Fig. 11. Thereare nearly no differences in the macroscopic structure ofthe attractor and only small deviations can be seen atfiner scales. C. Transition to the thermodynamic limit
In the piecewise-smooth system with a finite numberof spins the origin of chaos lies in the discontinuity cross-ings, whereas in the smooth system with an infinite num-ber of spins chaos comes from the nonlinearity in the ad-ditional magnetic force. Nevertheless, for an increasingnumber of spins the dynamics of the piecewise-smoothsystem converges on macroscopic scales to the dynamicsof the system in the thermodynamic limit. In the follow-ing this transition is studied in more detail.On the one hand, for increasing N the number of dis-continuities increases, but on the other hand simultane-ously the influence of the discontinuities goes to zero,because the jump in the magnetization ∆M at each dis-continuity vanishes ( ∆M →
0) and the saltation matrix X converges to the identity ( X → I ) for N → ∞ . To getan idea of the interplay between the increasing number ofboundaries and the decreasing influence of an individualspin flip, we consider a small segment of the attractorwith length q g , which is divided by n boundaries. For avery large number of spins, we can assume that the lo-cation of the boundaries is homogeneously distributed inthe small segment with length q g , which means that thedistance between two boundaries can be approximatelygiven by ∆q = q g /n . Note that for a Gaussian distributedlocal disorder fields of the RFIM the average distance ∆x still varies with the location of the attractor segment inphase space. In addition, we assume the velocity at theintersection to be 0 < v < ∞ for ∆M >
0, which is avalid assumption, because a positive velocity leads to anincreasing q and a positive ∆M . The short-time Lya-punov exponents can be determined by λ i = 1 ∆t ln | µ i ( XY ) | , (29)where ∆t = ∆q/v , µ i ( A ) denotes eigenvalues of a matrix A , and the matrices X and Y are the saltation matrixand the fundamental solution of the harmonic oscillatorfor a time step ∆t , respectively. They are given by X = (cid:18) v C∆M (cid:19) , Y = exp (cid:20)(cid:18) − − ζ (cid:19) ∆t (cid:21) . (30)For very large N , ∆t is very small and the matrix ex-ponential in Y can be approximated via a linear Taylorapproximation. In this case, the Lyapunov exponents canbe determined by λ / = Re (cid:34) − ζ ± (cid:115) C ∆M∆q + ζ − (cid:35) . (31)For C < ∆M∆q the short-time Lyapunov exponent is pos-itive, which means that, in general, chaotic behavior ispossible. Note that for an increasing number of spinsthe jump of the magnetization ∆M = 2 /N vanishesbut simultaneously the average distance ∆x between twojumps vanishes. The ratio ∆M∆q converges to a positive1 − − − − − − − V e l o c i t y v Position q − − − − − . . M ag n e t i z a t i o n M N → ∞ N = 20 000(a) (b) FIG. 11. Comparison of the chaotic attractor ( ), Poincar´e section ( ) and magnetization ( ). (a): The system in itsthermodynamic limit N → ∞ . (b): The piecewise-smooth system with N = 20 000 spins and one specific disorder realization { b i } . Both systems are evaluated for fixed coupling constant C = 3 . R = 1 . constant specifying the average density of discontinuitiesin the given attractor segment. A high number of jumpsand/or a high coupling constant C increases the short-time Lyapunov exponent, and therefore, the probabilityto observe chaos.For the smooth system in its thermodynamic limit N → ∞ we can find a similar condition by linearizingthe nonlinear system Eq. (19) around a location q ∗ on theattractor. In this case, the corresponding short-time Lya-punov exponents at q ∗ can be calculated from Eq. (29)by substituting the matrix product XY with the matrix Y ∞ = exp (cid:34)(cid:32) − C ∂M ( q ) ∂q (cid:12)(cid:12)(cid:12) q ∗ − ζ (cid:33) ∆t (cid:35) , (32)leading to the two Lyapunov exponents λ / = Re (cid:34) − ζ ± (cid:115) C ∂M ( q ) ∂q (cid:12)(cid:12)(cid:12)(cid:12) q ∗ + ζ − (cid:35) . (33)By comparing Eq. (31) and Eq. (33) it becomes clearthat the discrete magnetization jumps in one attractorsegment of the piecewise-smooth system translates intoa continuous increase of the magnetization in the smoothsystem and the short-time Lyapunov exponent dependson the slope of the magnetization in this segment. Byinserting the explicit expression for the derivative of themagnetization derived from the distribution of the localdisorder of the spins, we obtain the condition1 C < (cid:114) πR e − q ∗ R (34)for a positive short-time Lyapunov exponent at q ∗ .The behavior of the condition in Eq. (34) is illustratedin Fig. 12. For C = 0, there are two Lyapunov exponents λ / = − ζ and no chaos is possible. For increasing C atsome point the term under the square root in Eq. (33)becomes positive and by further increasing C the domi-nant exponent becomes positive. Thus, increasing C in-creases the probability for observing chaos, which is clearbecause a higher coupling constant C leads to a higherweighting of the nonlinear magnetic force. At q ∗ = 0,corresponding to the position with the maximal densityof boundaries, the dominant short-time Lyapunov expo-nent has its maximum, and the exponent decreases for anincreasing | q ∗ | . This is clear because the maximum slopeof the nonlinearity in Eq. (19), and therefore the largestinfluence of the magnetic force, can be found at q ∗ = 0,whereas for q ∗ → ±∞ the slope goes to zero and the de-pendence of the magnetization on the oscillator positionvanishes. The dependence of the short-time Lyapunov . . . . . . q πR e − q ∗
22 1 R /R / C q ∗ = 0 q ∗ = ± q ∗ = ± q ∗ . R of the local disorder canbe explained as follows. For a small R (large 1 /R inFig. 12) most of the spin flips occur around q ∗ = 0. Thereis a large change of the magnetization around the equi-librium position but only slight changes at other points,where the density for spin flips is much lower. As a con-sequence the short-time Lyapunov exponent is likely tobe positive near q ∗ = 0 and becomes smaller and nega-tive for increasing | q ∗ | . In contrast, for a large variance R the changes of the density for observing spin flips islow, and similarly the variations of the magnetization forvarying oscillator position are low. As a consequence,a positive short-time Lyapunov exponent and probablychaos can be found only for very high C but then in abroad region around q ∗ . D. Fractal dimensions of the chaotic attractor
Since the system can still produce chaos, even for aninfinite number of boundaries, it is natural to ask for thedynamic properties of a typical chaotic attractor, whichis shown, for example, in Fig. 11. Hence we are interestedin the behavior of the box counting and the Kaplan-Yorkedimension D BC and D KY [62, 63]. Thus we calculated themean values of both dimensions D BC and D KY of the at-tractor for a varying number of spins N by using 500 dif-ferent realizations of the local disorder { b i } at each valueof N . The results are shown in Fig. 13, where the cou-pling strength of the magnetization is chosen as C = 3 . D ∞ BC ≈ .
56 and D ∞ KY ≈ .
54 (dashed lines in Fig. 13). The limit valuesof the mean of both dimensions for the piecewise-smoothsystem are D ∗ BC ≈ .
52 and D ∗ KY ≈ .
49, calculated byusing the mean of the last five fractal dimensions valuesfrom N = 17 500 to N = 19 500. We also plotted D ∗ − D in dependency on N , which is shown in the inset in Fig.13. We find, that D BC as well as D KY converges exponen-tially to their limit values D ∗ BC and D ∗ KY . This supportsthe proposition, that the piecewise-smooth system with avery large number of spins behaves like a harmonic oscil-lator with an additional nonlinear smooth external forceand the piecewise-smooth character vanishes. Note, thatthe Kaplan-Yorke dimension should be lower than thebox counting dimension according to the theory of thedimensions of chaotic attractors [64], which is also fullyreflected by our simulations. Nevertheless the questionremains whether the variance of the fractal dimensionof the chaotic attractor vanishes for a large number ofspins. To answer this question, we take a look at thecoefficient of variation, also often called Self-AveragingParameter (SAP) [65], of the fractal dimension D of the . . N → ∞ (a) Number of spins N M e a n o ff r a c . d i m . D BCKY − − − (b) N D ∗ − D FIG. 13. (a): For an increasing number of spins the boxcounting (red crosses) and the Kaplan-Yorke dimension (bluesquares) converge to the corresponding values of the smoothsystem (dashed lines) in the thermodynamic limit ( N → ∞ ).Here the coupling strength equals C = 3 . N , which can be seen in a semi-log plot. Hence there is an exponential convergence to thelimit values. attractor, which is given bySAP[ D ] = D − D D . (35)Here, as before the bar X denotes an average of X overdifferent realization of the quenched local disorder. TheSAP of the dimension was calculated for 500 disorder re-alizations { b i } and a varying number of spins N and ispresented in a semi-log plot in Fig. 14. We find, thatthe SAP approaches to zero also roughly exponentiallyfor N → ∞ , hence, the system shows self-averaging withrespect to the box counting and the Kaplan-Yorke di-mension of the attractors. This means, that for a smallnumber of spins the fractal dimension strongly dependson the realization of the local disorder, whereas for alarge number of spins, this dependency vanishes. Thus,fractal dimensions are self-averaging quantities and canbe calculated from one typical disorder realization of alarge system. E. Magnetization
Besides the investigation of the dynamic properties ofthe system, the behavior of the magnetization of theRFIM shows some interesting behavior. We numericallycalculated the variance of the magnetizationVAR[ M ] = (cid:18) N (cid:88) i (cid:104) σ i (cid:105) (cid:19) − M (36)over 500 disorder realizations for two typical chaotic at-tractors with C = 2 . C = 3 .
5. Here, (cid:104) σ i (cid:105) denotes3 − − − Number of spins N S e l f a v e r ag i n g p a r a m e t e r S A P [ D ] BCKYFIG. 14. The self-averaging parameter of the box countingand Kaplan-Yorke dimension in dependence on the number ofspins appears to decrease exponentially to zero. Therefore thesystem shows self-averaging with respect to these dynamicalproperties for C = 3 . the time-average of the configuration of the i th spin and M denotes the average of the magnetization of the sys-tem over the disorder realizations. In our case due tothe symmetry in the distribution of the disorder with re-spect to the oscillator equilibrium, we have M = 0. Theresulting variance is presented in Fig. 15. We foundthat for the attractor at C = 3 . N − , this is fully in accordance to the expected behav-ior of the variance within the central limit theorem. Incontrast, the magnetization does not show self-averagingfor the attractor at C = 2 . C = 2 .
9, in general, many different attractors showup and it depends on the specific disorder realization,which asymptotic state is reached by the system. Foran increasing number of spins the system mainly endsup in one of two symmetric attractors, which are illus-trated in Fig. 16 for N = 20 000. The time-average ofthe magnetization for the blue attractor ( ) is greaterthan zero, while the time-average of the red attractor( ) is smaller than zero. The distribution of the mag-netization is roughly symmetric and has two maxima atthe positive and negative magnetization correspondingto the blue and red attractors (see Fig. 16). As a con-sequence, the variance of the magnetization does not goto zero even for a large number of spins and, in general,depends on the actual dynamics of the system. − − − − − − − . ∼ Number of spins N V a r i a n c e o f m ag . VA R [ M ] i.i.d. C = 2 . C = 3 . C = 3 . N − forincreasing number of spins similar to iid input of the RFIM.In contrast, for C = 2 . − − − q V e l o c i t y v − . . M P D F N = 5 000 N = 20 000FIG. 16. For C = 2 . N = 20 000 there exist two differenttypical attractors of the system illustrated by the red ( ) andblue ( ) dots. This explains the non self-averaging behaviorof the magnetization. (b): The non self-averaging behaviorof the magnetization is also reflected in the correspondinghistogram by the two main bars at M = ± . V. CONCLUSION
Motivated by the phenomenon of complex hysteresisin many dynamical systems, we studied the exemplarysystem of a harmonic oscillator coupled to a simplifiedRFIM, where the input and output of the RFIM is the os-cillator position and the magnetic force from the RFIM,respectively. We focused on the piecewise-smooth char-acter of the system and neglected spin-spin interactionsin the RFIM. In this case, each spin flips at a fixed oscil-lator position, which is determined by the local disorderparameter of the spins. These positions correspond toparallel boundaries in the phase space. At the bound-aries the magnetic force jumps, whereas between theboundaries the force remains constant and the systemis smooth.The dynamics of the system with only a small numberof spins is dominated by different grazing bifurcation sce-narios, which are typical for piecewise-smooth systems.Chaotic solutions and multistability can be found alreadyfor the oscillator coupled to only one spin. For a large4number of spins, the grazing bifurcation scenarios vanishand the dynamic behavior of the piecewise-smooth sys-tem is very similar to the dynamic behavior of the smoothsystem in the thermodynamic limit with infinitely manyspins. This is not obvious because the number of dis-continuities increases. However, on the other hand thechanges of the magnetization per spin flip decrease. As aresult, the system becomes smoother and in the thermo-dynamic limit the system can be described by a harmonicoscillator with a smooth nonlinear magnetic force. Thesmooth system is also able to show chaos. The typicalbox counting and the typical Kaplan-Yorke dimensionof the chaotic attractors of the piecewise-smooth systemconverge to the corresponding dimensions of the smoothsystem in its thermodynamic limit. The variance of theattractor dimension vanishes for an increasing number ofspins. This does not hold for the magnetization because there is a bi-stability between two symmetric attractorswith a positive or negative average magnetization.In future work we will focus on the case which includesspin-spin interactions. In this case hysteresis is possiblein the RFIM, that is, the internal state of the RFIMis not necessarily determined only by the instantaneousoscillator position but also by past values of the input.This means, that the system can still be treated as apiecewise-smooth dynamical system but the boundariesin phase space, which are associated with the discontinu-ities due to a spin flip, are no longer fixed but become ahistory-dependent dynamic quantity.
VI. ACKNOWLEDGEMENTS
We would like to thank Sven Schubert for helpful dis-cussions and valuable suggestions. [1] A. Brinkman, M. Huijben, M. van Zalk, J. Huijben,U. Zeitler, J. C. Maan, W. G. van der Wiel, G. Rijnders,D. H. A. Blank, and H. Hilgenkamp, Nature Materials , 493 EP (2007).[2] C. Song, J. Brandon, and C. A. Featherston, J. Mech.Eng. Sci. , 673 (2001).[3] D. J. Late, B. Liu, H. S. S. R. Matte, V. P. Dravid, andC. N. R. Rao, ACS Nano , 5635 (2012).[4] I. Urbanaviciute, T. D. Cornelissen, X. Meng, R. P. Si-jbesma, and M. Kemerink, Nature Communications ,4409 (2018).[5] S. Horike, S.and Shimomura and S. Kitagawa, NatureChemistry , 695 EP (2009).[6] S. Eckel, J. G. Lee, F. Jendrzejewski, N. Murray, C. W.Clark, C. J. Lobb, W. D. Phillips, M. Edwards, andG. K. Campbell, Nature , 200 EP (2014).[7] P. Crespo, R. Litr´an, T. C. Rojas, M. Multigner, J. M.de la Fuente, J. C. S´anchez-L´opez, M. A. Garc´ıa, A. Her-nando, S. Penad´es, and A. Fern´andez, Phys. Rev. Lett. , 087204 (2004).[8] G. Bertotti and I. D. Mayergoyz, The science of hystere-sis , Vol. 1-3 (Academic Press, New York, 2006).[9] F. Preisach, Zeitschrift fr Physik , 277 (1935).[10] P. Weiss, J. Phys. Theor. Appl. , 661 (1907).[11] R. Peierls, Math. Proc. Camb. Phil. Soc. , 477 (1936).[12] Y. Imry and S. Ma, Phys. Rev. Lett. , 1399 (1975).[13] J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl,B. W. Roberts, and J. D. Shore, Phys. Rev. Lett. ,3347 (1993).[14] T. G. Sorop, C. Untiedt, F. Luis, M. Kroll, M. Rasa, andL. J. de Jongh, Phys. Rev. B , 014402 (2003).[15] J. Ort´ın, J. Appl. Phys. , 1454 (1992).[16] M. P. Lilly, A. H. Wootters, and R. B. Hallock, Phys.Rev. Lett. , 4222 (1996).[17] O. Perkovic, K. A. Dahmen, and J. P. Sethna, Phys.Rev. Lett. , 4528 (1995).[18] T. Nattermann, in Spin Glasses and Random Fields , Se-ries on Directions in Condensed Matter Physics, Vol. 12,edited by A. P. Young (World Scientific, New Jersey,1997) pp. 277–298. [19] P. Shukla, Phys. Rev. E , 4725 (2000).[20] J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature , 242 (2001).[21] I. D. Mayergoyz and C. E. Korman, J. Appl. Phys. ,5478 (1994).[22] P. Rugkwamsook, C. E. Korman, G. Bertotti, andM. Pasquale, J. Appl. Phys. , 4361 (1999).[23] M. Dimian and I. D. Mayergoyz, Phys. Rev. E ,046124/1 (2004).[24] A. Adedoyin, M. Dimian, and P. Andrei, IEEE Trans.Magn. , 3934 (2009).[25] G. Radons, Phys. Rev. E , 061133 (2008).[26] G. Radons, Phys. Rev. E , 061134 (2008).[27] G. Radons, Phys. Rev. Lett. , 240602 (2008).[28] S. Schubert and G. Radons, Phys. Rev. E , 022117(2017).[29] H. Lamba, M. Grinfeld, S. McKee, and R. Simpson,IEEE Trans. Magn. , 2495 (1997).[30] A. Rezaei-Zare, M. Sanaye-Pasand, H. Mohseni,S. Farhangi, and R. Iravani, IEEE Trans. Power Delivery , 919 (2007).[31] G. Radons and A. Zienert, Eur. Phys. J. Special Topic , 1675 (2013).[32] P. A. Cook, Systems & Control Letters , 223 (1985).[33] K. H. Johansson, A. Rantzer, and K. H. ˚Astrm, Auto-matica , 539 (1999).[34] J. M. Goncalves, A. Megretski, and M. A. Dahleh, IEEETransactions on Automatic Control , 550 (2001).[35] M. Wiercigroch, R. D. Neilson, and M. A. Player,Physics Letters A , 91 (1999).[36] S. Theodossiades and S. Natsiavas, Journal of Sound andVibration , 287 (2000).[37] K. Popp and P. Shelter, Philosophical Transactions:Physical Sciences and Engineering , 89 (1990).[38] U. Galvanetto, Journal of Sound and Vibration , 653(2001).[39] M. di Bernardo, C. Budd, A. R. Champneys, andP. Kowalczyk, Piecewise-smooth Dynamical Systems:Theory and Applications , Applied Mathematical Sci-ences, Vol. 163 (Springer London, 2008). [40] M. S. Santina and A. R. Stubberud, in Handbookof Networked and Embedded Control Systems , editedby D. Hristu-Varsakelis and W. S. Levine (Birkhuser,Boston, MA, 2005) pp. 45–69.[41] E. Vives, M. L. Rosinberg, and G. Tarjus, Phys. Rev. B , 134424 (2005).[42] F. Salvat-Pujol, E. Vives, and M. L. Rosinberg, Phys.Rev. E , 061116 (2009).[43] A. A. Middleton, Phys. Rev. Lett. , 670 (1992).[44] D. Dhar, P. Shukla, and J. P. Sethna, Journal of PhysicsA: Mathematical and General , 5259 (1997).[45] A. V. Goldberg and R. E. Tarjan, J. ACM , 921 (1988).[46] A. K. Hartmann and K. D. Usadel, Physica A , 141(1995).[47] A. A. Middleton, Phys. Rev. Lett. , 017202 (2001).[48] I. Dukovski and J. Machta, Phys. Rev. B , 014413(2003).[49] P. E. Theodorakis and N. G. Fytas, Condensed MatterPhysics (2014).[50] U. Wolff, Phys. Rev. Lett. , 361 (1989).[51] A. K. Hartmann, Physica A: Statistical Mechanics andits Applications , 480 (1996).[52] M. C. Kuntz, O. Perkovic, K. A. Dahmen, B. W. Roberts,and J. P. Sethna, arXiv e-prints , cond-mat/9809122 (1998).[53] G. Bennetin, L. Galgani, A. Giorgilli, and J. Strelcyn,Meccanica , 9 (1980).[54] A. B. Nordmark, Journal of Sound and Vibration ,279 (1991).[55] P. C. M¨uller, Chaos, Solitons & Fractals , 1671 (1995).[56] H. Dankowicz and A. B. Nordmark, Physica D: NonlinearPhenomena , 280 (2000).[57] M. di Bernardo, C. J. Budd, and A. R. Champneys,Phys. Rev. Lett. , 2553 (2001).[58] M. di Bernardo, C. J. Budd, and A. R. Champneys,Physica D: Nonlinear Phenomena , 222 (2001).[59] C. Budd and F. Dux, Nonlinearity , 1191 (1994).[60] W. Chin, E. Ott, H. E. Nusse, and C. Grebogi, Phys.Rev. E , 4427 (1994).[61] S. Foale and S. R. Bishop, Nonlinear Dynamics , 285(1994).[62] H. G. E. Hentschel and I. Procaccia, Physica D: Nonlin-ear Phenomena , 435 (1983).[63] P. Grassberger and I. Procaccia, Physica D: NonlinearPhenomena , 189 (1983).[64] J. D. Farmer, E. Ott, and J. A. Yorke, Physica D: Non-linear Phenomena , 153 (1983).[65] A. Aharony and A. B. Harris, Phys. Rev. Lett.77