Can Drag Force Suppress Fermi Acceleration in a Bouncer Model?

Abstract

Some dynamical properties for a bouncer model—a classical particle of mass 𝑚 falling in the presence of a constant gravitational field 𝑔 and hitting elastically a periodically moving wall—in the presence of drag force that is assumed to be proportional to the particle's velocity are studied. The dynamics of the model is described in terms of a two-dimensional nonlinear mapping obtained via solution of the second Newton's law of motion. We characterize the behavior of the average velocity of the particle as function of the control parameters as well as the time. Our results show that the average velocity starts growing at first and then bends towards a regime of constant value, thus confirming that the introduction of drag force is a sufficient condition to suppress Fermi acceleration in the model.

1. Introduction

The bouncer model consists basically of a classical particle of mass 𝑚 that moves under the presence of a gravitational field 𝑔 and suffers elastic collisions with an infinitely heavy and periodically moving wall. It is well known in literature [1–3] that, depending on both the combination of control parameters and initial conditions, the particle can accumulate, along the orbit, large values of energy thus exhibiting a phenomenon called as Fermi acceleration [4].

The bouncer model has been studied along last decades under different versions [5, 6] and considering several approaches including theoretical and experimental. For the theoretical approach, results like fixed point stability [7, 8], manifold crossings [9], period doubling cascades, and fully chaotic behavior were discussed [10] both as function of the initial conditions and control parameters. For the case of experimental approach, many results have been obtained considering the presence of inelastic collisions [11–13] including the observance of locking phenomenon.

There exist many situations where one can observe Fermi acceleration in one-dimensional models. One of them is assuming that the perturbation is no longer periodic but rather stochastic [14]. This perturbation leads to a growth in the average velocity according to square root of the number of collisions of the particle with the randomly moving wall. It is more often to be studied, in such models, conditions where Fermi acceleration (FA) is feasible. Our approach in this paper is rather different. Instead of looking at conditions to produce FA, we will consider a condition to suppress such a phenomenon. Although many results are known in literature concerning the bouncer model, our approach is original in the sense we look at conditions to break down the unlimited energy growth.

In this paper, we consider the bouncer model with the particle experiencing a drag force that is proportional to its velocity. As usual in literature, we have described the model in terms of a two-dimensional nonlinear mapping for the variables velocity of the particle and time at the instant of the collisions with the moving wall. The mapping was obtained by analytical solution of Newton's law of motion. We show that the dissipation suppresses the phenomenon of Fermi acceleration.

The organization of the paper is as follows. In Section 2 we describe all the details used for the mapping construction. Section 3 presents and discusses our numerical results while final remarks and conclusions are drawn in Section 4.

2. The Model and the Mapping

In this section we describe with full details the model we are considering and the procedure used to obtain the nonlinear mapping. The system consists of a classical particle of mass 𝑚 which experiences a drag force that is proportional and contrary to the particle's velocity. The particle is suffering the action of a constant gravitational field 𝑔, and it collides elastically a time varying wall whose equation of motion is given by 𝑦_𝑤(𝑡)=𝜖cos(𝜔𝑡). Here 𝜖 is the amplitude of oscillation and 𝜔 denotes the angular frequency of the wall. As it is so common in literature [15–18] the dynamics of the model is described in terms of a two-dimensional mapping 𝑇(𝑣_𝑛,𝑡_𝑛)=(𝑣_𝑛+1,𝑡_𝑛+1) where 𝑣_𝑛 and 𝑡_𝑛 denote the posthit velocity of the particle and the instant of the 𝑛th collision with the moving wall, respectively.

To obtain the mapping that fully describes the dynamics of the system, we must first write down Newton's law of motion, which is given by 𝑚𝑑𝑣/𝑑𝑡=𝑚𝑔−𝜂^′𝑣 where positive 𝜂^′ is the coefficient of viscosity. Defining 𝜂=𝜂^′/𝑚 we can rewrite the differential equation in scalar form as 𝑑𝑣/𝑑𝑡=−𝑔−𝜂𝑣 where the ascendent motion has positive 𝑣, and 𝑔 is the modulus of gravitational acceleration. Integrating with respect to the time this equation assuming that at the time 𝑡=𝑡_𝑛 which is the time that the particle hits the moving wall, the particle is at the position 𝑦_𝑝(𝑡_𝑛)=𝜖cos(𝜔𝑡_𝑛) with initial velocity 𝑣(𝑡_𝑛)=𝑣_𝑛>0, and then we find that 𝑣(𝑡)=[𝑣_𝑛+𝑔𝜂]𝑒^{−𝜂(𝑡−𝑡𝑛)}−𝑔𝜂,𝑡≥𝑡_𝑛. (2.1)

