1 Introduction

Systems of coupled oscillations are largely studied on account of their wide possibility of application in many significant branches such as mechanics and medical and biological sciences. The corresponding mathematical problem is in no way easy to handle when all the effects are overlapped. Here, we propose a basic situation which will be discussed from the mathematical point of view.

The main question we deal with is the feasibility of in-phase or antiphase synchronization globally with respect to the time when no external forces (escapement) forcing the free oscillation are contemplated. Actually, in most of the models proposed in the literature, an escapement device is always used in the system, and it is somehow claimed to be the unique reason why the synchronization occurs. Thus, what the mathematical problem entails about the synchronization whenever the escapement is removed? This is our point.

We mainly take care of the mathematical path drawn by the equations of motion, aiming at developing the analytical scheme, even if in a simplified situation. We first formulate the mathematical model by allowing very general features of the mechanical phenomenon, admitting different sizes of masses, lengths of the pendula and including escapement conditioning the oscillations. This will supply a ground in order to make a brief comparison with some significant models proposed in the literature.

Instead of obtaining information directly via a numerical simulation approach, we rather aim at an analytical method of locating the eigenvalues linked with the damping of the system. The advantage is the prospect of recognizing some ranges of the parameters entering the phenomena which are able to predispose the system to the synchronization, i.e., this will be performed in Section 4.

On the other hand, some expected results are confirmed, as the unattainability of the in-phase synchronization in the absence of escapement (see Section 3) .

It must be said that most of the proposed models are used in almost exclusively numerical simulations in order to obtain information from the calculated data. In our mind, it is worth the effort investigating deeper the mathematical problem and assembling a definite basis for forthcoming numerical applications. This is the reason why we are determined to keep separated the two methods and to produce numerical simulations in a second time, by starting and taking advantage from the theoretical results here achieved.

The first step, introduced in the next section, is the mathematical formulation of the model, which will be achieved by the Lagrangian's method of deducing the equations of motion.

2 Mathematical model

The system we are going to study can be realized either by the device I or the device II shown in Fig. 1. For the device I, the apparatus consists of two pendula whose pivots (points A and B) are fasten on a horizontal and homogeneous beam with the mass $m_0$, and denote $P_0$ as the centre of mass. The beam is placed on a pair of rollers of the radius $R$. The massive bobs ($m_1$ and $m_2$) are suspended at the extremities $P_1$ and $P_2$ of the massless rods (lengths $\ell_1<R$ and $\ell_2<R$) and can oscillate on the vertical plane containing the beam. The distances $d_1$ and $d_2$ of $P_1$ and $P_2$ from $P_0$ are larger than $\ell_1$ and $\ell_2$, in order to avoid hits. A spring of stiffness $k$ whose mass and length at rest are negligible, is attached at one of the extremities of the beam.

In the device II, the pendula are suspended by means of a T-shaped support, with the negligible mass fixed on the beam which oscillates at a lower height.

If the physical quantities $m_0$, $m_1$, $m_2$, $\ell_1$, and $\ell_2$ are the same in the devices I and II, it is immediate to realize that the mathematical problem is identical (actually, the lengths $h$, $R$, $d_1$, and $d_2$ do not enter the equations of motion). Moreover, not even from the dynamical point of view, the two systems are different. Actually, active forces are the same, assuming that the rolling friction in the apparatus I is proportional to ${\dot \Psi}$, with $\Psi$ being the rotation angle of the roll, then, it is proportional to ${\dot x}$, with $x$ being the abscissa of $P_0$ (namely, $x-2R\Psi$ is constant), just like in the apparatus II.

2.1 Equations of motion

In this first part, the equations of motion are achieved by using a Lagrangian approach.

We find it suitable to draw synthetically the path of the mathematical scheme. Actually, one of the motivations inducing us to express a viewpoint about the problem is an accurate production of the mathematical scheme. As a matter of fact, some incongruencies can be remarked upon the proposed models, as we will shortly discuss in the next paragraph.

The cartesian coordinate system is fixed in order that the mechanism is contained in the vertical plane $y=0$, the beam swings along the $x$-axis, the $z$-axis is upward-vertically directed, and the origin $O$ corresponds to the fixed extremity of the spring.

We choose Lagrangian coordinates as ${q}=(x, \theta_1, \theta_2) $, where $\theta_1$ and $\theta_2$ are the amplitudes of oscillation with respect to the downward-vertical direction, and $x$ is the abscissa of $P_0$. The representative vector of the discrete system $(P_0, P_1, P_2) $ in terms of them is, for the device I,

and $0$ at the third position, $h$ replacing $2R$ for the device II. In both cases, the Lagrangian function ${\cal L}=T+U$, where $U$ is the potential of the elastic force and of gravity, i.e.,

(1)

where $m=\sum\limits_{j=0}^2 m_j$ is the total mass.

As for the friction forces, if the damping is formulated as ${\Phi}_{P_i}=-\beta_i{\dot P}_i (i=0, 1, $ and 2) , the Lagrangian component is

(2)

