FRACTAL BASINS OF ATTRACTION IN A BINARY QUASAR MODEL

In Fig.1(a) and Fig.1(b), we have presented the effect of parameter $R$ and effect of parameters $\text{M}_{n1},\text{M}_{d1},\text{M}_{n2}$, $\text{M}_{d2}$ on the positions of libration points, respectively. In Fig.1(a), we notice that the libration points are moving away from the origin as the value of $R$ increases in the interval [1, 3.5]. We find that L1 has less movement compared to L4, L5. The displacement of L2, L3 is larger than L1, L3 and L4. In Fig. 1(b), we are decreasing the mass of galaxy G1 and increasing the mass of galaxy G2. We observe shifts in all libration points towards the negative direction of the $x$-axis. The shift of L1 is comparatively more than the shift of L2, L3, L4 and L5. Thus the parameters have a considerable impact on the position of libration points.

In the equation (9), if we take $\dot{x}=0$ and $\dot{y}=0$ then the expression $2\Omega-C=0$ represents the zero-velocity curves. If we consider the expression $2\Omega-C\geq 0$ then it will represent the permissible region of motion where the movement of the charged particle take place.

In Fig. 4(a-f), the values of the fixed parameters are given. At these values, we find the position of all libration points. Five libration points (Li, i=1….5) are located. There are four different values of Ci’s corresponding to Li’s (Jacobi integral) as C4 and C5 are equal. The values of Ci’s are $\text{C1}=6.86825,\text{C2}=6.5988,\text{C3}=5.96756$ and $\text{C4}=\text{C5}=5.22324$. We can divide the whole configuration space $(x,y)$ into three parts: interior regions, outlying regions and forbidden regions of motion. There are five different cases:

• •

  If $C>C1$ then all necks are closed and therefore orbit can move around primaries or in outlying regions.

• •

  When $C2<C\leq C1$, no necks are open but interior regions around $m_{1}$ and $m_{2}$ get connected. Therefore motion is possible around both primaries, i.e., around interior regions or in exterior regions. However, the orbit can not escape from inner region to outer region or vice versa.

• •

  When $C3<C\leq C2$, one neck is open around $m_{2}$ allowing orbits to move close to $m_{2}$ and escape to exterior region.

• •

  When $C4<C\leq C3$, both necks are open; therefore, the movement of orbit can take place close to both primaries and escape to the outer region.

• •

  When $C\leq C4$, the orbit can move through the whole configuration space.

We have shown five different ((d) and (e) are of same cases) cases of zero velocity curves in Fig 2(a-f). In Fig. 2(a), we notice interior regions around $m_{1}$ and $m_{2}$ and also the exterior regions. The orbit can move near primaries and not around any libration points. In Fig. 2(b), we see that the interior regions around $m_{1}$ and $m_{2}$ get connected. So the orbit can move in the locality of primaries $m_{1}$ and $m_{2}$. The orbit can also move around $L_{1}$. In Fig. 2(c), there is one exit channel around $m_{2}$ and L2. Therefore the orbit can move close to $m_{1}$ and $m_{2}$, and through exit channel, it can escape to outlying exterior regions. The movement of orbit is possible around L1 and L2. Now, we decrease the value of C (Jacobi’s integral) gradually. In Fig. 2(d), at the value of $C=3.2$, there are two exit channels and also there are no forbidden regions around $m_{1}$, $m_{2}$, L1, L2 and L3. Thus we notice the increase in the region of motion, but still, the motion is not possible around the triangular liberation points. The same situation can be further observed in Fig. 2 (e). At the value $C=5.15$, the forbidden region of motion vanishes completely. Therefore the orbit can move anywhere in the configuration space $(x,y)$(see Fig. 2 (f)). Thus, we observe a significant impact on the geometry of zero-velocity curves due to variation of Jacobi’s constant.

We can determine different aspects of the dynamical system with the help of the NR-BoA. In the recent past, just a few researchers have applied this method in various dynamical system including different disturbing terms in the effective potential (for e.g. Zotos [2018], Suraj [2019], Zotos [2017], Kalvouridis et al. [2012], Sprott et al. [2015] and their references). We use NR method (multivariate form) to study the BoA related to libration points. To reveal the domain of convergence for a particular libration point, we examine a set of initial conditions. To solve the systems of bivariate function $f(\textbf{X})=0$, we apply the iterative scheme

where $J$ is the Jacobian matrix of $f(\textbf{X}_{n})$.
In this work, the system of differential equations are given by

With elementary calculations, we get the iterative formula for each coordinate as:

where ${x}_{n}$, and ${y}_{n}$ are $n$-th step of the NR method. Partial derivatives of $\Omega(x,y)$ are given in the form of subscripts. The algorithm for computing BoA, we have applied the following steps: {itemlist}