Since the expression of the velocity of the particle is known, we must integrate with respect to the time again the relation 𝑑𝑦(𝑡)/𝑑𝑡=𝑣(𝑡). Then we have that the position of the particle is given by 𝑦_𝑝(𝑡)=𝜖cos(𝜔𝑡_𝑛1)+𝜂(𝑣_𝑛+𝑔𝜂)(1−𝑒^{−𝜂(𝑡−𝑡𝑛)})−𝑔𝜂(𝑡−𝑡_𝑛),𝑡≥𝑡_𝑛. (2.2)

To find the instant of the collisions of the particle with the moving wall, we must take into account two different kinds of collisions that can occur, namely, (1) successive collisions—those collisions that the particle experiences coming from inside the collision zone 𝑦∈[−𝜖,𝜖]; (2) indirect collisions—those collisions where the particle comes from outside the collision zone.

Let us discuss the equations of the mapping for case (1) first. Matching the condition 𝑦_𝑝(𝑡_𝑐)=𝑦_𝑤(𝑡_𝑐) we obtain the instant of the collision (𝑡_𝑐=𝑡_𝑛+1) of the particle with the moving wall. Such condition gives rise to the following transcendental equation: 0=𝐺(𝑡_𝑐)=𝜖cos[𝜔(𝑡_𝑛+Δ𝑡_𝑐)]−𝜖cos(𝜔𝑡_𝑛𝑔)+𝜂Δ𝑡_𝑐−1𝜂(𝑣_𝑛+𝑔𝜂)(1−𝑒^{−𝜂Δ𝑡𝑐}), (2.3) where the interval of collision (Δ𝑡_𝑐=𝑡_𝑐−𝑡_𝑛)∈(0,2𝜋/𝜔], while exist, is the smallest solution of the equation 𝐺(𝑡_𝑐)=0 along Δ𝑡_𝑐∈(0,2𝜋/𝜔]. At the instant of the collision, we assume a noninertial system in which the moving wall is at rest, then the exchange of momentum can be obtained. Finally we return to the original inertial system, leading to a mapping of the type 𝑇_𝑀∶𝑣_𝑛+1=−[𝑣_𝑛𝑒^{−𝜂Δ𝑡𝑐}+𝑔𝜂(𝑒^{−𝜂Δ𝑡𝑐}−1)]−2𝜖𝜔sin(𝜔𝑡_𝑛+1𝑡),_𝑛+1=𝑡_𝑛+Δ𝑡_𝑐. (2.4)

