On the escape of a resonantly excited couple of particles from a potential well

The escape dynamics of a damped system of two coupled particles in a truncated potential well under biharmonic excitation are investigated. It is assumed that excitation frequencies are tuned to the modal natural frequency of the relative motion and to the modal frequency of the centre of mass on the bottom of the potential well. Although the escape is essentially a non-stationary process, the critical force strongly depends on the stationary amplitude of the relative vibrations within the pair of masses. The characteristic escape curve for the critical force moves up on the frequency-escape threshold plane with increasing relative vibrations, which can be interpreted as a stabilizing effect due to the high-frequency excitation. To obtain the results, new modelling techniques are suggested, including the reduction in the effect of the high-frequency excitation using a probability density function-based convolution approach and an energy-based approach for the description of the evolution of the slow variables. To validate the method, the coupled pair of particles is investigated with various model potentials.


Introduction
Escape from a potential well is a classic problem arising in various fields of natural science and engineering [1][2][3][4][5]. A wide range of physical phenomena from celestial mechanics to the dynamics of molecules has been investigated in the literature, including topics such as energy harvesting [6], transient resonance dynamics of oscillatory systems [7,8], the physics of Josephson junctions [9] and even the capsize phenomenon of ships and vessels [3,10]. Another important topic relating to the escape phenomenon is the dynamic pull-in in microelectromechanical systems (MEMS) [11][12][13][14][15][16].
The first investigations of the forced escape phenomenon started in 1940 by Kramers regarding the thermal activation of chemical reactions [17,18]. Although much research has taken place in the last 80 years, one still encounters important unresolved issues when considering the escape process [19].
For constant forcing, escape may occur because of the slow variation of system parameters and the consequent bifurcations of the steady-state regimes of the response [2,11,12].
Paper [3] addressed the escape dynamics of harmonically forced and damped particles in three model potential wells. To give an analytic expression for the escape criterion, the steady-state response of the particle was used, although this criterion is not exact and empirical corrections are performed when necessary.
A common feature of the escape phenomenon is numerically demonstrated in paper [3]. For the investigated potentials, the critical forcing amplitude curve has a sharp minimum value at a frequency value smaller than the frequency of small oscillations in the potential well. Similar escape curves are found in papers investigating the problem of dynamic pull-in in MEMS devices [14,15].
Further similar patterns were found when investigating the safe basins of attraction for several dynamical systems [20,21]. The sharp minimum value found in many papers allows the reader to formulate a conjecture about the mentioned property of the curve as a typical characteristic of the escape phenomena. The conjecture was further examined for different model potentials under harmonic forcing in the papers [22,23].
These papers addressed the escape dynamics of the forced particle in conditions of 1:1 resonance. The transient dynamics of the system can thus be explored in terms of the slow flow on the resonance manifold (RM). Exploration of special trajectories on the RM corresponding to different initial conditions (IC) allows for the identification of two distinct mechanisms of approaching the escape threshold while varying the system parameters. In the first mechanism, the maximum mechanism (MM), the phase trajectory only achieves the boundary of the potential well. In the second mechanism, the saddle mechanism (SM), the special trajectory achieves the saddle point on the RM. Then, the response amplitude drastically grows, and the particle escapes. The combination of these two mechanisms provides an understanding and prediction of the sharp minimum on the frequency-amplitude plane without any empirical corrections. This approach somewhat resembles recent studies of transient phenomena in systems of coupled oscillators in terms of limited phase trajectories (LPT) [24,25].
In the present paper, a more complex situation is considered. A system of two coupled particles (cf. Fig. 1a) is excited by a sum of two harmonic components. The first one corresponds (in the linear case) to the internal mode of the coupling corresponding to their relative oscillations. The second resonance is close to the 1:1 resonance for the centre of mass of the system. First, after describing the governing equations of the system (Sect. 2), the reduction in the degrees of freedom of the system follows (Sect. 3.1). Then, the arising equations are further simplified using the probability-based convolution to calculate the effective potential (Sect. 3.2). Subsequently, further simplification is achieved by the averaging technique to obtain equations for the slowflow variables (Sect. 3.3). After giving an explicit solution of the slow-flow equations, the presentation of some numerical results and the discussion of the results follows (Sect. 4). In Sect. 5, the results presented in the paper are summarised.

