Dynamics of dipole in a stationary non-homogeneous electromagnetic field

The non-relativistic equations of motion for a dipole in a stationary non-homogeneous electromagnetic field are derived and analysed. It is shown that they are Hamiltonian with respect to a certain degenerated Poisson structure. Described by them dynamics is complex because the motion of the centre of mass of the dipole is coupled with its rotational motion. The problem of the existence of linear in momenta first integrals which can be useful for the separation of rotational motion is discussed. The presence of such first integral appears to be related with a linear symmetry of electric and magnetic fields. Also results of search of quadratic in momenta first integrals for uniform and stationary electromagnetic fields are reported. Deriving equations of motion of a dipole in arbitrary stationary electromagnetic fields and analysis of described by them dynamics is important for the construction of electromagnetic traps for polar particles.


Dynamics of dipole in a stationary non-homogeneous electromagnetic field
Maria Przybylska 1* & Andrzej J. Maciejewski 2 The non-relativistic equations of motion for a dipole in a stationary non-homogeneous electromagnetic field are derived and analysed. It is shown that they are Hamiltonian with respect to a certain degenerated Poisson structure. Described by them dynamics is complex because the motion of the centre of mass of the dipole is coupled with its rotational motion. The problem of the existence of linear in momenta first integrals which can be useful for the separation of rotational motion is discussed. The presence of such first integral appears to be related with a linear symmetry of electric and magnetic fields. Also results of search of quadratic in momenta first integrals for uniform and stationary electromagnetic fields are reported. Deriving equations of motion of a dipole in arbitrary stationary electromagnetic fields and analysis of described by them dynamics is important for the construction of electromagnetic traps for polar particles.
The main aim of this paper is to give a general framework for a study dynamics of electric dipole in a stationary non-homogeneous electromagnetic field. In our considerations a dipole is modelled by two equal, and opposite sign point charges. The distance between them is fixed, and the masses of the charges are not necessarily equal. From the dynamical point of view, such a dipole is a degenerated charged rigid body-one of its principal moments of inertia vanishes. Our goal is to derive and investigate equations of motion in the most general case when the mass centre motion of the dipole is coupled with its rotation.
Motivations for this investigation come from our study of electromagnetic traps for charged and neutral particles [1][2][3][4][5] . More specifically, in 5 we investigated the motion of a small dipole in a special configuration of nonhomogeneous electromagnetic field that was proposed as a model of a trap for polar particles. These investigations showed that for a study similar traps we have to construct a general description of a dipole in an electromagnetic field of arbitrary form.
We showed in 4 that even in a stationary homogeneous electromagnetic field the motion of a dipole is complex and, in general, it is chaotic. If the field is not homogeneous, then having in mind possible physical applications, it makes sense to assume that the dipole is small. But then it is not clear what are the properties of reasonable approximations of the equations. In the commonly used approximation a dipole is just a neutral point mass which carry a dipole moment. However considering polar particles we have to take into account that they have finite dimensions and that their mass distribution in not spherically symmetric. Moreover, the charge and the mass distributions of a polar particle do not need to be necessarily proportional to each other. This is why in our considerations we assume that the dipole is small but we take into account that it can be inertially non-symmetric. That is, its mass centre does not necessarily coincide with its centre of the charge distribution. This asymmetry is described by the parameter δ = l(m 2 − m 1 )/(m 1 + m 2 ) , where m 1 and m 2 are the masses of the charges and l is the charges separation distance.
The Hamiltonian formalism gives appropriate tools and language for classical description of motion. However, considering problems where electromagnetic is involved the canonical description can be done in two ways, see 4 . For a single charge q of mass m in an electromagnetic field the canonical momentum conjugated to its position r is P = mṙ + qA(r) , where A(r) is a vector potential so that B(r) = ∇ × A(r) is the magnetic field. However, it is convenient to use variables r and linear momentum p = mṙ . Of course, in these variables the equations of motion are still Hamiltonian but with respect to a non-canonical symplectic structure which depends on the magnetic field, see, for example 4,6 . The advantage is that the system is gauge invariant. The equations of motion for a dipole can be obtained by a restriction of non-canonical Hamilton's equations for two charges. The constrain, the constant distance between charges, is holonomic so one can expect that the constrained system is Hamiltonian, and, it is. The other way is just to start from the Lagrangian formulation and use local generalised coordinates. However, the first approach needs not so well known theoretical tools, and the second gives highly complicated formulae. This is why we choose simple and direct method to derive the equations of motion for a dipole. We just use the Newton-Lorenz equations for a charge in an electromagnetic field and the basic laws of mechanics. We demonstrated it for a dipole in a stationary homogeneous electromagnetic field in 4 . In this article we show that this method works in a general case. The derived equations are Hamiltonian with respect to a degenerated Poisson bracket and they are starting point for our further analysis.
The key goal of out analysis was derivation of equations of motion for a small dipole. To this end we expand into the power series the right hand sides of general equations and truncate them at the first term proportional to the asymmetry parameter. We proved that these equations are Hamiltonian with respect to an appropriate Poisson bracket. This fact is not evident.
As we already mentioned the motion of the centre of mass of the dipole is coupled with its rotational motion. The centre of mass does not move with constant velocity. Nevertheless, if the electromagnetic field is stationary and homogeneous, then existence of first integrals linear in the momenta allows to separate equations for the rotational motion of the dipole, see 4 for details. This is why it is important to know whether the obtained equations of motion of the dipole admit first integrals linear in momenta. We analyse this problem in details and show among other things that if such a first integral exits, then the vector fields E(r) and B(r) have a certain linear symmetry field S(r) , that is they both commute with S(r) : [E(r), S(r)] = [B(r), S(r)] = 0 , where [·, ·] means the Lie bracket of vector fields.
The plan of this paper is as follows. In Section "Equations of motion" we derive equations of motion of an electric dipole in arbitrary time-independent electromagnetic fields using the Newton-Lorentz equations of translational and rotational motion. Motion of the dipole is governed by four vector equations which describe the motion of its centre of mass and rotations of the dipole depending on values of electric and magnetic fields at positions of two charges creating the dipole. These equations are Hamiltonian with respect to a certain degenerated Poisson structure with two Casimir functions. In Section "Equations of motion of a small dipole", we derive truncated equations of motion of the dipole under assumptions that its dimension is small with respect to the size of non-linearities of the electromagnetic field and in expansions terms up to the second order with respect to length of the dipole are preserved. It is interesting that these equations are also Hamiltonian with the Poisson structure that is the truncation of the Poisson structure from the previous section and that has the same two Casimir functions. Section "First integrals" contains results of the research of first integrals using the direct method. The most general results are obtained for the case when searched first integral is linear in linear momentum of the centre of mass of the dipole, coordinates of the dipole vector and angular velocity of rotational motion of the dipole. The forms of such first integrals and expressions for the electromagnetic field admitting such first integrals are given. For first integrals quadratic in the mentioned variables we were able to finish analysis only for stationary homogeneous electromagnetic fields and these results are given in Section "Final notes and comments". Also in this section we discuss limits of applicability of the presented model of an electric dipole. In Appendices 1 and 2, Lagrange equations for the cases of the dipole in electromagnetic fields without the smallness of the dipole assumption and with this assumption are given, respectively. We underline that in our considerations we do not assume that masses of two opposite point charges are the same as it is usually made.