We now discuss case (2), that is, indirect collisions. In the approach we have considered there are three steps to obtain the time the particle spent between collisions in the bouncer model. The first one is the elapsed time from the last collision up to the maximum height where the particle acquires null velocity. The second one is the time the particle spent from the maximum height up to the entrance of the collision zone. Finally the third time is the time the particle spent from the entrance of the collision zone up to experiences a collision with the moving wall again. With these steps in mind, we must first find the maximum height that the particle reaches, which is obtained by considering the case 𝑣(𝑡_𝑢)=0 where Δ𝑡_𝑢=𝑡_𝑢−𝑡_𝑛 corresponds to the time the particle spent from the last collision up to the maximum height. It is given by Δ𝑡_𝑢=−ln[𝑔/(𝜂𝑣_𝑛+𝑔)]/𝜂. At this time, the position of the particle is 𝑦_𝑝(𝑡_𝑢)=ℎ_max=𝜖cos(𝜔𝑡_𝑛𝑣)+^𝑛𝜂+𝑔𝜂²𝑔ln(𝜂𝑣_𝑛+𝑔). (2.5) From this maximum height (ℎ_max), the expression that describes the particle's position in the downwards motion is given by 𝑦_𝑝(𝑡)=𝜖cos(𝜔𝑡_𝑛𝑣)+^𝑛𝜂+𝑔𝜂²𝑔ln(𝜂𝑣_𝑛)+𝑔+𝑔𝜂²(1−𝑒^{−𝜂(𝑡−𝑡𝑢)}𝑔)−𝜂(𝑡−𝑡_𝑢),if𝑡≥𝑡_𝑢. (2.6) The next step is to obtain the time the particle spends until the entrance of the collision zone. The instant 𝑡_𝑑 is obtained by matching the condition 𝑦_𝑝(𝑡_𝑑)=𝜖, leading then to the following transcendental function 𝑍(𝑡)=𝑦_𝑝(𝑡)−𝜖 that must be solved: 0=𝑍(𝑡_𝑑)=𝜖cos(𝜔𝑡_𝑛𝑣)+^𝑛𝜂+𝑔𝜂²𝑔ln(𝜂𝑣_𝑛)+𝑔+𝑔𝜂²(1−𝑒^{−𝜂Δ𝑡𝑑}𝑔)−𝜂Δ𝑡_𝑑−𝜖, (2.7) where Δ𝑡_𝑑=𝑡_𝑑−𝑡_𝑢. Since we obtain the solution of the above equation, then we found that the velocity of the particle at the entrance of the collision zone is 𝑣_𝑒=𝑣(𝑡_𝑑), that is written as 𝑣_𝑒=𝑔/𝜂(exp(−𝜂Δ𝑡_𝑑)−1). From there, the expression of the particle velocity is given by 𝑣(𝑡)=𝑒^{−𝜂(𝑡−𝑡𝑑)}(𝑣_𝑒+𝑔𝜂𝑔)−𝜂if𝑡≥𝑡_𝑑. (2.8) To find the instant of the collision of the particle with the moving wall, we must solve the following equation 𝐹(𝑡_𝑐)=𝑦_𝑤(𝑡_𝑐)−𝑦_𝑝(𝑡_𝑐)=0 where 𝑦_𝑤(𝑡) corresponds to the moving wall position and 𝑦_𝑝(𝑡) is the particle position obtained by integration of (2.8). Then we found that 0=𝐹(𝑡_𝑐1)=𝜂(𝑣_𝑒+𝑔𝜂𝑔)−𝜂Δ𝑡_𝑐−1𝜂(𝑣_𝑒+𝑔𝜂)𝑒^{−𝜂Δ𝑡𝑐}−𝜖cos(𝜔(𝑡_𝑛+Δ𝑡_𝑢+Δ𝑡_𝑑+Δ𝑡_𝑐))+𝜖, (2.9) where Δ𝑡_𝑐=𝑡_𝑐−𝑡_𝑑 is the smallest positive root of the above equation for the interval Δ𝑡_𝑐∈[0,2𝜋/𝜔). Thus, in the instant of the collision we have 𝑇_𝐼=𝑣_𝑛+1=−[𝑒^{−𝜂Δ𝑡𝑐}(𝑣_𝑒+𝑔𝜂𝑔)−𝜂]−2𝜖𝜔sin(𝜔𝑡_𝑛+1𝑡),_𝑛+1=𝑡_𝑛+Δ𝑡_𝑢+Δ𝑡_𝑑+Δ𝑡_𝑐, (2.10) where Δ𝑡_𝑐 is obtained from the smallest solution of 𝐹(𝑡_𝑐)=0.