with $\beta=\sum\limits_{j=0}^2\beta_j$. However, if one assumes that the pendula run into damping only along the rotational direction ${e}_{\theta_i}=(\cos \theta_i, 0, \sin \theta_i)$, $i=1$ and 2, then the generalized friction force can be reduced to

(3)

The mechanism of escapement can be modelled by introducing a moment in the ${j}$-direction (i.e., it is orthogonal to the plane of the system), exerting a force $f_i(P_i, {\dot P}_i, t){e}_{\theta_i}$ on each $P_i$, $i=1$ and 2. The corresponding generalized force is

(4)

Assuming (2) , the aligns of motion $\frac{\rm d}{{\rm d}t}(\nabla_{\dot {q}}{\cal L})- \nabla_{q}{\cal L}={F}^{({q})}+{\Phi}^{({q})}$ are

(5)

2.2 Comparison with some models

By reviewing very briefly the mathematical formulation of some models existing in literature, we have the specific intention of remarking that

(i) If (3) is accepted to hold, then $\beta_0$ replaces $\beta$ in the first equation, and the terms $-\ell_i \beta_i{\dot x}\cos \theta_i$, $i=1$ and 2, have to be omitted in the second and third equation. However, the term $-\sum\limits_{j=1}^2\ell_j \beta_j{\dot \theta}_j\cos \theta_j$ of the first equation is present anyway.

(ii) Even if the escapement mechanism operates on the pendula vie the force (4) , the terms containing $f_1$ and $f_2$ are present in any case in the first equation of (5) concerning $x$, and they affect the motion of the beam.

In Ref. ^[1], where the apparatus I is tested, the escapement is formulated by means of a step function, depending on the amplitude of a threshold angle. In Ref. ^[2], the apparatus II is subject to the inversion of direction of the angular velocity (escapement mechanism) at a critical value of the angle. Also, in Ref. ^[3], the mathematical problem is formulated for two driven pendula (although the description of the experimental setup refers to a couple of metronomes), but (4) is expressed via a continous function.

The cited models are undoubtedly significant and useful for the exhibited experimental and numerical results. However, we remark the lack of the terms in the first equation of (5) , pointed out in (i) and (ii) just above. This aspect should not be unimportant, mainly when analytical results are pursued, as in our investigation.

An experimental device different from I and II described at the beginning of the section consists in placing two masses at each of the pivotal points, and let them oscillate horizontally, by means of a spring connecting the two points (hence, the beam is removed), for example, the interaction mechanism is studied in Ref. ^[4], and the analytical problem is essentially the same with (5) , since the centre of mass of the attachment points solves the first of such system with $k=0$.

Finally, in the system studied in Refs. ^[5] and ^[6], for example, the two pendula are coupled by a spring connecting some intermediate points ($Q_1$ and $Q_2$) of the two sticks supporting the weights (Kumamoto coupled pendula). The spring-coupled pendula are also proposed as a basic model for the neutrino oscillation.

Letting $d$ be the constant distance between the pivotal points and presuming that the length at rest of the spring is $d$, we write the Lagrangian function as $\frac{1}{2} m_1 \ell_1^2 {\dot \theta}_1^2+\frac{1}{2}m_2 {\dot \theta}_2^2 +m_1 g \ell_1 \cos \theta_1 +m_2 \ell_2 g \cos \theta_2 +\frac{1}{2}k( |Q_1(\theta_1) -Q_2(\theta_2) |-d)^2$, where

in which $\ell_{0, 1}$ and $\ell_{0, 2}$ are the distances between the pivots and the intermediate points.

The problem is two-dimensional, and the corresponding equations of the motion for the two angles, by assuming (3) , $\ell_1 f_1 =u(t)$, and $f_2=0$ (the external torque acting only on one pendulum, see Ref. ^[6]), are

which slightly differ from the dynamics formulated in Ref. ^[6], where $\ell_{0, 1}=\ell_{0, 2}$, and $d$ is neglected.

2.3 Mathematical problem for $\sigma$ and $\delta$

The main purpose is to investigate the existence of solutions to the system (5) such that one of the quantities

(6)

tends to zero when $t\rightarrow +\infty$. If $\delta\rightarrow 0$ and $\sigma\rightarrow 0$, then the system will proceed to the in-phase synchronization and antiphase synchronization, respectively. We stress that the mathematical investigation is focused on the solutions showing the in-phase or anti-phase synchronization globally with respect to the time.

Moving toward the more expressive coordinates (6) and setting ${y}=\left( \begin{array}{c} x \\ \sigma \\ \delta \end{array}\right)$, one has ${y}=L{q}$, where $L=$ $\left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 1 \\ 0 & 1 & -1 \end{array}\right)$ is the linear change of coordinates. Thus, writing again the equations of motion by taking account of $ \frac{\rm d}{{\rm d}t}(\nabla_{\dot {y}}{\tilde {\cal L}})- \nabla_{y}{\tilde {\cal L}} = L^{-1} ( \frac{\rm d}{{\rm d}t}(\nabla_{\dot {q}}{\cal L})-\nabla_{q}{\cal L})$, ${F}^{({y})}+{\Phi}^{({y})}=L^{-1}({F}^{({q})}+{\Phi}^{({q})})$ with ${\tilde {\cal L}}({y}, {\dot {y}})={\cal L}(L^{-1} {y}, L^{-1}{\dot {y}})$ and $L^{-1}={\left(\begin{array}{ccc} 1 & 0 & 0\\ 0 & 1/2 & 1/2 \\ 0 & 1/2 & -1/2 \end{array}\right)}$, we obtain