Equations of motion
In this section, we derive equations of non-relativistic motion of an electric dipole in a time independent electromagnetic field without assumption about the smallness of the dipole with respect to the size of non-linearities of the electromagnetic field. We show that it is possible to present them as Hamilton's equations with respect to a certain degenerated Poisson structure.
At first let us fix notation. A vector x ∈ R n will be considered as a one column matrix x = [x 1 , . . . , x n ] T . Thus, the scalar product x · y of two vectors is x · y = x T y = y T x . The gradient ∇f (x) of a scalar function f (x) is considered as a vector, and we denote f ′ (x) := df (x) := (∇f (x)) T . With this convention f ′ (x)a = a · ∇f (x) for a vector a ∈ R n . For a vector function we denote and The Hessian f ′′ (x) of a scalar function f (x) is For a vector function (1) we set For x, y ∈ R 3 the scalar and vector products we denote by x · y and x × y , respectively. For a vector x ∈ R 3 , the corresponding skew-symmetric matrix is The map x → x is the standard isomorphism between Lie algebras so(3, R) and R 3 with the vector product as algebra multiplication. It has the following properties where Id 3 is the three-dimensional identity matrix. Moreover, for an arbitrary 3 × 3 matrix M ∈ M(3, R) the following identity holds true. Let us recall that motion of a point with mass m, electric charge q and radius vector r in a stationary electromagnetic field is described by the Newton-Lorentz equation where E(r) and B(r) denote the electrostatic and magnetostatic fields, respectively.
We consider a dipole which is modelled by two charges +q and −q with the respective masses m 1 and m 2 . Their positions in an inertial frame are given by radius vectors r 1 and r 2 , see Fig. 1. The centre of mass of the dipole is localised at www.nature.com/scientificreports/ Here r 1 − r 2 = le , where e is the unity vector l = |r 1 − r 2 | is the distance between charges. Thus, the configuration space of the dipole is R 3 × S 2 , where S 2 ⊂ R 3 is the unite sphere. The force acting on the centre of mass of the dipole reads so its motion is described by Newton's equation Notice that the right hand side of this equation depends only on (r, e) and (ṙ,ė).
The total angular momentum of the dipole is where we denoted Thus e · ω = 0 , so we can write ė = ω × e. The torque of electrostatic forces acting on the dipole is Taking into account identities (12) the Newton equation (14)  In Appendix 1, we show that this system can be derived from appropriate Lagrange equations. System (23) has the energy integral (12) (15) L = m 1 r 1 ×ṙ 1 + m 2 r 2 ×ṙ 2 = mr ×ṙ + m 1 m 2 m l 2 e ×ė = mr ×ṙ + Iω, (20) L = mr ×r + Iω. Proof Facts that the bracket is bilinear, antisymmetric and satisfies the Leibniz identity follow directly from its definition. It is also simply to show that (x)H ′ 0 (x) = 0 and (x)H ′ 1 (x) = 0 . However checking the Jacobi identity is not trivial. Let us denote expressions for i, j, k ∈ {1, . . . , 12} . Taking into account the antisymmetry of the bracket it is enough to show that J i,j,k = 0 for 1 ≤ i < j < k ≤ 12 . From the form of matrix (x) it is easy to deduce that {r i , x j } = δ j,3+i for i = 1, 2, 3 and j ∈ {1, . . . , 12} . Thus, J i,j,k = 0 for i < 4 . Moreover, {e i , x j } = 0 for j < 10 , and {e i , g j } = ε i,j,k e k , so J i,j,k = 0 for i = 7, 8, 9. Direct calculations that for remaining values of indices the are 16 cases for which J i,j,k = 0 , however for these cases J i,j,k is a linear homogeneous polynomial of ∇ · B(r 1 ) and ∇ · B(r 2 ) . For example one can obtain Hence, if ∇ · B(r) = 0 , then J i,j,k = 0 for 1 ≤ i < j < k ≤ 12 , and this finishes our proof.