We can see that there is a relatively large number of control parameters, 4 in total, namely, 𝜖, 𝜂, 𝜔, 𝑔 and that the dynamics of the system does not depend on all of them. It is convenient to define dimensionless and more appropriate variables. Thus we define 𝑉=(𝑣_𝑛𝜔)/𝑔, 𝛿=𝜔/𝜂, 𝑎=(𝜖𝜔²)/𝑔, and measure the time in terms of the number of oscillations of the moving wall, 𝜙=𝜔𝑡. With this new set of variables, the mapping that describes the dynamics of the model is given by 𝑉𝑇=^𝑛+1=−[(𝑉^∗_𝑛+𝛿)𝑒^{−Δ𝜙𝑐/𝛿}−𝛿]−2𝑎sin(𝜙_𝑛+1𝜙),_𝑛+1=[𝜙_𝑛+Δ𝜙]mod(2𝜋), (2.11) where the variables 𝑉^∗_𝑛 and Δ𝜙 depend on the kind of collisions that occurs. For successive collisions we have that 𝑉^∗_𝑛=𝑉_𝑛 and Δ𝜙=Δ𝜙_𝑐 where Δ𝜙_𝑐 is obtained as the smallest positive solution of 𝐺(𝜙_𝑐)=0 for Δ𝜙_𝑐∈(0,2𝜋], with 0=𝐺(𝜙_𝑐)=𝑎cos(𝜙_𝑛)+𝛿(𝑉_𝑛+𝛿)(1−𝑒^{−Δ𝜙𝑐/𝛿})−𝛿Δ𝜙_𝑐−𝑎cos(𝜙_𝑛+Δ𝜙_𝑐). (2.12) For the case of indirect collisions, we have that 𝑉^∗_𝑛=𝑉_𝑒 where 𝑉_𝑒=𝛿[exp(−Δ𝜙_𝑑/𝛿)−1] and Δ𝜙=Δ𝜙_𝑢+Δ𝜙_𝑑+Δ𝜙_𝑐 with Δ𝜙_𝑢=−𝛿ln(𝛿/(𝑉_𝑛+𝛿)) corresponding to the upward elapsed time, and Δ𝜙_𝑑 is obtained as the smallest positive root of 𝑍(𝜙_𝑑)=0 with 0=𝑍(𝜙_𝑑)=𝑎cos(𝜙_𝑛)+𝛿𝑉_𝑛+𝛿²𝛿ln(𝑉_𝑛)+𝛿+𝛿²(1−𝑒^{−Δ𝜙𝑑/𝛿})−𝛿Δ𝜙_𝑑−𝑎, (2.13) where Δ𝜙_𝑑 is the elapsed time along the downward direction until the entrance of the collision zone. Finally Δ𝜙_𝑐 is the smallest positive solution of 𝐹(𝜙_𝑐)=0, with 0=𝐹(𝜙_𝑐)=𝑎+𝛿(𝑉_𝑒+𝛿)−𝛿Δ𝜙_𝑐−𝛿(𝑉_𝑒+𝛿)𝑒^{−Δ𝜙𝑐/𝛿}−𝑎cos(𝜙_𝑛+Δ𝜙_𝑢+Δ𝜙_𝑑+Δ𝜙_𝑐). (2.14)

After some straightforward algebra we found that the determinant of the Jacobian matrix is written as det𝜕𝑉_𝑛+1𝜕𝑉_𝑛𝜕𝑉_𝑛+1𝜕𝜙_𝑛𝜕𝜙_𝑛+1𝜕𝑉_𝑛𝜕𝜙_𝑛+1𝜕𝜙_𝑛𝑉=𝐴(^𝑛+𝑎sin(𝜙_𝑛)𝑉_𝑛+1+𝑎sin(𝜙_𝑛+1)), (2.15) where 𝐴 depends on the kind of collisions. For multiple collisions 𝐴=exp(−𝜙_𝑐/𝛿) and for indirect collisions 𝐴=𝛿exp(−[𝜙_𝑐+𝜙_𝑑]/𝛿)/(𝑉_𝑛+𝛿). We note that when 𝛿→∞, which is equivalent to 𝜂→0, then 𝐴→1 and therefore all results for the non dissipative case are recovered.

3. Numerical Results

In this section we will present and discuss our numerical results for the dissipative bouncer model. It is shown in Figure 1(a) the behavior of a chaotic attractor in phase space for the control parameters 𝛿=60 and 𝑎=10. To confirm that the orbit shown in Figure 1(a) is really chaotic, we should then evaluate the Lyapunov exponent and check whether the larger value is positive [19]. It is well known in literature that the Lyapunov exponent is a good tool that can be used to quantify the sensitivity to the initial conditions. They are defined as 𝜆_𝑗=lim_𝑛→∞1𝑛^𝑛_𝑖=1ln|Λ^𝑖_𝑗|,𝑗=1,2, (3.1) where Λ^𝑖_𝑗 are the eigenvalues of the matrix 𝑀=∏^𝑛_𝑖=1𝐽_𝑖(𝑉_𝑖,𝜙_𝑖) with 𝐽_𝑖 representing the Jacobian matrix of the mapping. In the appendix we show all the expressions of the Jacobian matrix for both multiple and indirect collisions.