(7)

where

(8)

and

(9)

It is worth noticing that

(i) If (3) holds instead of (2) , then $\beta_0$ replaces $\beta$ after the equals sign in the first equation, and each last term of the second and third equation (containing ${\dot x}$) has to be removed.

(ii) In case of the identical pendula ($m_1=m_2$, $\ell_1=\ell_2$), all the quantities in (9) with the superscript vanish, so that each term in (7) containing them must be cancelled.

(iii) Likewise, if also $f_1=f_2$, additional simplifications are evident.

(iv) Solving (7) explicitly with respect to the second order derivatives, one gets

(10)

where $\Theta= m-\frac{1}{2}((m_1+m_2) +\Psi_1 (m_1+m_2, m_2-m_1) )$, and

(11)

We notice that the escapement (4) is missing in the first equation governing the beam's motion. Furthermore, even in the general case of different masses and lengths, the motion of $\delta$ (the third equation) is more disentangled from the rest of the system than $\delta$. This property will be clearer afterward but actually can be prefigured here, if we imagine to replace in (10) the function (8) with the second order Taylor polynomial,

(12)

2.4 Equilibrium and stability

Clearly, ${q}={0}$ (i.e., $x=0$, and $\theta_1=\theta_2=0$) is an equilibrium point for the system (5) if and only if $f_1=f_2=0$ (see (4) ) at that position. In that case, since ${y}={0}$ if and only if ${q}={0}$, i.e., $x=0$, and $\sigma=0$, $\delta=0$ is an equilibrium configuration also for the system (7) or (10) .

Whenever ${F}^{(y)}+{\Phi}^{(y)}={0}$ (no escapement, no friction), the equilibrium point ${q}={0}$ (hence, ${y}={0}$) is Lyapunov stable by virtue of the Dirichlet's criterion, since it is an isolated minimum for the potential energy $V=\frac{1}{2}Kx^2-\sum\limits_{j=1}^2m_jg\ell_j \cos \theta_j$. Moreover, in presence of (2) or (3) , the equilibrium at the same position is asymptotically stable. Actually, we can write (2) as $D({q}){\dot {q}}$, with $D = \left( {\matrix{ { - \beta } \hfill & { - {\beta _1}{\ell _1}\cos {\theta _1}} \hfill & { - {\beta _2}{\ell _2}\cos {\theta _2}} \hfill \cr { - {\beta _1}{\ell _1}\cos {\theta _1}} \hfill & { - {\beta _1}\ell _1^2} \hfill & {8mm0} \hfill \cr { - {\beta _2}{\ell _2}\cos {\theta _2}} \hfill & {5mm0} \hfill & { - {\beta _2}\ell _2^2} \hfill \cr } } \right).$ Since $D$ is a negative-definite matrix, the energy balance $\frac{\rm d}{{\rm d}t}({\dot {q}}\cdot \nabla_{\dot {q}}{\cal L}-{\cal L})={\dot {q}}\cdot D{\dot {q}}< 0$ makes the energy a Lyapunov function tend to zero. In a similar way, one can proceed for the case (3) .

Nevertheless, if the escapement (4) is also operating, it must be said that this does not connote stability (not even equilibrium) of the system, whichever force ${F}^{({q})}$ we make use of.

3 Elimination of escapement, absence of damping and friction

In this paper, we are mainly involved in exploring the possibility of in-phase or antiphase synchronization settlements in absence of the escapement (4) . Now, we investigate the case ${F}^{({q})}={0}$ (i.e., no external devices are added to the system). Hence, the terms containing $f_1$, $f_2$, $f^{\pm}$ in (7) or (10) vanish. It must be said that such simplification entails a considerable advantage from the mathematical point of view, especially, if (4) is a step or discontinuous function, as it occurs in some mentioned models.

In our first investigation, even the friction forces are temporarily disregarded, the expected fact that the synchronization is unattainable will find confirmation.

Let us now examine the case, i.e., when $A_\beta^\pm$ and $B_\beta^\pm$ (see (9) ) are also removed, the equations of motion (7) are now simply $\frac{\rm d}{{\rm d}t}(\nabla_{\dot {y}}{\tilde {\cal L}})= \nabla_{y}{\tilde {\cal L}}$. Owing to the stability of the configuration ${y}={0}$, we are legitimated to replace (1) with the quadratic expansion $\frac{1}{2}({\dot {q}}\cdot {\overline A}{\dot {q}}-{q}\cdot {\overline V}{q})$, where ${\overline A}={\left(\begin{array}{ccc} m & m_1\ell_1 & m_2 \ell_2 \\ m_1 \ell_1 & m_1\ell_1^2 & 0 \\ m_2 \ell_2 & 0 & m_2 \ell_2^2\end{array}\right)}$ and ${\overline V}={\rm diag} (k, m_1\ell_1 g, m_2 \ell_2 g)$. In terms of ${\tilde {\cal L}}({y}, {\dot {y}})$, the linearized equations approximating (7) are