Description of the model
The setting investigated in the paper is depicted in Fig. 1a. Two bodies of mass m 1 and m 2 coupled with a linear spring of stiffness k and a linear damper with the damping coefficient c are located in the corresponding potential wells V 1 (x) and V 2 (x). m 1 is excited by the superposition of two harmonic forces with the amplitudes A 1 , A 2 , the excitation frequencies ω 1 , ω 2 and the phases β 1 , β 2 . The potentials act in a way that both bodies have the same linearized eigen- . For example, the truncated parabolic potential well has the following form: The equations of motion are given by where V (x) := ∂ V (x) ∂ x . Using the following dimensionless values Fig. 1 The basic model of the two particle system and the dimensionless potential U 0 (q) (see Eq. 9) with the notationq = d dτ q and U 0 (q) = d dq U 0 . Equation (2) can be rewritten in dimensionless form too, although the signs of the masses remain in the equations for clarity: where C := c ω 0 and K := k To decouple the equations from each other, a new pair of coordinates, (y 1 , y 2 ) is introduced for the position of the centre of gravity and for the relative displacement of the two bodies, respectively: Then, which substituted back into (4) results in with 1 m := 1 m 1 + 1 m 2 and μ := m 1 m 1 +m 2 .

Model reduction
The equations of motion are too complicated in their current form for further analytic investigations. In the following, some assumptions are introduced regarding the properties of the problem. To concentrate on the effects caused by the coupled bodies in contrast to a single body, the truncated, purely quadratic potential (1) is mainly investigated; however, most of the calculations are presented in a general form so that the procedure described below is applicable to arbitrary potentials.

Reduction in the degrees of freedom
In the following, F 2 F 1 and 2 1 will be assumed; furthermore, C is assumed to be not small and K 1. The approach described in [26,27] can be applied to simplify the equations. The asymptotic behaviour of the solution can be approximated with a small error by solving Eq. (8), neglecting the small coupling term in Eq. (7).
Since the action of the potential well on y 2 is always present in the coupling term only, the reduction of (8) to a damped, linear ODE of second order is possible for arbitrary potentials.

Example-truncated parabolic potential
In dimensionless form, the truncated parabolic potential takes the following form Neglecting the coupling term and solving the remaining linear ODE for y 2 would be possible here. However, in the above simple case, it turns out to be a slightly more accurate heuristic for the solution of (8), if only the case is considered when both of the particles are inside of the potential well. Then, Interesting effects appear when 2 is chosen close to the resonance frequency of the relative movements of the particles, i.e., Then, the effect of F 1 compared to F 2 is negligible, and the particular solution of equation (8) can be written as with the vibration amplitude and with the phase shift, γ 2 for which the explicit expression is omitted here, because it is not required for further analysis. Due to the non-small damping, the homogeneous solution y 2,h (τ ) decays rapidly, so its effect on the behaviour of the solution is also negligible.
Equation (7) can now be simplified by the elimination of the variable y 2 (τ ).
It should be noted that the method is applied only to the simplest case of truncated quadratic potential here, but according to [26,27], it is applicable for all kinds of nonlinear potentials as well when the coupling between the particles is much stronger than the effect of the potential on their relative movements.