Equations of motion of a small dipole
If the dimension of the dipole is small with respect to the size of non-linearities of the fields, then we can obtain approximate equations of motion just expanding the force and the torque acting on the dipole into power series and truncate these expansions at appropriate order. In our considerations we need expansions up to the second order with respect to l. For calculations it is convenient to introduce the following notation Then and so we get . In an electric field this force is non-zero only in an inhomogeneous electric field and we recognise famous expression for this force qlE ′ (r)e = (d · ∇)E(r) , see, for example, Section 4.1.3 in 7 . Next two terms are related with the action of a magnetic field on electric dipole and are not such frequent in literature but one can find them for example in 8 , see also additional explanations in Section 2.1 of 9 . Expressions for forces in the first row (39) agree with these in Eq. (33) in 8 remembering that moving electric dipole in a magnetic field also has magnetic moment µ = d ×ṙ . Expressions for forces given in the second line are present only for asymmetric dipoles with m 2 = m 1 and according to our knowledge are new. The total angular momentum of the dipole is given in (15) and the expansion of expression (17) for the torque to the second order with respect to l gives Here we also recognise two well known terms of torque for a dipole in an electric field de × E(r) = d × E(r) and in a non-homogeneous field dr × E ′ (r)e = r × (d · ∇)E(r) , see, for example, Section 4.1.3 in 7 , less known terms related to the presence of magnetic field, see 8 and new terms related to the asymmetry of the dipole given in the second and in the third lines of Eq. (40).
Taking into account equalities (39) and (20) the truncated Newton equation for the rotational motion (21) reads Using definitions of the linear p = mṙ and the angular momentum g = Iω of the dipole, the Newton equations can be rewritten as the system of first order differential equations In Appendix 2, we show that this system can be derived from an appropriate approximated Lagrange function given in (112). General properties of equations (42) in the special case of uniform and stationary electromagnetic (40) www.nature.com/scientificreports/ fields were analysed in 4 . In this case these equations can be reduced to a Hamiltonian system with two degrees of freedom which is non-integrable except two cases with specific values of parameters. Also results of numerical simulations for a water molecule and CdS nanocrystal physical dipoles were shown. Application of equations (42) to the special case of non-constant electromagnetic fields is given in our previous article 5 . Namely, in this paper we analysed numerically the dynamics of a small dipole in the superposition of quadrupolar magnetic field B(r) = B 1 2 (−x, −y, 2z) and electric field E(r) = E z e z − ∇� 3 r) , with sextupole potential , where U 0 is the trapping voltage and D is the characteristic trap size. We demonstrated that such a non-homogeneous electromagnetic field seems to be useful for trapping of polar particles. We shown this by numerical calculations of trajectories of the centre of mass various physical polar particles: the DAST nanocrystal, the cellulose nanocrystal and protein lac repressor. Simulations with different initial conditions give different mass centre trajectories but they remain bounded to a restricted volume. However the general properties of equations (42) have not been analyzed so far, and filling this gap is the goal of this article. We hope that this analysis enables a better understanding of the complex dynamics of the dipole in non-homogeneous electromagnetic fields and it will help for example in the construction of new electromagnetic traps for polar particles.
System (42) has two geometric first integrals and the energy integral Checking that H is a first integral needs some efforts. Namely, during calculations we have to use the following identities In the last formula X = B or X = B ′ (r)e . All of them can be checked directly and we left it to the reader. Let F and G be two smooth functions defined on M 12 = (R 3 ) 4 = R 12 with coordinates x = (r, p, e, g) . We define their bracket as with We have the following. Proof Except the Jacobi identity all statements of the lemma are easy to check. To prove the Jacobi identity we proceed like in Lemma 2.1. However, we have to use the fact that from Gauss's law for magnetism ∇ · B = 0 it follows that all first and second partial derivatives of the left hand side of this condition also vanish identically. Now, we can state the main result in this section. A direct proof of this theorem we left to the reader. We showed that system (42) restricted to level (47) M = (r, p, e, g) ∈ M 12 | e · e = 1, e · g = 0 www.nature.com/scientificreports/ is Hamiltonian with five degrees of freedom. It has one first integral H. Hence, for its integrability four additional commuting first integral are needed.