where

or explicitly

(13)

Clearly, the same set of equations can be obtained directly from (7) , by replacing (8) with the approximation (12) and by neglecting all second order terms. The fundamental frequencies (for both systems in ${q}$ and ${y}$) are found by solving det$(\lambda {\overline A}-{\overline V})=0$, leading to

(14)

with

(15)

If the two lengths $\ell_1=\ell_2=\ell$ are the same (but not necessarily the masses $m_1$ and $m_2$ are the same), $\Lambda=1$, $\rho=0$, and (14) can be reduced to

(16)

The quantity $\omega^2$ is an evident solution to it. This corresponds to the simplification of the third equation in (10) , i.e., $\ddot{\delta}= -g\delta/\ell$. The two remaining solutions are (say, $1$ with $-$ and $2$ with $+$)

(17)

We notice that $\omega_1\not = \omega_2$, since $(1+Y)^2>4Y (1-2\mu)$, with $\mu>0$. Moreover, both $\omega_1$ and $\omega_2$ are different from $\omega$, since $\mu< 1/2$, and from the relation, i.e.,

(18)

we deduce $\omega_1 < \omega$, and $\omega_2 >\omega$. Finally, we remark that

(19)

$\mu\approx 0$ means that the masses of the pendula are negligible with respect to the frame's one, however, when $\mu\approx 1/2$, the mass of the frame is negligible.

The generalized eigenvectors $(\lambda A_1 - V_1) {v}={0}$ correspond to $\lambda=\omega_j^2$, $j=1, 2$, and ${\overline \lambda}$ are, respectively, $(1, -2(\ell-g/\omega_j^2) ^{-1}, 0) ^{\rm T}$, $j=1, 2$, and $(0, (1-m_1/m_2) (1+m_1/m_2) ^{-1}, 1) ^{\rm T}$. Therefore, the solutions ${y}(t)$ starting from ${y}(0) =(x(0) , \sigma (0) , \delta(0) )^{\rm T}$, and ${\dot {y}}(0) =({\dot x}(0) , {\dot \sigma}(0) , {\dot \delta}(0) )^{\rm T}$ are

(20)

where $B$ is defined in Ref. [18], and

We write explicitly the solution in order to remark that

(i) The motion of $\delta$ is independent either of $x$, $\sigma$ and of $\mu$, even if $m_1\not =m_2$. As a consequence, none of the initial sets lead to the synchronization (i.e., $\delta(t)\rightarrow 0$), except for the matching of initial data ($\delta(0) =0$, and ${\dot \delta}(0) =0$). In that case, the pendula are exactly in-phase synchronized at any time.

(ii) The (not null) initial data $\delta(0) $ and ${\dot \delta}(0) $ produce an effect on the motion of $x$ and $\sigma$ if and only if the masses $m_1$ and $m_2$ are different.

The anti-phase synchronization ($\sigma(t)\rightarrow 0$) never occurs, except for the trivial case, i.e., $\sigma (0) =0$, $x(0) =0$, ${\dot x}(0) =0$, ${\dot \sigma}(0) =0$, and $m_1=m_2$ (whereas the date for $\delta$ are not null, otherwise, the system is at rest). For such data, $\sigma(t)\equiv 0$ (anti-phase at each time), and the frame is at rest.

(iii) The solution $x(t)$ is periodic if and only if $\frac{\sqrt{Y(1-2\mu)}}{1+Y}=\frac{q}{1+q^2}$, for some rational number $q\in (0, 1) $, if $m_1=m_2$, then $\sigma(t)$ is also periodic, otherwise, it must be added that the condition $\sqrt{\frac{Y}{1+Y}\frac{1+q^2}{q^2}}\in {\Bbb Q}$ in order to make $\sigma (t)$ periodic.

(iv) In order to have some understanding of the amplitudes of oscillation brought along each datum, we represent (18) by means of the parametric functions $B_\mu(Y)$, $Y\in (0, +\infty)$ with respect to the parameter $\mu\in (0, 1/2) $, and one can easily see that

when $0<\mu\leq 1/4$, the global maximum is $B_\mu(0) =\mu \ell$ (strictly decreasing), and when $1/4\leq \mu < 1/2$, the global maximum is $B_\mu(1-4\mu)=\sqrt{\frac{\mu}{2(1-2\mu)}}\frac{\ell}{2}$.

Consequently, the functions ${\mit \Phi}_\mu(Y)=\frac{2 B_\mu(Y)Y}{\mu\ell^2}$, ${\mit \Psi}^{(1) }_\mu(Y)=-\frac{2}{\ell}\frac{B_\mu(Y)\omega_1^2} {\omega_1^2 - \omega^2}$, and ${\mit \Psi}^{(2) }_\mu(Y)=\frac{2}{\ell}\frac{B_\mu(Y)\omega_2^2} {\omega_2^2 - \omega^2}$, where $\lambda_j(Y, \mu)$, $j=1, 2$, is (17) , verify