Initial conditions $({x}_{0},{y}_{0})$ are taken on the configuration plane $(x,y)$ and apply NR method. In present calculations, we have adopted a grid of $1024\times 1024$ initial conditions in an uniform way. These initial conditions are called nodes. The Min. and Max. values of ${x}$ and ${y}$ are chosen to view the complete picture of the BoA generated by the libration points.

The method is applied continuously till an accuracy of order $10^{-15}$ or the maximum number of iterations (500) is reached for each initial condition. The stopping condition $|\textbf{x}_{n+1}-\textbf{x}_{n}|\leq 10^{-15}$ or $N\leq 500$ is used.

For each initial condition, we record the number of iterations N to achieve the desired accuracy. For non-converging initial conditions, we further iterate for 10000 iterations using this method.

We fix different colours for each libration point, and a particular colour is assigned to the initial conditions according to its convergence towards a specific libration point.

After assigning all initial conditions one colour, we plot the colour coded graph, which is known as BoA. For all computations and simulations, we have used Mathematica 11.0 Wolfram [2017].

In 2016, A. Daza Alvar [2016] introduced a new tool to measure unpredictability of the basin of attraction. This new tool can quantify the uncertainty of BoA, known as basin entropy. We shall briefly discuss the algorithm for the computation of basin entropy: {itemlist}

First of all, we complete the process of plotting BoA. In the configuration plane $(x,y)$, 1024$\times$1024 initial conditions are taken, each having some colour as per its convergence towards libration points.

In this step, we divide the whole region into different non-overlapping boxes to completely cover up the entire area. Each box contains precisely 25 trajectories. We have considered approximately 42000 non-overlapping boxes for this computation.

We compute the probability of colour $j$ inside each box $i$ denoted as $p_{ij}$. The gibbs entropy for each box $i$ is computed as

where $m_{i}\in[1,\text{N}_{A}]$ is the number of colours inside the box $i$ and $\text{N}_{A}$ represents the number of libration points. The probabilities $\text{p}_{ij}$ are calculated as

Eventually, due to selection of non overlapping boxes N, the entropy of the whole grid is equal to the summation of entropy of each box $i$ of the grid.

Now, we define the basin entropy as

In similar way, we define boundary basin entropy as

where $\text{N}_{b}$ denotes the number of boxes having more than one colour. If the values of $\text{S}_{b}$ and $\text{S}_{bb}$ are greater than $\log{2}$, the BoA or boundaries along BoA is fractal in nature. Now, we study the effect of the parameters $\text{R},\text{M}_{n1},\text{M}_{d1},\text{M}_{n2}\ \text{and}\ \text{M}_{d2}$ on the BoA. We have considered two cases: In first case, the effect of parameter $R$ and in the second case, effect of parameters $\text{M}_{n1},\text{M}_{d1},\text{M}_{n2}\ \text{and}\ \text{M}_{d2}$.

We begin our investigations with the case, where $R$ varies in the interval $[1,3.5]$ and we take $\text{M}_{n1}=0.08,\text{M}_{d1}=2.0,\text{M}_{n2}=0.04$ and $\text{M}_{d2}=0.6$. In Fig. 3, the BoA for libration points are depicted by five different colours. Also, the BoA cover all of the configuration plane $(x,y)$. Most importantly, we notice that a slight change in the value of the parameter $R$, the geometry of the basins changes significantly. In particular, the change in Fig. 3 (a) to Fig. 3(b) can be seen. Here the region of liberation point L5 (green) in the upper part of Fig. 3 (a) shrunk significantly in Fig. 3 (b) due to little change in the value of $R$. A similar observation can be seen for the basin corresponding to liberation point L4 (pink) in the bottom of the figures. When the value of $R$ is changing from $1$ to till $3.5$, a major change is observed in Fig. 3(a) to Fig. 3 (f). Despite having the same number of initial conditions converging towards L4 and L5, there is a regular change in the shape of BoA for different values of $R$. The domain of convergence of libration point L1 extends towards infinity whereas the domain of convergence due to libration points L2, L3, L4 and L5 is finite for such values of $R$. Further, we observe that basins for libration points L2 and L3 have the shape like insect with legs and antennas. Also, the shape of BoA corresponding to libration points L4 and L5 appears as many wings of a butterfly. In Fig. 4 (a-f), we have shown Pie-Chart for the number of iterations needed for all initial conditions taken in the configuration plane $(x,y)$. In Fig. 4 (a), it is clear that the maximum number of iterations needed are 8, 9 and 10 $(8\to 229994,9\to 213628,10\to 163375)$. Likewise, we have recorded the number of iterations taken by initial conditions to converge for different values of $R$. We notice that approximately 95 % of the initial conditions converge after 25 iterations. A maximum number of initial conditions converge after 8 to10 number of iterations.

