Basins of Attraction in the Copenhagen Problem Where the Primaries Are Magnetic Dipoles

We deal with the Copenhagen problem where the two big bodies of equal masses are also magnetic dipoles and we study some aspects of the dynamics of a charged particle which moves in the electromagnetic field produced by the primaries. We investigate the equilibrium positions of the particle and their parametric variations, as well as the basins of attraction for various numerical methods and various values of the parameter λ.


Introduction
The Copenhagen problem is a particular case of the famous restricted three-body problem where the two big spherical and homogeneous primaries have equal masses and rotate with constant angular velocity in circular orbits around their centre of mass, while a small massless particle moves under the resultant Newtonian action of the two primaries.The problem was exhaustively studied by Stromgren and his collaborators in the decade of 1920's and has been revisited during the last ten years after the discovery of many exosolar planetary systems.A large part of these systems that have been identified since 1990, are two-body systems where the two primaries have almost equal masses.In this work we deal with a very interesting version of the Copenhagen problem in which we consider that the two big primaries are also magnetic dipoles and that the particle is a small charge which moves in the electromagnetic field produced by the two rotating dipoles.This new consideration besides the application to the aforementioned dynamical systems could also be applied to magnetic binary stars, the existence of which has been already confirmed.The dynamics of a charged particle in a magnetic binary system has been studied in the past ( [1,2]).The material of this work is organized in four sections.In the first section we extract the equations of motion of the charged particle, in the second section we study the equilibrium positions, their stability and their parametric variation, in the third section we investigate the basins of attraction of these equilibria by applying two well known numerical meth-ods, the Newton's method and an improved Modified Broyden's one.Finally, in the last section, we discuss the obtained results and we expose some of the major conclusions concerning our investigation.

The Problem
We consider a synodic coordinate system Oxyz centered at the mass centre O of the system, where the plane Oxy is the plane of the primaries' motion, and the x-axis is the axis of syzygies of the primaries which create dipoletype magnetic fields.We assume that the axes of the magnetic moments are perpendicular to the Oxy plane and that the primaries rotate about their centre of mass in circular orbits with constant angular velocity ω (Figure 1).That is The magnetic inductions of the dipoles are, where i A are the vector potentials.From Faraday's law we obtain A small particle of charge q and mass m is acted upon by a Lorentz's force created by the two rotating primaries where Under these assumptions, the normalized differential equations that describe the planar motion of S in the synodic coordinate system Oxyz, are, , where is the potential function, A x , A y are the x and y-components of the vector potential of the resultant electromagnetic field expressed in the synodic system, with where r 1 and r 2 are the distances of the particle from the primaries, The dynamical system is characterized by the parameter which is the ratio of the two magnetic moments.Since parameter λ connects the magnitudes of these moments, its role in the formation of the resultant electromagnetic field is very important.From Equation (1) we obtain a Jacobian-type integral of motion which has the form where C is a constant that is called the Jacobian constant or energy constant.

Equilibrium Locations and Their Parametric Variation
As it is known, in an equilibrium position on the plane of the primaries' revolution the following conditions hold, Therefore, the coordinates of these positions are the solutions of the system of non-linear algebraic equations, 0, 0 We have found that for our case where the magnetic moments are perpendicular on the plane of the primaries' revolution, all the equilibrium positions of the particle are located on this plane.Some of these points are dynamically equivalent, which means that are characterized by the same Jacobian constant C and by the same state of stability.Τhe total number of the existing equilibrium locations depends on the value of parameter λ.More precisely:  When .6 0 2    there are three collinear equilibrium positions (L 1 , L 2 , L 3 ). When 2.6 ≤ λ ≤ 10.8, there are three collinear positions (L 1 , L 2 , L 3 ) and four triangular positions that form two pairs (L 4 , L 5 ) and (L 6 , L 7 ).The points of each pair are symmetric with respect of the x-axis and are dynamically equivalent since they are characterized by the same Jacobian constant C and the same state of stability. When λ > 10.8 there are three collinear points (L 1 , L 2 , L 3 ) and one pair of dynamically equivalent triangular points (L 4 , L 5 ). Figure 2 shows the distribution of these equilibria on the Oxy plane for various λ.