First integrals
It seems that for general form of E(r) and B(r) the only first integrals of system (42) are H, H 0 and H 1 . In this section, we look for cases when the system admits an additional first integral which is polynomial of low degree with respect to momenta p and g . However, even for the case of first integrals which are linear with respect to p and g , without additional restrictions concerning coefficients of such a polynomial the problem is not tractable. Thus, we limited ourselves to searching of first integral F linear in variables p , g , and e with coefficients which are functions of r , that is Equating to zero Lie derivative we obtain a system of 62 partial differential equations for these unknown coefficients P(r) , Q(r) and R(r) , as well as for the electric and magnetic fields E(r) and B(r) . Fortunately partial differential equations for coefficients P(r) and R(r) separate from remaining equations. Solving them we obtain that R(r) = R is a constant vector, and Additionally, the following equation  For further consideration we assume that B(r) = 0 . So, by the above lemma, if the first integral (48) exists than it has the form As we will see it later vector S(r) := C + R × r plays an important role.
The Lie derivative of F is a polynomial of third degree in variables (p, e, g) of the form given by (52) where polynomials R ijk have the forms (48) F = P(r) · p + Q(r) · e + R(r) · g.

Lemma 4.2 If (55) is a first integral of system (42), then
Proof We deduce these two necessary conditions from equations R 010 = 0 and R 110 = 0 . Term R 010 , see (56a), can be rewritten as where we used symmetry of E ′ (r) . As R 010 = 0 for an arbitrary e we get first condition in (57). Now we rewrite four terms of R 110 , see (56b) in the following way Using property (8), and facts that TrB ′ (r) = ∇ · B(r) = 0 and Tr R = 0 , we rewrite this matrix in the form As R 110 = −p · Ue = 0 for arbitrary p and arbitrary e , the antisymmetric matrix U vanishes, and thus, also the corresponding vector This ends the proof.
For further considerations it is important to notice the following. The commutator [X(r), Y (r)] of arbitrary vector fields X(r), Y (r) in R 3 means the Lie bracket of vector fields defined as where X ′ (r) and Y ′ (r) mean Jacobian matrices ∂ j X i and ∂ j Y i for i, j = 1, 2, 3.
Proof When we substitute in the above definition X(r) = S(r) and notice that X ′ (r) = S ′ (r) = R , then we can rewrite this expression as and this ends the proof.
In other words, if the system admits a first integral (55), then S(r) is a symmetry field of electric and magnetic fields E(r) and B(r) . This property gives strong restrictions on these fields. Moreover, if δ = 0 , then all the conditions are fulfilled. So, we can summarize the above considerations with the following. For δ = 0 we have three additional conditions R 020 = 0 , R 120 = 0 and R 021 = 0 . The term R 020 can be rewritten as (56e) R 021 = S(r) · (g × e) × B ′ (r)e + (e · B(r)) R · (g × e) .  www.nature.com/scientificreports/ where Condition R 020 = 0 implies that the symmetric part of matrix G vanishes and this is equivalent to the following equality Polynomials R 021 and R 120 can be rewritten as Thus, conditions R 021 = 0 and R 120 = 0 and are equivalent to respectively. Each of these conditions defines three quadratic forms in components of e , and all of them have to vanish. So, altogether we have 36 additional conditions when δ = 0. Let us analyse the above conditions for given vectors C and R . They have form of a system of partial differential equations which must be fulfilled by fields E(r) and B(r) . We consider three complementary cases.
L1 Assume that R = 0 and C = 0 . By a proper orientation of the inertial frame we can achieve that C = (0, 0, C 3 ) , then by a rescaling we obtain that C 3 = 1 . In effect in this case we will assume that C = (0, 0, 1). L2 Assume that R = 0 and C = 0 . By the same arguments as in the previous case we will assume that R = (0, 0, 1). L3 Assume that R = 0 and C = 0 . By a proper orientation of the inertial frame we can achieve that R = (0, 0, 1) , and C = (0, C 2 , C 3 ) . A translation r → r + a transforms C → C + R × a . Taking a = (0, −C 2 , 0) we obtain C = (0, 0, C 3 ) . So, in this case we will assume that R = (0, 0, 1) and C = (0, 0, c) , c ∈ R.
In each case we perform analysis into two steps. First, we give necessary conditions which guarantee that in (52) terms R 010 and R 110 vanish. That is, by Lemma 4.2, we find E(r) and B(r) which satisfy (57). If δ = 0 , then these conditions are also sufficient for the existence of first integral (55). In the second step we assume that E(r) and B(r) satisfy derived conditions and we find those which satisfy R 020 = 0 , R 210 = 0 and R 012 = 0 , or, equivalently, which satisfy conditions (62), (63a) and (63b), respectively.     www.nature.com/scientificreports/ Proof We can assume that fields E(r) and B(r) have the form given by the previous lemma. For R = 0 and C = (0, 0, 1) condition (62) reads E ′′ 3 (r) = 0 , but it is fulfilled because by the remark before this lemma E 3 (r) is a constant function.
For the above choice of C and R condition R 021 = 0 , see (63a), reads and it implies that B 1 (r) = B 1 and B 2 (r) = B 2 are constant functions. Condition R 120 = 0 , see (63b), has the form and it does not give new restrictions on B(r) because we already proved that B 1 (r) and B 2 (r) are constant functions. This ends the proof.
As it is well known, if a Hamiltonian system generated by Hamiltonian H has a first integral F then Hamiltonian vector fields X H and X F commute. For the system considered in this paper the Poisson structure is not linear, see (46). Thus, even a very simple first integral F can generate complicated non-linear vector field. Amazingly, the first integral F given in (55) generates Hamiltonian vector field X F which is affine with respect to coordinates (r, p, e, g). F given by (55) is a first integral of the system (42), then the Hamiltonian vector field X F has the form Proof The partial derivatives of function F with respect all variables are following From definition the vector field generated by F is X F = P(x)∇F(x) . As we also split components of vector X F accordingly, that is we set X F = (X F,r , X F,p , X F,e , X F,g ) . Using explicit form of J(x) , see (46), we easily obtain It is more complicated to calculate X F,r . We have Now we regroup and simplify several terms. Notice that where the last equality follows from fact that the cross product satisfies the Jacobi identity. Moreover Finally, applying identity (8) with M = B ′ (r) and taking into account that ∇ · B(r) = TrB ′ (r) = 0 we get With the above simplification we obtain B ρ (ρ, ζ ) = 0 and B z (ρ, ζ ) = c ρ B ϕ (ρ, ζ ).