(a)

(b)

Figure 1 

(a) represents a chaotic attractor for the dissipative bouncer model. The control parameters used were 𝛿=60 and 𝑎=10. (b) corresponds to the convergence of the positive Lyapunov exponent for five different initial conditions for the same control parameters used in (a).

After implementing the algorithm in a computer code, we found that the Lyapunov exponent for the attractor shown in Figure 1(a) is positive as can be seen in Figure 1(b). We have evolved in time five different initial conditions for the control parameters 𝛿=60 and 𝑎=10 and saw that both converge for long time for a constant value of 𝜆=2.443±0.003.

We have also obtained the behavior of the positive Lyapunov exponent for the dissipative bouncer model as a function of the control parameters 𝑎 and 𝛿, as we can see in Figures 2(a) and 2(b). Each point in Figure 2 was obtained as an average of five different initial conditions along the chaotic attractor and evolved up to 6⋅10⁶ collisions of the particle with the moving wall. We can see that the positive Lyapunov exponent grows as a function of 𝑎 as well as a function of 𝛿.

(a)

(b)

Figure 2 

Positive Lyapunov exponent as function of the control parameter: (a) 𝑎 and (b) 𝛿.

Let us now discuss the behavior of the average velocity for chaotic orbits. To obtain the average velocity we consider two different kinds of averages: (1) average over the orbit (sometimes also called as Birkoff's average) and (2) average over an ensemble of initial conditions. The average over the orbit is obtained by evolving for a long time a single initial condition so that we have 𝑉_𝑖=1𝑛^𝑛_𝑗=1𝑉_𝑖,𝑗, (3.2) where 𝑗 corresponds to the number of collisions and 𝑖 corresponds to a sample in an ensemble of 𝑀 different initial conditions. After that, we can average over different initial conditions leading to 1𝑉=𝑀^𝑀_𝑖=1𝑉_𝑖, (3.3) where 𝑀 is the number of different initial conditions used. It is then shown in Figures 3(a) and 3(b) the behavior of 𝑉×𝑛 for different control parameters. We can see that the velocity grows at first and then approaches a regime of saturation for long enough time (similar results for nondissipative dynamics have also been observed [20–22]). We consider now the behavior of the velocity for long time, thus corresponding to the constant plateau. Such plateaus are clear evidence that the phenomenon of Fermi acceleration has been suppressed [14].

(a)

(b)

Figure 3 

Behavior of 𝑉×𝑛 for different control parameter, as shown in the figure.

After a careful analysis of Figures 3(a) and 3(b) we can see that as both the control parameters rise, the value of the average velocity for long time, we call it now as 𝑉_sat, also rises. Such a behavior allows us to suppose that 𝑉_sat∝𝛿^𝛼1𝑎^𝛼2, (3.4) where 𝛼₁ and 𝛼₂ are critical exponents that must be obtained.

If we obtain the behavior of 𝑉_sat as function of both 𝑎 and 𝛿, we can then obtain the critical exponents 𝛼₁ and 𝛼₂. It is shown in Figures 4(a) and 4(b) the behavior of 𝑉_sat as a function of the control parameter 𝑎 and 𝛿. After adjusting a power law fitting, we obtain that the critical exponents are 𝛼₁=0.299±0.003 and 𝛼₂=(0.747±0.007)≅3/4.

(a)

(b)

Figure 4 

Behavior of 𝑉sat as function of (a) 𝑎 and (b) 𝛿. After a power law fitting we obtain that 𝛼1=0.299±0.003 and 𝛼2=(0.747±0.007)≅3/4.

Let us now discuss the consequences of the results obtained above. If we increase the value of the control parameter 𝑎, it corresponds to enlarge the ratio of accelerations, then allowing the particle to acquire more energy upon collision with the moving wall. This result is then confirmed since the critical exponent obtained was 𝛼₂=0.747±0.007. Similar result can be obtained as function of the control parameter 𝛿. However, we must emphasize that as the control parameter 𝛿 rises, the straight of the dissipation decreases since 𝛿∝1/𝜂. So in the limit of very large values of 𝛿, almost no dissipation is present. For such a regime, it is then expected that Fermi acceleration should holds. Our results confirm such a supposition since the critical exponent 𝛼₁=0.299±0.003 leading then to the phenomenon of Fermi acceleration for the regime of 𝛿→∞. However for finite values of 𝛿, suppression of Fermi acceleration is observed. This result enlarges the condition where FA is suppressed since it is already known [14] that inelastic collision also suppress Fermi acceleration.

4. Conclusions

We have studied the bouncer model with the presence of a drag force that is proportional to particle's velocity. The model was described by using a two-dimensional nonlinear map, obtained by the solution of a Newton's law of motion. We show that chaotic attractors are present in the system and we have characterized them by use of the positive Lyapunov exponent. We have shown that the introduction of a drag force leads the particle to experience a limited energy gain. Then as to answer the question proposed in the title, the introduction of a drag force proportional to the particle's velocity can indeed be assumed as a mechanism to suppress Fermi acceleration. Such suppression is topologically demostrated by the determinant, of the Jacobian matrix. Since the pre-factor A in the determinant of the Jacobian matrix is less than the unity, attractors are present in the phase space; thus unlimited energy growth is not feasible anymore.

Appendix

We present in this appendix all the expressions of the partial derivatives of the Jacobian matrix. We start with the expressions obtained for the multiple collisions 𝜕𝑉_𝑛+1𝜕𝑉_𝑛=𝑒^{−Δ𝜙𝑐/𝛿}[𝜕𝜙_𝑐𝜕𝑉_𝑛𝑉(1+^𝑛𝛿)−1]−2𝑎cos(𝜙_𝑛+1)𝜕𝜙_𝑛+1𝜕𝑉_𝑛,𝜕𝑉_𝑛+1𝜕𝜙_𝑛=𝑒^{−Δ𝜙𝑐/𝛿}𝜕Δ𝜙_𝑐𝜕𝜙_𝑛𝑉(1+^𝑛𝛿)−2𝑎cos(𝜙_𝑛+1)𝜕𝜙_𝑛+1𝜕𝜙_𝑛,𝜕𝜙_𝑛+1𝜕𝑉_𝑛=𝜕𝜙_𝑐𝜕𝑉_𝑛,𝜕𝜙_𝑛+1𝜕𝜙_𝑛=𝜕𝜙_𝑐𝜕𝜙_𝑛+1, (A.1) where 𝜕𝜙_𝑐/𝜕𝑉_𝑛 and 𝜕𝜙_𝑐/𝜕𝜙_𝑛 are obtained, respectively, via implicit derivative of 𝐺(Δ𝜙_𝑐) with respect to 𝑉_𝑛 and 𝜙_𝑛. Their expressions are 𝜕𝜙_𝑐𝜕𝑉_𝑛=𝛿(𝑒^{−Δ𝜙𝑐/𝛿}−1)𝑒^{−Δ𝜙𝑐/𝛿}(𝑉_𝑛+𝛿)−𝛿+𝑎sin(𝜙_𝑛+1),𝜕𝜙_𝑐𝜕𝜙_𝑛=𝑎sin(𝜙_𝑛)−𝑎sin(𝜙_𝑛+1)(𝑉_𝑛+𝛿)𝑒^{−Δ𝜙𝑐/𝛿}−𝛿+𝑎sin(𝜙_𝑛.+1) (A.2)

We now obtain the expression of the Jacobian terms for the indirect collisions 𝜕𝑉_𝑛+1𝜕𝑉_𝑛=𝑒^{−Δ𝜙𝑐/𝛿}[𝜕Δ𝜙_𝑐𝜕𝑉_𝑛𝑉(1+^𝑒𝛿)−𝜕𝑉_𝑒𝜕𝑉_𝑛]−2𝑎cos(𝜙_𝑛+1)𝜕𝜙_𝑛+1𝜕𝑉_𝑛,𝜕𝑉_𝑛+1𝜕𝜙_𝑛=𝑒^{−Δ𝜙𝑐/𝛿}[𝜕Δ𝜙_𝑐𝜕𝜙_𝑛𝑉(1+^𝑒𝛿)−𝜕𝑉_𝑒𝜕𝜙_𝑛]−2𝑎cos(𝜙_𝑛+1)𝜕𝜙_𝑛+1𝜕𝜙_𝑛,𝜕𝜙_𝑛+1𝜕𝑉_𝑛=𝜕Δ𝜙_𝑢𝜕𝑉_𝑛+𝜕Δ𝜙_𝑑𝜕𝑉_𝑛+𝜕Δ𝜙_𝑐𝜕𝑉_𝑛,𝜕𝜙_𝑛+1𝜕𝜙_𝑛=1+𝜕Δ𝜙_𝑑𝜕𝜙_𝑛+𝜕Δ𝜙_𝑐𝜕𝜙_𝑛, (A.3) where 𝜕𝜙_𝑑/𝜕𝑉_𝑛, 𝜕𝜙_𝑑/𝜕𝜙_𝑛, 𝜕𝜙_𝑐/𝜕𝑉_𝑛, and 𝜕𝜙_𝑐/𝜕𝜙_𝑛 are obtained by implicit derivative of 𝑍(Δ𝜙_𝑐) and 𝐹(Δ𝜙_𝑐) with respect to 𝑉_𝑛 and 𝜙_𝑛. Their expressions are 𝜕Δ𝜙_𝑑𝜕𝑉_𝑛=𝜕Δ𝜙_𝑢/𝜕𝑉_𝑛−1𝑒^{−Δ𝜙𝑑/𝛿},−1𝜕Δ𝜙_𝑑𝜕𝜙_𝑛=𝑎sin(𝜙_𝑛)𝛿(𝑒^{−Δ𝜙𝑑/𝛿},−1)𝜕Δ𝜙_𝑐𝜕𝑉_𝑛=𝛿(𝜕𝑉_𝑒/𝜕𝑉_𝑛)(1−𝑒^{−Δ𝜙𝑐/𝛿})+𝑎sin(𝜙_𝑛+1)((𝜕Δ𝜙_𝑢/𝜕𝑉_𝑛)(𝜕Δ𝜙_𝑑/𝜕𝑉_𝑛))𝛿−𝑒^{−Δ𝜙𝑐/𝛿}(𝑉_𝑒+𝛿)−𝑎sin(𝜙_𝑛+1),𝜕Δ𝜙_𝑐𝜕𝜙_𝑛=𝛿(𝜕𝑉_𝑒/𝜕𝜙_𝑛)(1−𝑒^{−Δ𝜙𝑐/𝛿})+𝑎sin(𝜙_𝑛+1)(𝜕Δ𝜙_𝑑/𝜕𝜙_𝑛+1)𝛿−𝑒^{−Δ𝜙𝑐/𝛿}(𝑉_𝑒+𝛿)−𝑎sin(𝜙_𝑛+1), (A.4) where𝜕Δ𝜙_𝑢/𝜕𝑉_𝑛, 𝜕𝑉_𝑒/𝜕𝑉_𝑛, and 𝜕𝑉_𝑒/𝜕𝜙_𝑛 are the implicit derivative of Δ𝜙_𝑢 with respect to 𝑉_𝑛 and the the implicit derivative of 𝑉_𝑒 with respect to 𝑉_𝑛 and 𝜙_𝑛, respectively. They are written as 𝜕Δ𝜙_𝑢𝜕𝑉_𝑛=𝛿𝑉_𝑛,+𝛿𝜕𝑉_𝑒𝜕𝑉_𝑛=−𝑒^{−Δ𝜙𝑑/𝛿}𝜕Δ𝜙_𝑑𝜕𝑉_𝑛,𝜕𝑉_𝑒𝜕𝜙_𝑛=−𝑒^{−Δ𝜙𝑑/𝛿}𝜕Δ𝜙_𝑑𝜕𝜙_𝑛. (A.5)

Acknowledgments

F. A. de Souza and L. E. A. Simões thank FAPESP, Brazilian agency. E. D. Leonel gratefully acknowledges support to CNPq, FAPESP and FUNDUNESP, Brazilian agencies.