Methods and Algorithms for Finding the
Basins of Attraction

Basins of Attraction
As we have mentioned before, in order to locate the equilibrium positions of the particle, we have to solve numerically the non-linear algebraic system (4).This is achieved by applying an iterative scheme, provided that an initial approximation is given.This starting value represents a point on plane (x, y).The iterator stops when some predetermined accuracy is reached for an equilibrium point which can be considered as an "attractor" of the method.We call the set of the initial points that lead to the discrete equilibrium points or the dynamically equivalent points, "basins of attraction" (or "basins of convergence", or "attracting domains").These regions have been studied in the past in some problems of Celestial Dynamics, such as the restricted three-body problem [3], the ring problem of (N + 1)-bodies [4][5][6][7][8], and the Hill's problem with radiation and oblateness [9].The technique to find these regions is based on a double scanning process of the Oxy plane.Here, we have considered the intervals for the vertical and the horizontal scanning respectively and a step equal to 0.005.We have created in this way a dense grid, the nodes of which are the initial values of the algorithm in use.In all the considered cases we have used 10 −8 as the accuracy criterion to terminate the iterative process.We have stored in separate files all the initial points that lead to dynamically equivalent equilibrium points.At the same time, for each initial point, we have recorded the number of the iterations required to find the equilibrium position with the aforementioned accuracy.Naturally, the number of the iterations required to locate an equilibrium position depends on the predetermined accuracy, while the number of the points that constitute the set of the attracting domain of dynamically equivalent set of equilibrium positions depends on the step-size of the scanning process, and the method used.At this point we notice that many numerical methods solving systems of nonlinear algebraic equations have been proposed and described in the past (see for example [10,11]).In what follows we briefly describe the two algorithms used, that is, the old classic Newton's algorithm and a modified version of Broyden's iteration scheme.

Newton's Algorithm
For our particular problem this very popular and simple algorithm takes the form, , , are the values of x and y at the n-th approximation of the iteration process.Relations (5) can also be written in the form,

Modified Broyden's Method
Broyden's method is based on the secant's one (it is a generalization of it) and belongs to a class of methods called quasi-Newton.We suppose that and we compute the Jacobian matrix and its inverse We calculate the next approximation   1 x in the same manner as in Newton's method and then we determine each step, until the criterion of convergence is satisfied.
In this method 1 n A  has replaced the Jacobian matrix   that appears in Newton-Raphson's method.This is the main difference between the two methods.Broyden's method has convergence that is called "superlinear"; this implies that On the contrary, Newton's method has a "quadratic" convergence.For our problem, for Broyden's method, where 1 n A  is given by the formula: where . This method has been applied by Gousidou and Kalvouridis [5][6][7][8].
An improvement of Broyden's method is based on the Sherman-Morrison's matrix inversion formula.This formula computes the inverse matrix  , eliminating in this way, the need for a matrix inversion in each iterative step.So, we have not to inverse the matrix 1 i A  , which appears in Newton's method and not need to solve a system in each iteration, as it is needed in Broyden's method.Its convergence is also "superlinear", but with less number of operations than the original Broyden's method, since the operations that are needed are those for the multiplication of a matrix by a vector, as it is resulted from the formula: (Faires and Burden [12]).
ments number of rts of the