Lemma 4.13 If
(86) (87) (89) www.nature.com/scientificreports/ The last equality occurs because terms in square brackets in the above formula vanish since F is by assumption first integral, so conditions (57) and (63b) are fulfilled. We calculate the other half of components the last equality occurs because the term in the square backed vanishes according to condition (63a).

Final notes and comments
We showed that equations of motion of a small dipole are Hamiltonian with respect to a certain Poisson structure. It is not trivial as a truncation of Hamilton's equation does not give Hamiltonian equations in general. Fact that the existence of a first integral linear in momenta is connected with certain symmetry of the system is well known. For a single charged particle motion in stationary electromagnetic fields the results of group symmetry analysis are known, see, for example 12,13 . However system considered in this paper is much more complicated -it has five degrees of freedom and the form of electromagnetic field is not specified. This is why symmetry analysis of the system can be perform only in restricted cases. We show that if equations of motion for a small dipole admit a first integral linear with respect variables (p, g, e) , then they possess a linear non-homogeneous Hamiltonian symmetry field. The proof of Lemma 4.13 shows that this fact is not obvious. This symmetry condition restricts strongly fields E(r) and B(r) -they have to commute with linear vector field S(r) = C + R × r . Thanks to this we were able to classify electromagnetic fields which admit the prescribed first integral.
As we already mentioned, it is hopeless to investigate integrability of the system without explicit specification of fields E(r) and B(r) . As we shown in 4 , even if we assume that these fields are constant the system is not integrable in general. It appears that the relative orientation of these fields is important. As integrability analysis presented in 4 concerns the reduced form of equations (42) here we give several remarks concerning the integrability of unreduced equations of motion.
For constant E and B first integrals linear in p , e and g always exit. It is enough to take R = 0 and arbitrary C in formula (55), then all conditions R i,j,k = 0 with R i,j,k given by (56a)-(56e) are fulfilled. Thus we have three first integrals given by components of Integrals H, H 0 , H 1 and Q are functionally independent and pairwise commute but for the integrability of the system still one first integral is missing. In our previous work 4 we used integrals Q for the reduction of the equations of motion of the dipole to a Hamiltonian system with two degrees of freedom.
If δ = 0 , then for constant electric and magnetic fields conditions R i,j,k = 0 with R i,j,k given by (56a)-(56e) reduce to two vector conditions R × E = 0 and R × B = 0 . In the case when vectors R , B and E are parallel, substitution C = 0 in formula (55) (94) S = (B · Q) B · (g + r × Q) + mr · (Q × (E × B)). In the case when B = 0 and E = [E 1 , E 2 , E 3 ] T = const first integrals Q simplify just to linear momenta of the centre of mass Q = p , moreover also angular momentum K = r × p and function G = E · g are preserved during evolution of the system. Among them functionally independent is the set of first integrals H 0 , H 1 , H, p 1 , p 2 , p 3 , K 1 , K 2 , G . This system is integrable because first integrals H 0 , H 1 , H, p 1 , p 2 , p 3 , G commute each other. In fact it is maximally super-integrable with two additional first integrals K 1 , K 2 . Non-vanishing Poisson brackets of first integrals are following In this case translational motion of the centre of mass separates from rotational motion and the centre of mass moves with constant momentum p(t) = p 0 along straight line r(t) = p 0 m t + r 0 . Equations of motion for rotational motion simplify to This system has three obvious first integrals It has an additional first integral so it is integrable on the common level of geometric integrals H 0 = 1 and H 1 = 0 . To see this explicitly let us parametrise vector e by spherical coordinates