$0<\mu \le 1/4$, the global maximum is ${{\Psi }_{\mu }}(1/(1-4Y))=\frac{1}{\sqrt{\mu (1-2\mu )}}\frac{1}{2\ell }$;

$1/4\le \mu <1/2$, the supremum is $2/\ell $ (strictly increasing);

$\Psi _\mu ^{(1)}(0) = 0,\quad \Psi _\mu ^{(1)}(1) = 1/2,\quad \mathop {\lim }\limits_{Y \to + \infty } \Psi _\mu ^{(1)}(Y) = 1,\quad \Psi _\mu ^{(1)}$ is strictly increasing for any $\mu \in (0,1/2)$

As for the last point, the quantities $Y$ and $\mu$ are considered to be independent of each other by assuming, as an instance, the length $\ell$ and the mass of the frame $m$ as fixed and modifying $m_1$, $m_2$, and $k$.

Finally, we study the case of different lengths of the pendula. The circumstance is in the whole analogous to the previous case, except for the lesser simplicity of solving (14) , in order to achieve expressions similar to (20) . With respect to (15) , we have $\Lambda=1$ if and only if $\ell_1=\ell_2$, whereas $\rho=0$ whether $\ell_1=\ell_2$ or $m_1=m_2$. Hence, $\Lambda>1$ or $\rho\not =0$ causes somehow a "disturbance" to the exact solutions we found for $\ell_1=\ell_2$. In regards of that, we only point out that, letting $P(\lambda)$ be the third-degree polynomial of (14) , we have $P({\overline\lambda})=(1-\Lambda^2) (1-2\mu -Y)-2\Lambda^2\rho$. Let us have, for example, $m_1\leq m_2$, assuming $\ell_1<\ell_2$, $Y>1-2\mu$ (resp., $\ell_1>\ell_2$, $Y< 1-2\mu$), then, $P({\hat \lambda})=0$ for ${\hat \lambda}<{\overline\lambda}$ (resp., ~$>$) (see (15) ), i.e., the fundamental frequency decreases (resp., ~increases) by comparison with the case $\ell_1=\ell_2$. Similar remarks can be done about the other two solutions (17) , for which it is

(21)

4 Elimination of escapement, presence of damping and friction

Let us start from the system (10) . Defining $y={\dot x}$, $u={\dot \sigma}$, $v={\dot \delta}$ and setting the system as ${\dot {x}}={\cal F}({x})$ with ${x}=(x, \sigma, \delta, y, u, v)^{\rm T}$, we consider the linear approximation ${\dot {x}}=(J_{x}{\cal F}\vert_{{x}={0}}){x}$, where $J_{x}$ is the Jacobian matrix. Eliminating the terms containing the escapements $f_1$ and $f_2$, and reverting to the explicit expressions of the constant quantities (9) , one gets

(22)

Remark 1 In the presence of escapement (4) , the terms $\sum\limits_{j=1}^2\frac{1}{m_j\ell_j} (f_{j, \sigma}\sigma+f_{j, \delta}\delta+f_{j, u}u+f_{j, v}v)$ and $\frac{1}{m_1\ell_1}(f_{1, \sigma}\sigma+f_{1, \delta}\delta+f_{1, u}u+f_{1, v}v)- \frac{1}{m_2\ell_2}(f_{2, \sigma}\sigma+f_{2, \delta} \delta+f_{2, u}u+f_{2, v}v)$, where $f_{i, \zeta}=\frac{\partial f_i}{\partial \zeta}$ ($i=1, 2$, and $\zeta=\sigma, \delta, u$, and $v$) are calculated at the equilibrium ${x}={0}$, must be added to the fifth and sixth equations, respectively.

This time, the motion of $\delta$ (the last equation in (22) ) is not independent of the other variables if simply $\ell_1=\ell_2$. Actually, this equation is completely disentangled from the rest when also $m_1=m_2$, and $\beta_1=\beta_2$.

4.1 Localization of eigenvalues

The characteristic polynomial associated with the linear system (22) is

(23)

The terms with the odd exponent of $\lambda$ are due only to the friction contribution (2) . Actually, when $\beta_j=0$ for each $j=0, 1$, and 2, the characteristic problem (23) is equivalent to the one formulated in (14) .

Besides, the presence of a Lyapunov function which guarantees the asymptotical stability (see Subsection 2.1) , is worthwhile to assert the following

Property 1 The real part of each root of the polynomial (23) is negative.

Proof It is sufficient to make use of the Routh-Hurwitz criterion (see, for instance, Ref. ^[7]). writing (23) as $\sum\limits_{n=1}^6 a_n \lambda^n=0$, $a_6=1-2\mu$, $\cdots$, $a_0=\frac{g^2}{\ell_1\ell_2}\frac{k}{m}$ it can be checked, even though calculations last long, that the chain of seven numbers required for the mentioned criterion

where $b_1=a_4a_5-a_3a_6$, and $b_2=a_2a_5-a_1a_6$, consisting of all positive numbers. Since the sequence has no sign change, all the roots of the sixth degree polynomial (23) have negative real parts.