In Fig. 5 (a-f) the relationship between numbers of iterations needed and the initial conditions considered in the configuration plane $(x,y)$ has been established. The blue tone is used to represent the iterations needed for each initial condition. We notice that the initial conditions along the boundaries of the basins require more iterations than other initial conditions. In Table 1, we have mentioned the details of convergence of initial conditions towards the libration points L1, L2, L3, L4 and L5. The value of basin entropy $\text{S}_{\text{b}}$ and boundary basin entropy $\text{S}_{\text{bb}}$ is also given in Table 1 for different values of $R$. Based on Table 1, Fig. 3, Fig. 4 and Fig. 5, we sum up the following things:

• •

  Due to changes in the value of $R$, if we move from Fig. 3(a) to Fig. 3(f), we observe a significant change in the geometry of BoA.

• •

  The region occupied by the initial conditions converging towards L4 and L5 are approximately equal for each case, although their geometry is changing in every case (see Table 1).

• •

  In all cases, the maximum number of iterations required for the convergence of $95\%$ of the initial conditions are below 25. To achieve the desired accuracy, i.e. of order $10^{-15}$, the maximum probable number of iterations lies between 8 to 10 (See Fig. 4(a-f)).

• •

  From Table 1, we see that the value of basin entropy $\text{S}_{\text{b}}$ and boundary basin entropy $\text{S}_{\text{bb}}$ is greater than $\log{(2)}$ and therefore the existence of fractal is verified. The value of boundary basin entropy $\text{S}_{\text{bb}}$ is higher than the entropy of other regions indicates the existence of higher fractal regions along the boundaries. From Table 1 and Fig. 5, we conclude that the initial conditions which lie within the fractal regions require more iterations to converge.

We continue our investigations with another case where the masses of a nucleus $\text{M}_{n1}$ and the disk $\text{M}_{d1}$of a galaxy (G1) is decreasing and the masses of nucleus $\text{M}_{n2}$ and the disk $\text{M}_{d2}$ of galaxy G2 is increasing. We have displayed BoA for six values in Fig. 6(a-f). In all cases, the value of $\text{R}=2$ is fixed. In Fig. 7(a-f), the Pie-chart related to different parameters is displayed. Fig.8(a-f) comprises of the different graphs representing the iterations needed for each initial conditions.

Further, in Fig. 6(a-f), we notice that the BoA for L2 is decreasing, and for L3, it is increasing due to variation in parameters. We observe in Fig. 6(a) few small regions in the left part which disappeared in the next Fig. 6(b) due to change in parameters. We notice such changes or other when we move from Fig. 6(a) to Fig. 6(f). Similar to the previous case, the BoA corresponding to L1 extends to infinity. The BoA of L4 and L5 remain more or less the same in all cases. The shape of BoA of libration points L4 and L5 looks like many butterfly wings. We also observe that the BoA of libration points L2 and L3 have the shape like insects with legs and antennas.

Similar to Table 1, Table 2 comprises the number of initial conditions converging towards different libration points. The value of basin entropy and boundary basin entropy is also given in Table 2. By Table 2, it is evident that the almost entire region in the BoA (except one case) is fractal. The value of $\text{S}_{\text{bb}}$ is more than the value of $\text{S}_{\text{b}}$ indicates the presence of more fractal regions along the boundaries.

The number of iterations needed for convergence towards libration points for all initial conditions varies from 5 to 30. The distribution of iterations is shown using Pie-chart in Fig. 7 (a-f). The maximum probable iterations are 7 or 8 or 9 in almost all cases. Similar to the case I, 95% initial conditions converge in 25 iterations. The number of iteration needed to convergence for each initial condition is displayed in Fig. 8(a-f). The number of required iterations is presented using a blue tone, i.e., deep blue tone represents more iteration than faded blue tone. This figure gives us a clear idea of the regions where there is a need for more number of required iterations. We summaries followings from Table 2, Fig. 6, Fig. 7 and Fig. 8.

The total area covered by BoA corresponding to L2, L3, L4 and L5 are finite whereas the area covered by a BoA of libration point L1 extends to infinity. The configuration plane $(x,y)$ is filled by BoA completely.

All initial conditions are converging towards one of the five libration points with the accuracy $10^{-15}$. We do not find any non-converging points.

Except for one value of a parameter, the existence of fractal is verified by the values of basin entropy given in Table 2. However, in this particular case, the boundaries of BoA are fractal. In all other cases, there is a presence of highly fractal regions along the boundaries (See Table 2 and Fig. 8(a-f)).

In all cases, the maximum number of iterations needed for the convergence of $95\%$ of the initial conditions are below 15. The maximum probable number of iterations needed to achieve the desired accuracy lies between 7 to10.