{S, Q}
The canonical momenta conjugated with (q 1 , q 2 ) are p 1 = Iq 1 and p 1 = I sin(q 1 ) 2q 2 . By a proper choice of the inertial frame we achieve that E = (0, 0, E) . In these canonical variables the Hamiltonian function reads Clearly q 2 is a cyclic variable so p 2 is a first integral. Fixing its value we obtain a system with one degree of freedom. It can be solved explicitly in terms of elliptic function. In fact, the Hamiltonian (99) describes a spherical pendulum for which explicit solutions can be found for example in par. 55 of 14 .
One can continue search of first integrals of system (42) which are polynomial of degree two with respect to variables p , g and e . Vanishing of time derivative of such first integral gives a non-homogeneous polynomial of degree four of variables (p, e, g) . Analysis of obtained 22 conditions which are homogeneous polynomials of degree 1, 2, 3 or 4 in these variables has appeared to be much more complicated because system (42) possesses already certain quadratic first integrals: Hamiltonian H and Casimir functions H 0 and H 1 and we were not able to complete it.
Considered equations of motion were obtained in non-relativistic framework under assumption that velocities of material points are very small in comparison with the velocity of light.
We did not consider the energy losses resulting from the fact that a charge moving with acceleration emits electromagnetic waves which carries the energy, the momentum and the angular momentum. The change of energy of an accelerated point charge in time is described by means of famous Larmor's formula where µ 0 is the vacuum permeability, q is the charge of the particle and a is its acceleration, see, for example, Sec 11.2.1 in 7 . Since the dipole is created by two charges of opposite signs which radiate independently its energy loss is To neglect the radiation phenomenon of an accelerating dipole, the translational acceleration of its centre of mass and the angular velocity related to the rotation of the dipole around its centre of mass must be sufficiently small.
Loss of the energy, the momentum and the angular momentum means damping of motion of the point charge. Using conservation of energy and Larmor's formula (100) one can obtain the so-called Abraham-Lorentz formula for the radiation damping force where ε 0 is the vacuum permittivity, see, for example, Sec. 11.2.2 in 7 , Sec. 75 in 15 or 16 , which we should add to the right-hand side of the Newton-Lorentz equation (9) if we want to take into account radiation phenomena.
(102) www.nature.com/scientificreports/ Although that the Abraham-Lorentz formula is commonly accepted non-relativistic expression for the radiation damping force acting on a point charge in an external electromagnetic field it suffers from certain nonphysical solutions: runaway solutions and preacceleration, for more details see, for example, Sec. 11.2.2 in 7 or 16 . One of the most frequently applied remedies is to use the Landau-Lifshitz formula as an alternative to the Abraham-Lorentz expression. The Landau-Lifshitz formula is obtained by expressing the time derivative of acceleration in (102) by means of the acting external forces (neglecting the radiation damping force), thus by electromagnetic field and its derivatives, for details see Sec. 75 in 15 . These difficulties also persist in the relativistic version of the Abraham-Lorentz formula called the Lorentz-Abraham-Dirac formula and there is a relativistic version of the Landau-Lifshitz formula described for example in Sec. 76 of 15 .