The criterion places the eigenvalues in the left half plane of the complex plane. In order to discriminate the occurrences of real roots, or conjugate pairs of complex roots of (23) one way could calculate the discriminant (see, for instance, Ref. ^[8]), which gives additional information on the nature of the roots, real or complex, although calculations for the sixth degree polynomial (23) are quite complex.

Our definitive aim is to infer some information about the qualitative behaviour of the system (22) , by means of locationing as much as possible the portion on the complex plane where the solutions of (23) lie.

Remark 2 It must be said that the Gershgorin circle theorem used to bound the spectrum is not$, $ in our case$, $ especially powerful$.$ It can be easily seen that the Gershgorin discs where the eigenvalues are confined do not keep them away from the origin of the complex plane$.$ This will be a crucial point in our analysis.

4.1.1 Case of identical pendula

(23) will be now discussed in the simpler case of identical pendula. We will assume

(24)

where the subscript p is necessary in order not to confuse with the quantities defined in (1) and (2) . Nevertheless, we presume that most of the shown results are still valid in the general case of different pendula.

We avoid to write again the system (22) in case of the assumption (24) , since the simplifications are evident. As we touch upon, the assumption (24) makes the motion of $\delta$ independent of the rest of the system. In fact, the most evident advantage of (24) is the reduction of the equation for $\delta$, i.e., $\ddot{\delta}+\dfrac{\beta_{\rm p}}{m_{\rm p}} {\dot \delta}+\dfrac{g}{\ell}\delta=0$ with negative eigenvalues, i.e.,

(25)

giving the solution

(26)

The characteristic equation (23) can be now written as

(27)

where the factorization is related to the uncoupling of $\delta$ from the system. We will focus our attention on the factor between square brackets, which gives the eigenvalues related to the motion of $x$ and $\sigma$.

If, in addition to $\eta$ (see (25) ), $\omega$ (see (16) ), and $Y=\frac{k}{m}\frac{\ell}{g}$ (see (15) ), we define

(28)

then the fourth degree polynomial in square brackets, (27) , can be written as

(29)

We notice that $\eta$, $X$, and $Y$ are dimensionless quantities (whereas the unit of $\omega$ is s$^{-1}$).

At this point, we use the Eneström-Kakeya (E-K) theorem^{[9, 11]}, in the following version:

Property 2 Let $p_n(\lambda)=a_0 + a_1 \lambda + \dots + a_{n-1}\lambda^{n-1}+a_n\lambda^n$ be a polynomial with $a_j>0$ for any $j=0, \cdots, n$}. Then, all the zeros of $p_n$ are contained in the annulus of the complex $z$-plane,

(30)

By writing (29) as $\sum\limits_{n=1}^4 a_n\lambda^n=0$ (even though the coefficients are different from the one used in the proof of Property 1) , the quantities needed for (30) are

(31)

It is immediate to check that $\frac{a_0}{a_1}\leq \frac{a_2}{a_3}$ and $\frac{a_1}{a_2}\leq \frac{a_3}{a_4}$ for any data $\eta$, $\omega$, $X$, and $Y$, hence,

(32)

and the positive quadrant ${\cal Q}=\{ (X, Y) | X>0, Y>0\}$ is splitted in the following four regions:

(33)

Sketching a physical reading, we see that a state $(X, Y)\in {\cal Q}$ in the right part of the quadrant ($X\gg Y$) exhibits a predominance of the damping force on the elastic force, both acted on the beam. On the contrary, in the left upper part of the quadrant ($Y\gg X$), the effects are reversed.

4.1.2 Locating synchronization regions

If one opts for fixing $\omega$ (length of pendula) and $\eta$ (damping of pendula), and let $X$ (friction of the beam), and $Y$ (elastic constant of the spring) vary, one can depict the four regions ${\cal Z}_i$ ($i=1, 2, 3$, and 4) on the quadrant ${\cal Q}$. In finding them, we see that the comparison of (31) leads to the following quadratic conditions in the variables $X$ and $Y$:

(34)

involving the construction of arcs of conics. It is easy to check that for $\eta^2>4/3$, the first three conditions define three regions on the positive quadrant $(X, Y)$ delimited by hyperbolae, while for $\eta^2< 4/3$, the first and the third conics are real ellipses ($\eta^2=4/3$, two intersecting lines). The fourth condition refers to a parabola attaining its vertex for some $X< 0$. The case $\eta\leq 1$, which will be examined deeper, is plotted in Fig. 2, where the curves are numbered in the same order with that in (34) .

We make use now of the localization carried out by the E-K Theorem in order to compare the eigenvalues governing $\delta$ (see (25) ) with those governing $\sigma$ (solutions to (27) ). We will hereafter focus on the case $\eta\leq 1$, which is physically more consistent, asserting that the complementary case can be conceptually treated in the same way.

Having in mind (25) , we are in the case of conjugate complex eigenvalues for $\delta$ with $-\frac{1}{2}\eta \omega$ as the real part and $\omega$ as the modulus. The main question is how the eigenvalues (25) are located with respect to the annulus (30) delimited by the radii (32) . We prove the following.