Effective potential: PDF Approach
After the approximation, applied above, there are three different dimensionless time scales in equation (7): • the very fast time scale includes the fast parametric excitation U 0 (y 1 (τ ), 2 τ + β 2 ) and the fast force excitation F 2 sin( 2 τ + β 2 ). In the next step, this time scale will be eliminated; • in addition, the fast time scale including the motion of the centre of mass of the bodies, y 1 (τ ) and the excitation F 1 sin( 1 τ + β 1 ) is also present in the dynamics of the system; • modulations in the vibration amplitude of the centre of gravity, A (Eq. 32) and its phase difference from the excitation F 1 sin( 1 τ + β 1 ) can be described by the slowest time scale.
To distinguish the last two time scales from each other, the variable t := ετ is altered to denote the slow time, whereas the fast time is denoted by τ .
To eliminate the fastest time scale, some changes in the modelling approach are helpful. Instead of always considering the the exact position of the bodies, one suggests a statistical approach in which a particle has a probability to be found at an exact position. The movement of the bodies around their common centre of gravity, y 1 , is sinusoidal, so the corresponding cumulative distribution function (CDF) is given by the arcsine distribution The resulting probability density function (PDF) is After appropriate scaling and normalisation of this PDF, one obtains the density distribution for the bodies depending on y 0 and y 1 (t) An example for ρ(x, y 1 ; y 0 ) is depicted in Fig. 1b. The resulting force acting on the center of gravity can now be calculated through the following integration: (18) from which the effective potentialŨ 0 (y 1 ; y 0 ) is obtained through integration with the conditioñ U 0 (−∞; y 0 ) ! = 0. In general (for arbitrary potentials), S(y 1 ) cannot be calculated in an analytically closed form. In the following, the calculations are left generic, but for a specific one, the used quantities are determined numerically. (For the truncated quadratic potential, the integrals can be evaluated analytically, but it is quite elaborate and does not provide any significant advantage.) A possible way to handle the problem is piecewise spline fitting on a fine-resolved discrete data set. This step is not par-ticularly complicated when using MATLAB ® 1 . For functions calculated in this way, see Fig. 2. 3.3 Obtaining differential equations for the "slow-flow" variables One may be interested in the behaviour of the solution if the excitation frequency is close to the natural frequency determined by the potential well. Using the above assumption, Eq. (7) may be rewritten as follows: Here, S and S * also depend on the parameters y 0 , m 1 and m 2 . Equation (19) is handled as follows. It is close to a conservative onë with F 1 = O(ε) and F 2 = O(ε). The unperturbed equation enables the conversion to the variables characterizing the total energy of the system and the corresponding phase: The choice of the sign in the last equation is determined by the sign ofẏ 1 .
Although F(τ ) is not small due to the term F 2 , its effect on the energy of the system is small because 2 1. The governing equations of the new variables in the perturbed system are slow: These equations are suitable for averaging with respect to the explicit time [28,29]. The corresponding equations of the first-order approximation can be formally written as follows: where denotes averaging with respect to explicit time. However, to evaluate these equations, the func-tionŨ 0 (y 1 ) has to be expressed explicitly in terms of the new variables (E, θ). In general, it is impossible. Thus, an approximate heuristic approach is applied for obtaining analytic predictions. Let us suppose that the solution of (20) can be approximated by almost harmonic oscillations:

Here, A(t) is the amplitude of the oscillations and (t)
is the phase shift. Note that these variables depend on the energy of the system, which changes slowly according to (23)(24). According to the transformation (21-21), we can replace the total energy with the level of the potential energy at the maximal deflection, i.e. at y 1 = A: This means that the time derivative of the energy can be calculated to the first-order approximation through the time derivative of the amplitude: However, the maximal velocity is achieved as the system passes through the minimum of the potential energy. Hence, the amplitude of the velocity can be calculated as follows: Inserting Eqs. (30) and (31) into (25), the following governing equation for the amplitude is obtained: Requiring the consistency of the transformation (27) and (28), the following equation for the phase shift is obtained: Assuming that the low excitation frequency is close to the eigenfrequency of the system on the bottom of the potential well, a small discrepancy can be introduced: Then, Eq. (33) can be rewritten as follows: In addition, according to the introduced non-dimensional parametersŨ 0 (y)| y→0 = y. Hence, the second and third terms in this equation can be combined: −μF(τ ) sin( 1 τ + β 1 (t)), In general, the second term in this equation is not small. Hence, averaging of this equation is not completely justified. Nevertheless, we replace Eq. (36) with the averaged one and check the quality of this approximation by comparing its results with the results of the direct numeric simulations of the full system. The averaged Eq. (36) is: with where ω(A) is the amplitude-dependent eigenfrequency of the potential. Based on [1], the time period of a particle's oscillation with mass m and with the given total energy E in a potential wellŨ (x) without excitation, i.e. having a conservative system, can be calculated as follows: Here, x 1 and x 2 are the energy-dependent maximal displacements in the potential with the given total energy of the particle E. Then, .
The introduced term ω * (A) describes the effect of nonlinearity on the frequency of free oscillations in the potential well.
For several examples of the potential, this function is displayed in Figs. 2d and 7d. The terms describing the averaged effect of the external excitation can be easily calculated: Inserting (41) into Eqs. (32) and (37), the final form of the first-order approximation is obtained: with the initial conditions The initial values of y 1,0 andẏ 1,0 can be converted to find 0 and A 0 . To achieve this, consider the solution of the linearized differential equation of the center of mass neglecting the small right-hand sidë Then, the solution is Comparing these equations to (27 -28) and using the initial conditions (46-47) for τ = 0, one obtains Due to the linearization, the conversion is exact only in the case of a purely quadratic potential; still, for smaller deviations from the quadratic potential, it gives a reasonably good approximation of the initial conditions.

Solution by integrating factors
It is possible to write the solution to Eqs. (43)-(44) in a closed form with the use of integrating factors. Let us define Dividing Eq. (44) by Eq. (43) and reorganising the terms, one obtains An appropriate integrating factor μ(A) is given by . can be observed by the LPT's straight alignment. The effective potential has, however, a nonlinear force that requires a minimum amount of exciting force to let the escape happen Multiplying (53) by μ(A), the partial derivatives of the first integral C(A, ) can be determined The first integral is then which is due to its complexity being hardly solvable for any given U 0 . Although this is an analytic solution in a closed form, its evaluation is even more elaborate than direct numerical integration. Next, we present the numerical results by direct simulation.