General Remarks and Com
Newton-Raphson's method gives a greater converging points per equilibrium point, in comparison to Broyden's method which gives a few number of converging points.This is shown in the bar charts of Figure 3 for the three considered cases with λ = 2, 10 and 12.
The general structure of the obtained basins of attracttions generally consists of some very dense parts that mainly evolve around the respective equilibrium point and of dispersed points that lie on the boundaries of the dense regions of this equilibrium or other ones.Figures 4-7 illustrate the results obtained from our investigation.The two methods provide results showing similarities and differences that mostly arise from the philosophy and the geometric origin of the individual method.
Fundamental similarities:  of the central pa In all diagrams, the shapes corresponding attracting regions that evolve in the neighborhood of the dynamically equivalent equilibrium points are similar. The boundaries of these central parts are not clearly defined.They look as if they are being unraveled or unweaved.If one chooses launching points near these boundaries, very slight differences can give rise to very different outcomes.Substantial differences: e The central parts of th Raphson's method are "compact" and homogeneous, that is, all their points lead to the equilibrium positions of a particular set and therefore show the deterministic aspect of these areas for the particular numerical method.On the contrary, the corresponding regions obtained by means of modified Broyden's method are less extended and less cohesive than the previous ones. The initial approximations of plane (x, y) which do not lead to a target in the Newton-Raphson method, are very few.On the contrary, in the modified Broyden's method, a great number of non-convergent initial approximations exists.Most of them are located at a relatively great distance from the primaries and are distributed in some extended areas of the diagram.Newton-Raphson's method "strikes" the target very st, after a very few iterations.On the contrary, Broymethod needs many more iterations to achieve 1i Copyright © 2012 SciRes.AM convergence.This is depicted in the bar charts of 3 tained for the three considered cases with λ = 2, 10 and 12.

5.2.
We have extended our investigation by studying th which the attrac or set is resolved into smaller regions, according to the number of the iterations needed to reach the equilibrium  The general observations that concern all sets of points and both methods can be summarized as follows:  The subset of the "launching" points of (x, y) plane corresponding to the class interval 1 -10 iterations for any equilibrium position of any equilibrium set, consists of a small region which occupies the central part    ng terval 40 itera ns, it merely consists of dispersed points lying either on the boundaries ted the attracting areas for the equilibrticle in a magnetic binary system us-ues of parameter λ nd the at acting domain of ea dynamically equivalent set of equilibria consists of a dense Regardi class in > tio of the dense regions of the previously mentioned class interval, or between the dense regions of the other equilibrium points or sets.Naturally, a decreasing density is observed as the number of iterations increases.

Conclusion
We have investiga ria of a charged pa ing both Newton's and Modified Broyden's methods.In this system the number of equilibria depends on the val-sub-domain in the close neighbourhood of each equilibrium point of this set and of dispersed points distributed a tr ch at the boundaries of the former domain of the same or other equilibrium sets.These boundaries are not clear and in some cases present a fractal structure.The application of the two methods reveals similarities and differences that exist between them.The dense part of an attracting region calculated with Newton's method is more "compact" than the respective one calculated with Modified Broyden's method.Another difference concerns the dispersed points.When these points are calculated with Newton's method they are distributed in an extended area of the diagram and mainly on the boundaries of the "compact" regions.In Modified Broyden's method these points are very few and in many cases form small point clusters beyond or near the central domain.We also observe that Newton's method "strikes" the final target very fast.For the vast majority of the initial approximations, the method needs only 1 -10 iterations.On the contrary, Modified Broyden's method needs often many more iterations to achieve convergence.Figures 3-7 and Table 1 justify our remarks.The presented material confirms that Newton's method is more efficient than Modified Broyden's one (for this type problem) as regards all the indices examined, such as the speed of convergence and the total number of launching points of the xy-plane that lead to the desired final destination.

P 1 and P 2 BFigure 1 .
Figure 1.The configuration of the primaries-magnetic dipoles and the small body S.

Figure 2 .
Figure 2. Distribution of the equilibrium positions for various values of λ.

Figure 3 .
Figure 3. Bar charts showing the number of convergi points per equilibrium poin t in each nu-

Figure 7 .
Figure 7. λ = 10.Modified Broyden's method.Basins of attraction of the equilibrium 1 ; (b) L 3 ; (c) L 4 .osition of this particular set and of some dispersed points: (a) L of the attracting domain around each equilibrium p points, which appear near these areas or between the dense regions of other equilibria or sets.Regarding Modified Broyden's method this class appears as a sparse cluster of points that are scattered among the points of the second class.The subset of the points corresponding to the class interval 11 -20 iterations o gions or points which complement the central regions of the first class interval and of a large number of isolated points in the attracting area of this particular equilibrium point or set.These points are spread in a widely extended area on (x, y) plane. The points corresponding to the class interval 21 -40 iterations of each set evolve in a simila