Proposition 1 For $\eta\leq 1$, the two roots (25) cannot lie in the half-plane Re $z<-\rho_{\rm M}$.

Proof The real part of the two roots (25) is $-\frac{1}{2}\eta \omega$. They belong to the half-plane Re $z<\rho_{\rm M}$ if and only if (see also (31) ) $\eta > 2\frac{\eta X+Y+1}{X+(1-2\mu)\eta}$ in ${\cal Z}_1 \cup {\cal Z}_3$, $\eta > 2\frac{X+(1-2\mu)\eta}{1-2\mu}$ in ${\cal Z}_2 \cup {\cal Z}_4$. However, the two conditions define empty regions in ${\cal Q}$, since they are equivalent to $Y+\frac{1}{2} \eta X < \frac{1}{2} (1-2\mu)\eta^2-1< 0$ and $X<-\frac{1}{2} (1-2\mu)\eta< 0$, respectively (we recall that $\mu < 1/2$, see (15) ).

Remark 3 The eigenvalues (25) belong to the semicircumference $|z|=\omega$, Re $z < 0$. The more specific question whether they lie in the semicircle $|z|\leq \rho_{\rm M}$, Re $z>0$ can be easily solved, by comparing $\omega$ with (32) , the second couple. (25) is within the semicircle if and only if $X>(1-2\mu)(1-\eta)$ when $\rho_{\rm M}= a_2/a_3$ (see (31) and if only if $Y>(1-\eta)X+(1-2\mu)\eta -1$ when $\rho_{\rm M}=a_3/a_4$. Graphically, each of the two regions ${\cal Z}_1\cup {\cal Z}_3$ and ${\cal Z}_2\cup {\cal Z}_4$ (see Fig. 4) is splitted by a straight line, and the required condition is true only on one side.

Our first conclusion is that the system cannot establish a status where the difference $\delta(t)$ decays to zero more rapidly than the sum $\sigma(t)$. Thus, the in-phase synchronization onset is inhibited. The result is consistent with the experimental detection, starting from Huygens, on the grounds that the anti-phase synchronization is indeed prevailing on the in-phase one in this sort of phenomenon.

We finally discuss the possibility for the system of establishing a status of anti-phase synchronization, still keeping $\eta\leq 1$. Making use once again of the localization (32) , we compare the real part of (25) with $\rho_{\rm m}$. Whenever $\eta \omega < 2\rho_{\rm m}$, the decay of $\sigma$ is expected to be faster than the one of $\delta$, so that the anti-phase synchronization is facilitated. Thus, we are going to check whether (see (31) )

(35a)

(35b)

The two conditions (35a) and (35b) define two half-planes above the straight lines $(2-\eta^2) X+\eta Y - \eta (1-4\mu)=0$ and $\eta X + (\eta^2-2) Y+2\mu \eta^2=0$, respectively.

As for the condition (35a), it can be easily seen as follows:

If $\mu \leq \frac{1}{2} \frac{2-\eta^2}{4-\eta^2}$, then (35a) is valid in ${\cal Z}_1 \cup {\cal Z}_2$ except for a triangular lower region cut off by the straight line (see Fig. 4). If $\frac{1}{2} \frac{2-\eta^2}{4-\eta^2}<\mu<\frac{1}{2}$, then (35a) is valid in all the set.

By recalling $\mu=m_{\rm p}/m$, we see that, $\eta$ being equal, the smaller is the mass of the pendula with respect to the mass of the system, the larger is the region which excludes the possibility of anti-phase synchronization, i.e., the removed zone.

The condition (35b) eliminates all the lower part from ${\cal Z}_3 \cup {\cal Z}_4$ and selects the region bounded from below by the straight line $\eta X + (\eta^2-2) Y+2\mu \eta^2=0$ and from above by the parabola related to the fourth condition in (34) . Moreover, if $\mu \geq \frac{1}{2} \frac{2-\eta^2}{4-\eta^2}$, then the selected region is bounded on the left hand side by a segment on the $Y$-axis. If $\mu <\frac{1}{2} \frac{2-\eta^2}{4-\eta^2}$, then the selected region is disjointed from the $Y$-axis. The different cases are summarized in Fig. 4.

Let ${\cal A}\subset {\cal Q}$ be the subset where (35a) and (35b) are both fulfilled. Whenever a state $(X, Y)\in {\cal A}$ is such that the four roots of (29) are all real, then, they exceed $\frac{1}{2} \eta\omega$, therefore, the decay of $\delta(t)$ is lengthier than the decay of $\sigma$, and the system shows a tendency to an anti-phase arrangement.

However, if some or all of the roots of (29) are not real, the condition $(X, Y)\in {\cal A}$ does not guarantee by itself that the real part of such solutions is greater than $\frac{1}{2} \eta \omega$. In order to overcome this problem, let us explain a method, without delving into the detail. Whenever the solutions to (29) are, as an instance, all complex, one can start from the conditions, i.e.,

(36)

where $a_2$, $a_3$, and $a_4$ are the same with those in (31) , and Re $\lambda$ is the real part of any root of (29) . The estimation (36) can be proved by setting the complex roots as $\xi_r\pm {\rm i} \kappa_r$ ($r=1, 2$), making use of the known relations between roots and coefficients, i.e., $\xi_1^2+\kappa_1^2+\xi_2^2+\kappa_2^2+ 4\xi_1 \xi_2 = a_2/a_4$, $2(\xi_1+\xi_2) =-a_3/a_4$ and finally recalling that $\xi_r^2+\kappa_r^2\geq \rho_{\rm m}^2$ ($r=1, 2$). In a mixed occurrence of two real roots and two complex conjugate roots, one argues in a similar way.

(36) can be used in order to separate the real parts away from zero. Actually, if we require $\frac{a_2}{a_4}\leq 2\rho_{\rm m}^2$, then Re $\lambda\leq B$, with $B$ being the negative solution of (36) , where "=" replaces "$\geq$". Finally, by demanding $-\frac{1}{2}\eta \omega\geq B$, we achieve that the real parts of the solutions to (29) exceed $\frac{1}{2} \eta \omega$. The conditions $\frac{a_2}{a_4}\leq 2\rho_{\rm m}^2$ and $-\frac{1}{2}\eta \omega\geq B$, once they have been expressed in terms of $X$ and $Y$ via (31) and (32) , will determine the regions of ${\cal Q}$ where the anti-phase synchronization is facilitated.

5 Conclusions

In the frame of the study of coupled oscillations of two pendula, we intend to pursue two objectives as follows:

(a) To formulate the model in the experimental context as general as possible.

(b) To develop the corresponding mathematical problem in a simplified case, drawing the attention to the fact that some classical results about the localization of the spectrum of a matrix can allow us to predict the qualitative behaviour of the system.

The results achieved in Section 3 confirm the expectation that the synchronization is not possible whenever friction and dissipative forces are absent, even if the features of the pendula (length, mass) are different. In Section 4, the mathematical analysis is consistent with the observed tendency of the system, showing a greater aptitude to the anti-phase rather than the in-phase synchronization (see Fig. 3). The main theoretical results are summarized by plotting the regions where synchronization is facilitated, as shown in Fig. 4.

We select the parameters $X$ (see (28) ) and $Y$ (see (15) ) in order to plot on the quadrant ${\cal Q}$ a certain number of regions in each of which the system develops in a different way. The parameters $\mu$ (defined in (15) ) and $\eta$ (see (25) ) are considered to be constant, but different choices can be made in order to represent the states (for instance, fixing the friction of the $\beta_0$ but varying the damping of the pendula $\beta_{\rm p}$).

The sections of the complex plane where the spectrum is confined and the regions on ${\cal Q}$ are not computed via the numerical simulation, however, they are predicted by the analysis of the spectrum of the linearized system (22) .

Within specific ranges of the parameters, the tendency of the system to evolve towards the anti-phase synchronization, rather than the in-phase one, is predicted by the present analysis. This conclusion is in step with the real development of the phenomenon, starting from the Huygens' observation in the seventeenth century of the 180 out of phase swings (see Ref. ^[5])).

Generally speaking, our purpose is to highlight that the method, beyond the specific circumstance which has been exerted to, can be extended to more general situations where the roles of certain parameters are exchanged or some restrictions are relaxed. At the same time, the analysis is appropriated, in our mind, to be combined with a simple numerical approach. Even if we restrict our investigation to a pure mathematical problem, the explicit results depicted as regions on the plane of the physical parameters facilitate a possible validation via the computer.

Both the mentioned points (generalization and matching via computer) are now topics for our current research. More precisely, the in-depth investigations of the problem are as follows:

(i) Extending the method to the case of different pendula, by examining (23) via the spectrum analysis performed in the simpler case, or by elaborating a formula similar to (21) for evaluating the effects of a variance in the features of the two pendula.

(ii) Estimating the time of decay of the motion and discard the situations where a very short time from the starting time to the almost rest state would not produce interesting cases.

(iii) Verifying the analytical outcome by means of simulations via the computer, either for calculating the spectrum and for tracing the profiles of $\delta$ and $\sigma$.

(iv) Locating the eigenvalues in more restricted regions, by making use of some generalization of the E-K Theorem (see Ref. ^[10]), which confines the spectrum in specific circular sectors of the complex plane.

(v) Checking whether the decrease of $\beta_i$ ($i=1, 2, 3$) to zero in (22) will lead to the solutions described in Section 3.

(vi) Adding the effects (4) of an escapement.

As for the point (i), the difficulty comes from the non-possibility of factorizing the characteristic equation (23) as in (27) , so that the eigenvalues for $\delta$ cannot longer be separated from the rest of the spectrum. At this point, the small coefficients ($\ell_1-\ell_2$, etc.), which join the equation for $\delta$ (the last equation in (22) ) with $x$ and $\sigma$, will play a significant role.

Moreover, the point (vi) renders the analytical problem much complex and deeply different from the present one. As we remark, the new formulation requires a non-trivial discussion of equilibrium and stability and of the existence and the regularity of solutions, where the difficulty arises from the typical (but experimentally adequate) discontinuous profile of (4) .