Numerical results and discussion
Equations (43)-(44) are easy to solve numerically. The structure of the resonance manifold is shown in Fig. 3 for the parameter values m 1 = m 2 = 1 and y 0 = 0.2. Figures 4 and 5 show the critical force over the excitation frequency diagram for the case of the truncated quadratic potential with a relative amplitude of the bodies of y 0 = 0.2 for zero initial conditions, which corresponds to the LPT. Given the initial conditions q 1 ,q 1 , q 2 andq 2 for the particles, one can determine y 1,0 anḋ y 1,0 and consequently A 0 and 0 as well. This allows us to investigate the escape scenarios for nonzero initial conditions as well. For nonzero initial conditions, the diagram of the critical force over the excitation frequency changes its shape. Figure 6 shows a typical scenario when the escape cannot occur, although a single body would escape for any positive excitation force, as the excitation takes place exactly on the natural frequency of the system. However, the nonlinear boundary of the effective potential does not allow such an escape scenario. Figure 7 shows the graphs of some effective quantities for the reduced system in the case of asymmetrical mass distribution, m 1 = 2m 2 = 2 compared to the original graphs derived from the truncated quadratic potential. transferred to the system is positive and thus increases the amplitude, consequently leading to escape if the phase of the excitation does not differ by more than ± π 2 from the phase of the motion. However, a phase difference larger than ± π 2 means that the system is delivering work to the exciter, and so A decreases.
Due to the progressive or degressive characteristic of the force resulting from the underlying potential, ω(A) is either a monotonously increasing or decreasing function. Therefore, for increasing A, the body swings faster or slower, accumulating an increasing phase difference from the phase of the excitation. If this difference reaches π 2 or − π 2 , the energy transferring from the excitation to the energy of the body turns over to the energy transferring from the body to the exciter, and no escape occurs, at least for the moment. In the case of the truncated quadratic potential, the degressivity of the effective eigenfrequency has a very remarkable explanation. Once one of the bodies gets out of the potential well, a restoring force is no longer acting on it; thus, the restoring force on the center of mass drops to half of the original value (for m 1 = m 2 ). This results in softening spring characteristics as the bodies spend increasingly more time outside of the potential well for the increasing oscillation amplitude of the center of mass of the bodies, A.

Some further remarks
It can also be noted that the convolution operator used to determine the effective potential increases the regularity of the effective potential function. The kernel of the integral, the PDF of the arcsin distribution (16), is a Hölder continuous function with the exponent of the Hölder condition α = 0.5.
In the limiting case, when y 0 → 0, the problem reduces itself to a single body problem, which has been investigated by many authors. This agrees with the model, since the integral kernel tends to a Dirac distribution for y 0 → 0 by having a PDF, whose integral is by definition 1 and whose support's length tends to 0.

Conclusions
The dynamics of a couple of strongly coupled particles in a potential well under biharmonic resonant excitation have been analytically and numerically investigated. Under the assumption of non-small damping, the system's dynamics have been reduced to one degree of freedom describing movements of the system's centre of mass in the modified effective potential well, which takes into account the smoothened probabilistic relative vibrations of the particles. The subsequent application of the averaging procedure for the system being close to a conservative one has enabled very accurate prediction of the escape behaviour, especially for small values of the critical excitation force. Both the analytic and numerical results clearly show the stabilizing effect of the high-frequency relative movements of the pair of particles on their escape behaviour. This effect increases with increasing amplitude of the relative vibrations. Although the obtained analytic results demonstrate excellent agreement with the results of direct numerical simulations of the full system, the investigated problem and the applied approach open a whole class of theoretical and practical questions that must be addressed in the future.
Do two bodies with high-frequency excitation always stabilise the escape behaviour? If not, for what kind of potentials does this not occur? What happens to potentials where ω(A) is not monotonous? To what extent could the probability-based modelling approach be applied to describe high-frequency excited physical systems? How large is the uncertainty of the solution caused by this approach? Is there a way to synthesize potentials to shape the form of the F crit ( 1 )-curve?
The authors hope that these open questions can attract the attention of readers to the fascinating problems of escape dynamics for groups of coupled particles and even more complex structures, opening a broad field for analytic investigations and creative applications.
Funding Open Access funding enabled and organized by Projekt DEAL.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.