Appendix 1: Lagrange equations for two opposite sign point charges with a fixed distance between them
The Lagrange function of two point masses m 1 and m 2 with charges +q and −q localised at r 1 and r 2 ,respectively, has the form where �(r) and A(r) are scalar and vector potentials of the electromagnetic field defined by the formulae E(r) = −∇�(r) and B(r) = ∇ × A(r).
We assume that the distance l between charges r 1 − r 2 = le is fixed thus we have to take into account the presence of the constrain e 2 = 1 . Moreover, from the above Lagrangian we subtract the total time derivative Thus as result we consider the following Lagrange function Taking into account definitions of the centre of mass radius vector r and the dipole vector e given in (12) expression for L can be rewritten as where we also used identity m 1 (δ + l) + m 2 (δ − l) = 0 . To obtain the first Lagrange equation we need to calculate derivatives

As result we obtain
We used � ′ (r i ) = −E(r i ) and identities for i = 1, 2.
(105) where we used the well known formula for the vector triple product. When we recall definitions of linear and angular momenta, then we immediately obtain system (23) from the above Lagrange equations (107) and (110).

Appendix 2: Lagrange equations for a small dipole
The length of dipole l is small, so it is justifiable to expand the scalar and vector potentials into power series with respect to l up to a certain order. In further considerations it is necessary to keep terms up to the second order.
r · A ′ (r)e − e · A ′ (r)ṙ = −ṙ · (e × B(r)), e · A ′ (r)e − e · A ′ (r)ė = −ė · (e × B(r)), r · A ′′ (r)(e, e) − e · A ′′ (r)(ṙ, e) = −ṙ e × (B ′ (r)e) .  www.nature.com/scientificreports/ 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:// creat iveco mmons. org/ licen ses/ by/4. 0/.