The Dipole Segment Model for Axisymmetrical Elongated Asteroids

The gravitational potential of the dipole segment at a certain position ${\boldsymbol{r}}={[x,y,z]}^{{\rm{T}}}$ can be expressed in a simple form as

$$\begin{eqnarray}&&U=-\displaystyle \frac{{{Gm}}_{1}}{{r}_{1}}-\displaystyle \frac{{{Gm}}_{2}}{{r}_{2}}-G{\int }_{l}\displaystyle \frac{{{dM}}_{2}}{r},\end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, and ${r}_{i}=\parallel {{\boldsymbol{r}}}_{i}\parallel (i=1,2$) is the distance to the corresponding primary. As shown in Figure 1, the two position vectors are written as

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\boldsymbol{r}}}_{1}={[x+{l}_{1},y,z]}^{{\rm{T}}}\\ {{\boldsymbol{r}}}_{2}={[x-{l}_{2},y,z]}^{{\rm{T}}}\end{array}\right..\end{eqnarray} \tag{ 3 }$$

The line integral in the third term on the right represents the gravitational attraction of the segment, which is given with

$$\begin{eqnarray}&&{\int }_{l}\displaystyle \frac{{{dM}}_{2}}{r}=\rho {\int }_{l}\displaystyle \frac{{dl}}{r}=\displaystyle \frac{{M}_{2}}{l}\mathrm{ln}\left(\displaystyle \frac{{r}_{1}+{r}_{2}+l}{{r}_{1}+{r}_{2}-l}\right).\end{eqnarray} \tag{ 4 }$$

Substituting the above equation and mass ratios into Equation (2) yields

$$\begin{eqnarray}&&U=-{GM}\left[\displaystyle \frac{(1-{\mu }_{1})(1-{\mu }_{2})}{{r}_{1}}+\displaystyle \frac{{\mu }_{1}(1-{\mu }_{2})}{{r}_{2}}+\displaystyle \frac{{\mu }_{2}}{l}\right.\\ &&\left.\qquad \times \,\mathrm{ln}\left(\displaystyle \frac{{r}_{1}+{r}_{2}+l}{{r}_{1}+{r}_{2}-l}\right)\right].\end{eqnarray} \tag{ 5 }$$

The potential function depends on the total system mass, segment length, position vectors, and mass ratios. Particularly, the potential function is symmetrical with plane oxy and plane oxz. Therefore, application of the dipole segment model is limited to approximating an actual asteroid without considering nonsymmetrical mass distributions along the oy and oz directions.

Assume the dipole segment is rotating with a constant angular velocity. The vector equations of motion of a test particle that is attracted by the dipole segment are

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}+2{\boldsymbol{\omega }}\times \dot{{\boldsymbol{r}}}+{\boldsymbol{\omega }}\times {\boldsymbol{\omega }}\times {\boldsymbol{r}}=-{\rm{\nabla }}U\end{eqnarray} \tag{ 6 }$$

in the body-fixed frame oxyz shown in Figure 1. As the centrifugal term is conservative, the above equation is reorganized by introducing an effective potential V, such that

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}+2{\boldsymbol{\omega }}\times \dot{{\boldsymbol{r}}}=-{\rm{\nabla }}V,\qquad V=U-\displaystyle \frac{1}{2}{\parallel {\boldsymbol{\omega }}\times {\boldsymbol{r}}\parallel }^{2}.\end{eqnarray} \tag{ 7 }$$

For convenience, the equations of motion are nondimensionalized during numerical computations. The length unit LU is taken to be l and the mass unit MU is M. The time unit TU is equal to $\omega -1$, resulting in a 2π TU segment rotation period. Accordingly, the nondimensional effective potential is

$$\begin{eqnarray}V & = & -\displaystyle \frac{{x}^{2}+{y}^{2}}{2}-\kappa \left[\displaystyle \frac{(1-{\mu }_{1})(1-{\mu }_{2})}{{r}_{1}}+\displaystyle \frac{{\mu }_{1}(1-{\mu }_{2})}{{r}_{2}}\right.\\ & & +\left.{\mu }_{2}\mathrm{ln}\left(\displaystyle \frac{{r}_{1}+{r}_{2}+1}{{r}_{1}+{r}_{2}-1}\right)\right],\end{eqnarray} \tag{ 8 }$$

where the scalar parameter κ is referred to as the "force ratio", representing the ratio between the gravitational and centrifugal accelerations with

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{{GM}}{{\omega }^{2}{l}^{3}}.\end{eqnarray} \tag{ 9 }$$

Note that the variables x, y, r₁, and r₂ in Equation (8) are already in dimensionless units. The nondimensional position vectors are written by introducing a new auxiliary parameter ν = l₂/l, yielding

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\boldsymbol{r}}}_{1}={[x+1-\nu ,y,z]}^{{\rm{T}}}\\ {{\boldsymbol{r}}}_{2}={[x-\nu ,y,z]}^{{\rm{T}}}\end{array}\right..\end{eqnarray} \tag{ 10 }$$

Using the center of mass equation

$$\begin{eqnarray}&&{r}_{c}=\displaystyle \frac{{m}_{1}(-{l}_{1})+\rho {l}_{1}\left(-\tfrac{{l}_{1}}{2}\right)+\rho {l}_{2}\cdot \tfrac{{l}_{2}}{2}+{m}_{2}{l}_{2}}{M}=0,\end{eqnarray} \tag{ 11 }$$

one can obtain the following relationship:

$$\begin{eqnarray}&&\nu =1+{\mu }_{1}{\mu }_{2}-{\mu }_{1}-\displaystyle \frac{{\mu }_{2}}{2}.\end{eqnarray} \tag{ 12 }$$

Here, ν is always within the range [0.5, 1). For example, if both μ₁ and μ₂ are 0.5, the value of ν is 0.5, which is consistent with the physical intuition that the mass center should be at the midpoint of the segment. The case of ν = 1 is excluded, corresponding to the case of a single point mass with μ₁ = μ₂ = 0.

Notice that Equation (7) is expressed in the body-fixed frame and yields an autonomous system in a six-dimensional phase space [x, y, z, $\dot{x}$, $\dot{y}$, ż]^T. Therefore, an energy integral can be obtained, one that corresponds to the Jacobi constant

$$\begin{eqnarray}&&C({\boldsymbol{r}},\dot{{\boldsymbol{r}}})=\displaystyle \frac{\dot{{\boldsymbol{r}}}\cdot \dot{{\boldsymbol{r}}}}{2}+V.\end{eqnarray} \tag{ 13 }$$

The Jacobi constant determines the region of allowable motion for a test particle. Substituting Equations (8), (10), and (12) into Equation (7) yields the equation of motion for each position component,

$$\begin{eqnarray}&&\ddot{x}-2\dot{y}=-\displaystyle \frac{\partial V}{\partial x},\ddot{y}+2\dot{x}=-\displaystyle \frac{\partial V}{\partial y},\ddot{z}=-\displaystyle \frac{\partial V}{\partial z},\end{eqnarray} \tag{ 14 }$$

where the partial derivatives of the effective potential function in the right-hand terms of the above equations are

$$\begin{eqnarray}\displaystyle \frac{\partial V}{\partial x} & = & -x+\kappa \\ & & \times \,\left[\begin{array}{l}\tfrac{(1-{\mu }_{1})(1-{\mu }_{2})}{{r}_{1}^{3}}(x+1-\nu )+\tfrac{{\mu }_{1}(1-{\mu }_{2})}{{r}_{2}^{3}}(x-\nu )...\\ -\tfrac{2{\mu }_{2}}{{r}_{1}{r}_{2}}\cdot \tfrac{(x-\nu )({r}_{1}+{r}_{2})+{r}_{2}}{1-{({r}_{1}+{r}_{2})}^{2}}\end{array}\right]\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial V}{\partial y} & = & -y+\kappa y\\ & & \times \left[\tfrac{(1-{\mu }_{1})(1-{\mu }_{2})}{{r}_{1}^{3}}+\tfrac{{\mu }_{1}(1-{\mu }_{2})}{{r}_{2}^{3}}-\tfrac{2{\mu }_{2}}{{r}_{1}{r}_{2}}\cdot \tfrac{({r}_{1}+{r}_{2})}{1-{({r}_{1}+{r}_{2})}^{2}}\right]\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial V}{\partial z} & = & \kappa z\\ & & \times \left[\tfrac{(1-{\mu }_{1})(1-{\mu }_{2})}{{r}_{1}^{3}}+\tfrac{{\mu }_{1}(1-{\mu }_{2})}{{r}_{2}^{3}}-\tfrac{2{\mu }_{2}}{{r}_{1}{r}_{2}}\cdot \tfrac{({r}_{1}+{r}_{2})}{1-{({r}_{1}+{r}_{2})}^{2}}\right].\end{eqnarray} \tag{ 17 }$$

In the potential approximation for irregular asteroids, the angular velocity ${\boldsymbol{\omega }}$ of the simple model is typically the same as the actual system as well as the total system mass (Bartczak & Breiter 2003; Elipe & Lara 2003; Bartczak et al. 2006). Then, the value for three independent parameters of the dipole segment, that is, [κ, μ₁, μ₂]^T, remains to be assigned. From Equation (9) there is a one-to-one correspondence between the force ratio κ and the characteristic length l. Accordingly, the three free parameters can also be chosen as [l, μ₁, μ₂]^T.

The value for the force ratio κ is a key parameter in analyzing the irregular gravitational field of an asteroid. In the following discussion, a method to determine the force ratio is explored. Assume the constant-density polyhedral model of the asteroid is a sufficiently accurate representation of the gravity potential. Results from the polyhedron serve as reference values. When a test particle is gravitationally attracted by an irregular asteroid within the polyhedral model, the potential in Equation (2) is changed to an integral over the volume as

$$\begin{eqnarray}&&U=-G{\iiint }_{\mathrm{body}}{r}^{-1}{{dM}}_{{\rm{P}}}=-G{\rho }_{{\rm{P}}}\cdot {\iiint }_{\mathrm{body}}\displaystyle \frac{{{dV}}_{{\rm{P}}}}{r},\end{eqnarray} \tag{ 18 }$$

where the subscript P denotes the polyhedron. The parameter ${\rho }_{{\rm{P}}}$ is a constant density and can be extracted from the integral, ${M}_{{\rm{P}}}$ is the total system mass, and r is the same distance as in Equation (2). The integral of the volume ${V}_{{\rm{P}}}$ (Werner & Scheeres 1997) is difficult to handle. By substituting Equation (18) into Equation (6), one can obtain the particle equations of motion within the polyhedral model:

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}+2{\boldsymbol{\omega }}\times \dot{{\boldsymbol{r}}}+{\boldsymbol{\omega }}\times {\boldsymbol{\omega }}\times {\boldsymbol{r}}=G{\rho }_{{\rm{P}}}\cdot {\rm{\nabla }}\left({\iiint }_{\mathrm{body}}\displaystyle \frac{{{dV}}_{{\rm{P}}}}{r}\right).\end{eqnarray} \tag{ 19 }$$

The length and time units

$$\begin{eqnarray}&&\left\{\begin{array}{l}{T}_{{\rm{P}}}={\omega }^{-1}\\ {{\rm{L}}}_{{\rm{P}}}=\sqrt[3]{\tfrac{3{V}_{{\rm{P}}}}{4\pi }}\end{array}\right.\end{eqnarray} \tag{ 20 }$$

are employed to nondimensionalize Equation (19). The length unit ${L}_{{\rm{P}}}$ corresponds to the radius of an equivalent spheroid with the same volume as the asteroid. The dimensionless form for the equations of motion in Equation (19) is

$$\begin{eqnarray}&&\ddot{{\boldsymbol{R}}}+2{\boldsymbol{\Omega }}\times \dot{{\boldsymbol{R}}}+{\boldsymbol{\Omega }}\times {\boldsymbol{\Omega }}\times {\boldsymbol{R}}={\kappa }_{{\rm{P}}}\cdot {\rm{\nabla }}\left({\iiint }_{\mathrm{body}}\displaystyle \frac{{{dV}}_{{\rm{P}}}}{R}\right),\end{eqnarray} \tag{ 21 }$$

where ${\boldsymbol{R}}$ is the dimensionless position vector ${\boldsymbol{R}}={\boldsymbol{r}}/{L}_{{\rm{P}}}$, and ${\boldsymbol{\Omega }}={[0,0,1]}^{{\rm{T}}}$ denotes the angular velocity. The dimensionless parameter ${\kappa }_{{\rm{P}}}$, similar to the force ratio in Equation (9), is

$$\begin{eqnarray}&&{\kappa }_{{\rm{P}}}=\displaystyle \frac{G{\rho }_{{\rm{P}}}}{{\omega }^{2}}.\end{eqnarray} \tag{ 22 }$$

For an actual asteroid, its rotating period can be determined with observation data. The bulk density can be obtained from the estimated volume and system mass of the asteroid based on the albedo or variations of the light curve (Magri et al. 2011). When we use the dipole segment to approximate an elongated asteroid, the system parameters are set to be M = ${M}_{{\rm{P}}}$ and l = ε ${L}_{{\rm{P}}}$, where ε is a positive coefficient. With these conditions, one can obtain

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{{GM}}{{\omega }^{2}{l}^{3}}=\displaystyle \frac{G{\rho }_{{\rm{P}}}{V}_{{\rm{P}}}}{{\omega }^{2}{(\varepsilon {L}_{{\rm{P}}})}^{3}}=\displaystyle \frac{4\pi }{3{\varepsilon }^{3}}\cdot \displaystyle \frac{G{\rho }_{{\rm{P}}}}{{\omega }^{2}}=\displaystyle \frac{4\pi }{3{\varepsilon }^{3}}\cdot {\kappa }_{{\rm{P}}}.\end{eqnarray} \tag{ 23 }$$

In the above equation, a relationship in terms of the force ratio has been established between the simplified and polyhedral models. It also reveals the relationship between the characteristic length of the simplified model and the equivalent spheroid of an asteroid. It can be treated as an improvement of previous efforts made by Elipe & Lara (2003), Hirabayashi et al. (2010), and Zeng et al. (2015).

One of the simplest cases for the dipole segment model is a degenerated, single point mass case, as displayed in Figure 1, which corresponds to μ₁ = μ₂ = 0. There is an infinite number of equilibrium points in the stationary orbit with radius r = $\sqrt[3]{{GM}/{\omega }^{2}}$. Furthermore, if ${\mu }_{1}\ne 0$ and μ₂ is still zero, the dipole segment is transformed into the classical dipole model. Moving from a single point to two rotating point masses with a fixed distance, the number of system equilibria becomes finite. Five or three points can be appear, depending on the parameters of κ and μ₁. The stability for those equilibria can be found in Prieto-Llanos & Gómez-Tierno (1994) and Hirabayashi et al. (2010). If μ₂ = 1, there are four equilibria for the traditional segment model (Romanov & Doedel 2014). A detailed discussion on the linear and nonlinear stability of the four equilibrium points refers to Riaguas et al. Bifurcations of equilibria for the above three conventional cases of the dipole segment model are briefly summarized as a reference. The number and locations of equilibria for the general dipole segment model will be studied in detail as follows.

Due to the symmetrical property with respect to the plane oxy, all of the equilibria of the dipole segment are located in the equatorial plane. The collinear equilibria can be obtained with the following conditions:

$$\begin{eqnarray}&&{\left[\ddot{{\boldsymbol{r}}},\dot{{\boldsymbol{r}}},\displaystyle \frac{\partial V}{\partial x},y,z\right]}^{{\rm{T}}}=0,\end{eqnarray} \tag{ 24 }$$

whereas the triangular (or noncollinear) equilibria must fulfill

$$\begin{eqnarray}&&{\left.{\left[\ddot{{\boldsymbol{r}}},\dot{{\boldsymbol{r}}},\displaystyle \frac{\partial V}{\partial x},\displaystyle \frac{\partial V}{\partial y},z\right]}^{{\rm{T}}}\right|}_{y\ne 0}=0.\end{eqnarray} \tag{ 25 }$$

From the above two equations, Equation (17) is always zero with z = 0. When y = 0, the collinear equilibrium points can be achieved by solving Equation (15) as a seventh-order polynomial function of x:

$$\begin{eqnarray}&&\sum _{i=0}^{7}{c}_{i}{x}_{\mathrm{CE}}^{i}=0,\end{eqnarray} \tag{ 26 }$$

where the subscript CE denotes the collinear equilibrium points. The explicit expressions of the coefficients c_i are relatively lengthy. Therefore, only the first three coefficients are demonstrated, yielding

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{c}_{7}={({s}_{1}+{s}_{2})}^{2}\\ {c}_{6}=(2-6\nu ){({s}_{1}+{s}_{2})}^{2}+2{s}_{1}({s}_{1}+{s}_{2})\\ {c}_{5}=(15{\nu }^{2}-10\nu +6){({s}_{1}+{s}_{2})}^{2}-10\nu {s}_{1}^{2}\\ \ \ -5{s}_{2}^{2}-(10\nu +6){s}_{1}{s}_{2}-1\end{array}\right.\\ &&{s}_{1}=\mathrm{sign}(x+1-\nu );{s}_{2}=\mathrm{sign}(x-\nu ).\end{eqnarray*}$$

If there is an inner equilibrium point on the segment similar to the L₁ point of CRTBP, the sum of the two distances r₁ and r₂ between the point and the two primaries is exactly unity in dimensionless units. It results in a zero value of the term 1 − (r₁ + r₂)² in the denominator of Equations (15)–(17), corresponding to a singular case. Hence, no inner equilibrium point exists for the dipole segment model. There are only two collinear equilibria for the general dipole segment with ${\mu }_{2}\in (0,1)$. The locations of these two collinear points lie within the respective ranges

$$\begin{eqnarray}&&{x}_{\mathrm{CE}1}\in (-\infty ,-1+\nu );\qquad {x}_{\mathrm{CE}2}\in (\nu ,+\infty ).\end{eqnarray} \tag{ 27 }$$

The triangular equilibria of Equation (25) can be deduced from the set of equations

$$\begin{eqnarray}&&\left\{\begin{array}{l}\tfrac{(1-{\mu }_{1})(1-{\mu }_{2})}{{r}_{1}^{3}}-\tfrac{2{\mu }_{2}}{{r}_{1}[1-{({r}_{1}+{r}_{2})}^{2}]}=\displaystyle \frac{\nu }{\kappa }\\ \tfrac{{\mu }_{1}(1-{\mu }_{2})}{{r}_{2}^{3}}-\tfrac{2{\mu }_{2}}{{r}_{2}[1-{({r}_{1}+{r}_{2})}^{2}]}=\displaystyle \frac{1-\nu }{\kappa }\end{array}\right..\end{eqnarray} \tag{ 28 }$$

It is difficult to analytically solve the above equations for the general dipole segment model. Nevertheless, its solution can be examined from particular cases. For instance, if ${\mu }_{1}\ne 0$ and μ₂ = 0, corresponding to the classical dipole model, Equation (28) can be simplified by considering Equation (12) as r₁ = r₂ = κ³. By combining Equation (10), the solution for the above equations can be given as

$$\begin{eqnarray}&&{[x,y,z]}^{{\rm{T}}}={\left[\displaystyle \frac{1}{2}-{\mu }_{1},\pm \sqrt{{\kappa }^{\tfrac{2}{3}}-\displaystyle \frac{1}{4}},0\right]}^{{\rm{T}}},\end{eqnarray} \tag{ 29 }$$

which is the same as the result reported by Prieto-Llanos & Gómez-Tierno (1994).

Moreover, if μ₂ = 1 (resulting in μ₁ = 0 and ν = 1/2), corresponding to the traditional segment, one can easily obtain r₁ = r₂ from Equation (28). For the segment, the coordinates of triangular equilibria are in the form [0, y, 0]^T. By substituting r₁ = r₂ into any one of the Equation (28), the solution is that the y coordinate satisfies the following equation:

$$\begin{eqnarray}&&{y}^{2}\sqrt{\displaystyle \frac{1}{4}+{y}^{2}}-\kappa =0,\end{eqnarray} \tag{ 30 }$$

which is consistent with the result of Romanov & Doedel (2014). Note that the dimensionless units between Romanov & Doedel (2014) and our study are different. The length of the segment is two in Romanov & Doedel (2014) and unity for us.

The triangular equilibria for the general dipole segment can be obtained by numerically solving Equation (28) via a Newton method. As an example, consider the symmetrical case with ${\boldsymbol{\kappa }}=1$ and μ₁ = μ₂ = 0.5. This choice of parameters indicates that there is no cohesive force in the model and the dipole segment is symmetrical with all three coordinate planes. Figure 2(a) illustrates the equilibrium points and section view of zero-velocity surfaces that are obtained numerically. There are four equilibria in total: E₁ [−1.1499, 0, 0]^T, E₂ [+1.1499, 0, 0]^T, E₃ [0, −0.9199, 0]^T, and E₄ [0, +0.9199, 0]^T. Zero-velocity surfaces indicate the motion boundaries with C = V in Equation (13). As the variables are in dimensionless units, the length of the dipole segment in Figure 2 is exactly unitary. A comparison of locations of equilibria with the dipole (μ₂ = 0.0) and segment (μ₂ = 1.0) are given in Table 1. As E_i (i = 1, 3) are symmetrical with E_i+1, only E₁ and E₃ are listed.

Zoom In Zoom Out Reset image size

Table 2 illustrates the displacement direction of the equilibrium points with respect to the central body as the three parameters [κ, μ₁, μ₂]^T are varied for the case ${\mu }_{2}\in (0,1)$. The symbol $\nearrow$ denotes that the corresponding parameter is increasing, while $\equiv$ is used for parameters that are fixed. Directional arrows indicate the moving direction of coordinates of the equilibria along with the variation of parameters. For instance, for fixed values of κ and μ₁, the collinear equilibria E₁ and E₂ move inward by increasing μ₂. The triangular equilibria in contrast move outward with the increase of μ₂.

Table 2.  Variation Trends of Coordinates for the Equilibria with the System Parameters for ${\mu }_{2}\in (0,1)$

	x(E1)	x(E2)	x(E3)	y(E3)
$\kappa \nearrow ,{\mu }_{1}\equiv ,{\mu }_{2}\equiv$	$\leftarrow$	$\to$	$\equiv$	$\downarrow$
$\kappa \equiv ,\,{\mu }_{1}\nearrow ,{\mu }_{2}\equiv$	$\leftarrow$	$\leftarrow$	$\leftarrow$	$\uparrow$
$\kappa \equiv ,\,{\mu }_{1}\equiv ,\,{\mu }_{2}\nearrow$	$\to$	$\leftarrow$	$\equiv$	$\downarrow$

Download table as:  ASCIITypeset image

It can also be confirmed that, for fixed mass ratios, all of the four equilibria move away from the central body with the increase of the force ratio κ. By increasing μ₁ with fixed values of κ and μ₂, the system center of mass moves toward the midpoint of the segment. Meanwhile, the x coordinates of all equilibria move along the -x direction, and the triangular points move outward along the axis oy. Figure 2(b) portrays the distribution of the equilibrium points for the special case with μ₁ = 0. The geometric shape of the near-realm, zero-velocity surface is similar to that of the asteroid Gaspra, which is very different from the symmetrical case of Figure 2(a).

The two triangular points of the classical dipole model (μ₂ = 0) vanish when the value of κ is not greater than 0.125 (Prieto-Llanos & Gómez-Tierno 1994). However, such a property is not found for the dipole segment with ${\mu }_{2}\in (0,1)$ through our numerical simulations. Figure 3 gives an example by illustrating the case of μ₂ = 0.8 for each constant value of μ₁. The lowest boundary value of κ is taken to be $2\times {10}^{-3}$ (far less than 0.125) to avoid misleading by numerical precision. The upper boundary value of κ is 0.5 with a step of $2\times {10}^{-3}$ to be gradually decreased to its end of $2\times {10}^{-3}$ in a number of total 250. Seven discretized values of μ₁ are selected with [0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5]^T. Each curve shows the variational trend of the triangular equilibrium point E₄ along with the decreasing of κ. Since the two triangular equilibrium points are symmetrical with respect to the axis ox, the other point E₃ is not given. The arrow shows the variational direction of the equilibrium point where the circle marker "o" denotes E₄.

Zoom In Zoom Out Reset image size

From Figure 3, the y coordinate of the point E₄ decreases along with the decrease of κ. Nevertheless, for all sample values of μ₁, the y coordinate of E₄ is not zero even with $\kappa =2\times {10}^{-3}$. Thus, the introduction of the parameter ${\mu }_{2}\in (0,1)$ makes the bifurcation property of the dipole segment different from the classical dipole model. In other words, when ${\mu }_{2}\ne 0$, the triangular equilibrium points do not vanish even for a very small value of κ, which is much less than the critical value of 0.125 for the rotating dipole model. It indicates that the existence of the segment between the two point masses prevents the vanishing of the triangular points.

Figure 4 shows another variational trend of E₄ by decreasing the value of μ₂ from 1.0 to 0.0 in a step of 0.1. Sample values of μ₁ are the same as given in Figure 3. The force ratio κ is taken to be 0.1 in a value less than 0.125, which is used to evaluate the bifurcation. In fact, there are topological transitions for the central simplified model corresponding to the above choice of μ₂, that is, from traditional segment (μ₂ = 1.0) to the general dipole segment and finally to the classical dipole (μ₂ = 0.0). Therefore, all seven sample curves start from the same initial position of E₄ right above the origin on the oy axis. With the decrease of μ₁, the x coordinate of E₄ increases, which is consistent with Table 2 but not the focus here. Seven circle markers "o" representing the equilibrium point are finally located on the axis ox when μ₂ = 0.0. Such a bifurcation caused by the topological transition of the model indicates that the triangular equilibrium point E₄ (as well as E₃) vanishes and a new collinear equilibrium point between the two point masses emerges. It also tells us that the characteristic of the equilibrium point changes from the triangular point to the collinear point. By considering all equilibrium points for the dipole segment, this bifurcation leads to the variation of the number of equilibria. When ${\mu }_{2}\in (0,1]$, there are in total four equilibria, with two triangular points and two collinear points. When the value of μ₂ decreases to exactly zero, three equilibrium points are obtained, located on the axis ox as collinear equilibrium points.

Zoom In Zoom Out Reset image size

As a simplified model to approximate the potential distribution of elongated asteroids, two limiting cases with $\kappa \to 0\,(\omega \to \infty$) and $\kappa \to \infty$ are excluded, corresponding to the asteroid rotating very fast and the nonrotating case. Moreover, it seems impossible to cover all cases through numerical simulations, due to additional parameters of the dipole segment model compared to the traditional segment (Romanov & Doedel 2014). Therefore, the aforementioned results summarized in Figures 1–4 and Tables 1 and 2 are only used to show a part of the complexity of the new model. More detailed discussions are expected to be made in further studies.

It is emphasized that we focus on the scenario of ${\mu }_{2}\in (0,1)$ in this section by excluding other topological cases of the dipole segment model specified in Figure 1. After obtaining the location for the equilibria and its variation trends, the stability of these points has to be determined, which is imp-ortant for understanding local topological structures. Since we already ensure that all equilibria are in the equatorial plane, let [x₀, y₀]^T denote the position vector of an equilibrium point. By considering a small displacement [ξ, η]^T from the equilibrium point such that x = x₀ + ξ and y = y₀ + η, we substitute them into Equation (14) to obtain a set of linear differential equations:

$$\begin{eqnarray}&&\left\{\begin{array}{l}\ddot{\xi }-2\omega \dot{\eta }=-{V}_{{xx}}\xi -{V}_{{xy}}\eta \\ \ddot{\eta }+2\omega \dot{\xi }=-{V}_{{xy}}\xi -{V}_{{yy}}\eta \end{array}\right..\end{eqnarray} \tag{ 31 }$$

Similar to the CRTBP, the above system can be transformed into first-order nonlinear equations by introducing ${\boldsymbol{\chi }}={[\xi ,\eta ,\dot{\xi },\dot{\eta }]}^{{\rm{T}}}$ as

$$\begin{eqnarray}\dot{{\boldsymbol{\chi }}}={\boldsymbol{\Phi }}\cdot {\boldsymbol{\chi }},\qquad {\boldsymbol{\Phi }}=\left[\begin{array}{cccc}0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ -{V}_{{xx}} & -{V}_{{xy}} & 0 & 2\\ -{V}_{{xy}} & -{V}_{{yy}} & -2 & 0\end{array}\right],\end{eqnarray} \tag{ 32 }$$

where ${\boldsymbol{\Phi }}$ is the linearized coefficient matrix, and V_ij are the elements of the Jacobian matrix of the effective potential with respect to the variables x and y. The eigenvalues ${\lambda }_{i}$ (i = 1, 2, 3, 4) determine the stability of the equilibria. An equilibrium point is stable if all values of Re(${\lambda }_{i}$) are negative. It is unstable when some of the eigenvalues have positive real parts. If all values for Re(${\lambda }_{i}$) are exactly equal to zero, the equilibrium point is also stable without considering the condition of multiple roots (Szebehely 1967).

The explicit forms of V_ij are very lengthy and are omitted for the sake of brevity. To verify the V_ij derivatives are correct, the two special cases in Table 1 are first examined. For the classical dipole model with κ = 1.0 and μ₂ = 0.0, the critical value of μ₁ that corresponds to a linearly stable E₃ (E₄) is approximately 0.03852, which coincides with the triangular points within the CRTBP model (Murray & Dermott 1999). As for the traditional segment with μ₂ = 1.0, it is found that the critical value of κ is 4.5481, which is consistent with the result of Elipe & Lara (2003). In a statement from Elipe & Lara (2003, last paragraph of the second section), the noncollinear equilibria are reported as unstable when κ is greater than 4.5481. However, the noncollinear equilibria E₃ (E₄) are actually unstable when κ is less than 4.5481. With such premises, the present work may be correctly compared with the work from Elipe & Lara (2003).

Figure 5(a) shows the linearly stable area for the two triangular equilibria E₃ (E₄). Each set of parameters [κ, μ₁, μ₂]^T covered by the surface corresponds to unstable equilibrium points. The maximum critical value of κ is 25.85 with μ₁ = 0.5 and μ₂ = 0.0. Such a case corresponds to the symmetrical dipole model without the segment. It indicates that for the case of μ₁ = 0.5 and μ₂ = 0.0, the two triangular equilibria are linearly stable when $\kappa \geqslant 25.85$ and unstable when κ < 25.85. To clearly show the variational trend of the stability, the stable area is projected on the plane κ−μ₂ as illustrated in Figure 5(b), while the projection view on the plane κ−μ₁ is given in Figure 5(c).

Zoom In Zoom Out Reset image size

There are four sample curves with μ₁ = [0.0, 0.14, 0.22, 0.5]^T in Figure 5(b). For a fixed value of μ₁, the variation of κ for μ₂ increasing from 10⁻⁵ to 1.0 is complex. When μ₁ < 0.22, the critical value of κ first increases to its maximum and then decreases to the final value of 4.5481. Its global shape is similar to a parabolic curve. It tells us that the maximum critical value of κ is greater than the segment or a single dipole in the same total system mass when μ₁ < 0.22. However, when ${\mu }_{1}\geqslant 0.22$, the critical value of κ decreases monotonically to the final value of 4.5481 along with the increase of μ₂. It should be noted that the accuracy of the critical value for μ₁ in Figure 5(b), μ₁ = 0.22, can improve with more numerical iterations.

When the value of μ₂ is fixed, the variation of the stable area for the two triangular points is summarized in Figure 5(c). As for the case of μ₂ = 1 × 10⁻⁵, the region under the curve corresponds to unstable cases of equilibrium points, which is similar to that of only the dipole model provided by Prieto-Llanos & Gómez-Tierno (1994). By increasing μ₁, the critical value of κ increases up to its maximum of ∼25.85 for μ₁ = 0.5. Notice that the case of μ₂ = 0 is excluded from this study.

When μ₂ = 0.64 and μ₁ = 0, the critical value of κ increases for the increase of μ₂. If we further increase the value of μ₂ to be greater than 0.64, the value of κ decreases until the final value of 4.5481. Combined with Figure 5(b), the critical value of κ decreases monotonically to the final value of 4.5481 by increasing μ₂ when μ₁ = 0.22. For example, when μ₂ = 0.5, the respective values of κ are 6.505 and 13.532, corresponding to the cases of μ₁ = [0.0, 0.5]^T. Particularly, when μ₂ = 0.8, the critical value for κ is nearly constant, resulting in a rectangular unstable area. Accordingly, the influence of the two point masses on the stability of the triangular points is very small when μ₂ = 0.8. In conclusion, the triangular equilibrium points of the system are conditionally stable, depending on the system parameters [μ₁, μ₂, κ]^T. According to numerical simulations, the two collinear equilibria are always unstable with two pairs of real roots in this planar case.

When using a dipole segment model to approximate the potential distribution of asteroids, five parameters [M, ω, μ₁, μ₂, l]^T are to be determined. In this study, both the system mass M and its angular velocity ω are set to be the same as the target asteroid. The other three parameters [μ₁, μ₂, l]^T are required to be obtained. According to Equation (9), the value of κ is dominated by the system mass M, the rotational angular velocity ω, and the characteristic length l. With the above choice, the force ratio κ can be solved after obtaining the segment length l. Recalling Equation (23), the required three parameters [μ₁, μ₂, l]^T or [μ₁, μ₂, κ]^T are equivalent to determine [μ₁, μ₂, ε]^T.

The optimized parameters of the dipole segment are obtained by minimizing the mean percent error of

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{d}}U=\sum _{i=1}^{N}\left[100\times \displaystyle \frac{\parallel {({\rm{\nabla }}{U}_{{\rm{P}}})}_{i}-{\rm{\nabla }}{U}_{i}\parallel }{\parallel {({\rm{\nabla }}{U}_{{\rm{P}}})}_{i}\parallel }\right]\end{eqnarray} \tag{ 33 }$$

in the "test zone." Recall that the subscript P denotes the polyhedron. The subscript i indicates the sequence number of test points with a total number of N. The percent error on the gravitational potential can be defined in a similar manner. Note that ${\rm{\nabla }}U$ is a vector but ΔdU is a scalar.

The next challenge is how to select the test zone to effectively evaluate the above percent error. A group of test points in a cuboid centered at the mass center of the asteroid can be taken to evaluate Equation (33). Particularly, test points inside the asteroid surface should be removed, as the simplified model can only be used to approximate the exterior potential. Even the mascons model cannot agree with its corresponding polyhedral model inside the body (Chanut et al. 2015). To reduce the computational effort, two steps are adopted in this preliminary study. The first is that the cuboid is replaced by a planar rectangle in the equatorial plane, indicating the test zone is reduced from 3D to the 2D case. The second is to introduce an inner box to exclude all test points at the inner part of the asteroid, which will be discussed in detail at the end of this section.

The outline boundary of the closed test zone (i.e., the planar rectangle in the equatorial plane) is selected as two or three times the longest axis of the asteroid. If the outline boundary is too far from the asteroid, its irregularity is nearly negligible, such as eight times the central body (Elipe & Lara 2003), where the asteroid can be treated as a point mass. For example, the magnitude of the gravitational acceleration near the surface of 1996 HW1 is ∼0.8 mm s⁻², as shown in Figure 6. When a test particle is 4 km away from the center of mass, the acceleration decreases to ∼0.05 mm s⁻², which is already 10 times smaller than the surface area. One main principle for the determination of the outline boundary is it must cover the exterior equilibrium points. As a general reference, the outline boundary can be set as four or five times the longest axis of the central body, but this is not further examined in the current work as it is not a fundamental problem.

Zoom In Zoom Out Reset image size

The inner box is introduced to avoid computation when a test point is inside the asteroid, referred to as "inner box removal." There could be a few choices for the inner box, including but not limited to a circle, an ellipse, or a rectangle covering the projection of the asteroid in the equatorial plane. Here, a rectangular area with the dimension of the body along the ox and oy directions is arbitrarily adopted, as portrayed in Figure 6, which is more easily applicable than a circular or an elliptical area. Although such a choice loses a few effective test points in the test zone, it reduces the computational complexity by avoiding checking the relative position between a test point and the asteroid surface. Eventually, the test zone for the acceleration error is a square area with the inner box removed (by excluding test points in the rectangle area contacting the asteroid surface in Figure 6). The square area is specified by the outline boundary, where the rectangular area is determined according to the dimension of the body along the ox and oy directions.

The gravitational accelerations generated by the polyhedral model on 100 × 100 points within the test zone in the equatorial plane are first computed and supply a reference. With the same coordinates of each point above, the accelerations are calculated that are produced by the dipole segment model with the jth set of system parameters SP_j = [μ₁, μ₂, ε]_j (1 $\leqslant$ j $\leqslant$ N_j). Here, the number N_j is determined by the discretizing numbers of the three parameters. By excluding points in the inner box contacting the asteroid surface, the scalar acceleration error ΔdU_j in Equation (31) is computed. By comparing all acceleration errors ΔdU_j (1 $\leqslant$ j $\leqslant$ N_j), easily implemented in Matlab by using the "min" function, the group of system parameters [μ₁, μ₂, ε]^T corresponding to the minimum ΔdU is obtained as the optimal solution.

To further improve the precision of system parameters, more discretized numbers near the optimal solution can be applied up to the required tolerance of parameters. If converged solutions are not unique (e.g., two groups of system parameters are nearly in the same value of ΔdU, which has not been found through our simulations), more test points can be added to further examine which is the best solution. In fact, two additional problems deserve to be considered. The first is the influence of the total number of grid points (or the grid density) on the optimal system parameters in a specified outline boundary with inner box removal. The second is how the system parameters of the dipole segment differ between two scenarios: with the inner box removal and with computing test points just outside the real shape of the asteroid in the same outline boundary. Due to the length of this paper, only the second problem is discussed in the following study.

It is emphasized that in the following numerical simulations uncertainties of the asteroid shape due to observations are not considered. All benchmark accelerations at test points are derived from a certain polyhedral model. To determine the optimal solution SP* = [μ₁, μ₂, ε]*, the enumeration method is utilized. Specifically, an investigated region is given for each parameter. For example, μ₁ is evaluated at [0, 0.5] with a step of ${10}^{-3}$ for a total number of 501, μ₂ at (0, 1.0) with a step of ${10}^{-3}$ for 999, and ε at (0, 5] with a step of ${10}^{-3}$ for 5000. Thus, the value of N_j is $2.5025\times {10}^{9}$. It means that, to obtain the required SP* in a precision of ${10}^{-3}$ for each parameter, total calculation times of $2.5025\times {10}^{9}$ regarding Equation (33) are implemented, where the number of test points is 100 × 100 (equivalently $N=1\times {10}^{4}$) if there is no inner box removal. By comparing all $2.5025\times {10}^{9}$ results, the optimal SP* is obtained.

If a higher precision is required for the parameter of SP*, the value of N_j is further increased, resulting in an increase of the calculation time. Therefore, a trade-off between the precision of the converged solution and the calculation time is necessarily made. Additionally, a technique to reduce the calculation effort can be applied that narrows the investigated region of each parameter step by step. For instance, an initial step of ${10}^{-2}$ can be applied to discretize each parameter where the upper boundary of ε can be even higher. For the converged solution at this step, taking [μ₁, μ₂, ε]^T = [0.34, 0.37, 2.16]^T as an example, investigated regions for them can be set as [0.3, 0.38], [0.33, 0.41], and [1.76, 2.56] with a discretized step of ${10}^{-3}$. Such an iterative procedure can be implemented until the converged solution is obtained with the required precision.

The near-Earth asteroid 1996 HW1 (hereafter referred to as HW1) is taken to evaluate the effectiveness of the potential appro-ximation method. A relatively low-fidelity polyhedron (Magri et al. 2011) is adopted with 1392 vertices and 2780 facets. Its overall dimensions are roughly 3.9 × 1.6 × 1.5 km. The bulk density of this asteroid is estimated to be $3.56\times {10}^{3}$ kg m⁻³, resulting in a total system mass of $1.54\times {10}^{13}$ kg. The rotating period is about 8.757 hr. The exterior equilibrium points of HW1 are illustrated in Figure 7(a) along with zero-velocity curves in the equatorial plane.

Zoom In Zoom Out Reset image size

The four equilibria in Figure 7(a) are E₁: [−3.212, −0.133, −0.002] km, E₂: [3.268, −0.084, −0.001] km, E₃: [0.150, −2.808, 0.001] km, and E₄: [0.181, 2.826, 0.000] km. These results indicate that HW1 is nearly symmetrical with respect to the equatorial plane since all z coordinates are at the level of 0.1% km. Moreover, it is also axisymmetrical with respect to the axis ox by orienting to the negative oy direction. The above locations are summarized in Table 3 by omitting the z coordinates. One reason for this omission is that the z coordinates of equilibrium points for HW1 are relatively small. Another is that equilibrium points of the dipole segment model are all in the equatorial plane. Note that the marker "■" denotes no value for the parameter.

Table 3.  Four Exterior Equilibria of 1996 HW1 for Different Models

	[μ1, μ2]	E1	E2	E3	E4
${\boldsymbol{A}}$	■, ■	[−3.212, −0.133]	[3.268, −0.084]	[0.150, −2.808]	[0.181, 2.826]
${\boldsymbol{B}}$	[0.358, 0.388]	[−3.202, 0.0]	[3.268, 0.0]	[0.153, −2.810]	[0.153, 2.810]
${\boldsymbol{C}}$	[0.461, 0.061]	[−3.602, 0.0]	[3.694, 0.0]	[0.118, −2.484]	[0.118, 2.484]
${\boldsymbol{S}}$	[0.362, 0.350]	[−3.201, 0.0]	[3.269, 0.0]	[0.161, −2.808]	[0.161, 2.808]
${\boldsymbol{A}}$. Polyhedral model (ε = ■, κ = ■, l = ■).
${\boldsymbol{B}}$. Dipole segment with inner box removal (ε = 2.187, κ = 2.394, l = 2.212 km).
${\boldsymbol{C}}$. Dipole segment without inner box removal (ε = 3.276, κ = 0.712, l = 3.314 km).
${\boldsymbol{S}}$. Dipole segment with real-shape removal (ε = 2.156, κ = 2.498, l = 2.181 km).

Download table as:  ASCIITypeset image

By using the method proposed in Section 4.1, three sets of parameters are obtained and summarized in Table 3 besides case A corresponding to the polyhedral model. Case B is with the inner box removal by excluding 20 × 47 = 940 points, while case C considers all of the 100 × 100 test points inside the outline boundary. All test points in the inner box are illustrated in Figure 7(b) with dot marks where some of them are covered by the polyhedral model of HW1. Case S is another referenced scenario with real-shape removal, where 210 test points out of the asteroid HW1 but in the inner box are also considered when calculating Equation (33). Those test points (previously in dot marks in the inner box) are marked with circular markers, which can be easily recognized. The scalar parameters ε and κ in Equation (23) are listed in Table 3 as well as μ₁, μ₂, and l. The characteristic length l of case ${\boldsymbol{B}}$ is 2.212 km, which is a little greater than half of the longest axis of HW1. The value l for case ${\boldsymbol{C}}$ is 3.314 km, which is nearly the same as the asteroid's longest axis. Case ${\boldsymbol{S}}$ is similar to Case ${\boldsymbol{B}}$ in a characteristic length of 2.181 km.

Let us first check the position error of equilibrium points for the two cases. The position error of the system equilibria is defined as

$$\begin{eqnarray}&&{R}_{\mathrm{error}}=\displaystyle \frac{\sqrt{{(x-{x}_{{\rm{P}}})}^{2}+{(y-{y}_{{\rm{P}}})}^{2}}}{{L}_{\mathrm{ast}}},\end{eqnarray} \tag{ 34 }$$

where the term ${L}_{\mathrm{ast}}$ in the denominator is the length of the longest axis of the target asteroid. For example, the longest axis of HW1, ${L}_{\mathrm{ast}}$, is approximately 3.904 km. Thus, the position errors of E_i (i = 1, 2, 3, 4) for case ${\boldsymbol{B}}$, case ${\boldsymbol{C}}$, and case ${\boldsymbol{S}}$ should be

$$\begin{eqnarray*}&&\mathrm{Case}\_{\boldsymbol{B}}:\qquad {[3.42 \% ,2.15 \% ,0.09 \% ,0.83 \% ]}^{{\rm{T}}}\\ &&\mathrm{Case}\_{\boldsymbol{C}}:\qquad {[10.55 \% ,11.12 \% ,8.34 \% ,8.91 \% ]}^{{\rm{T}}}\\ &&\mathrm{Case}\_{\boldsymbol{S}}:\qquad {[6.12 \% ,3.85 \% ,0.50 \% ,1.23 \% ]}^{{\rm{T}}}.\end{eqnarray*}$$

The position error for case ${\boldsymbol{B}}$ falls within the 1% level, which is much less than that for case ${\boldsymbol{C}}$. The locations of the corresponding equilibrium points are displayed in Figure 8(a) and Figure 8(b) for the ${\boldsymbol{B}}$ and ${\boldsymbol{C}}$ cases. The four equilibria marked by "×" are those generated by the polyhedral model (case ${\boldsymbol{A}}$), while the equilibria marked by "o" correspond to the dipole segment model. Figure 8(e) shows the distribution of equilibrium points generated by the dipole segment model with real-shape removal. Each position error of the equilibrium point from case S is larger than that from case B. The additional test points taking into account the real-shape removal case do not reduce the position error further compared to the case with inner box removal.

Zoom In Zoom Out Reset image size

Respective percent errors of acceleration in the equatorial plane are also illustrated in Figure 8 for both cases ${\boldsymbol{B}}$ and ${\boldsymbol{C}}$. As a first remark, we note that the acceleration error decreases with increasing distance away from the central body. In Figure 8(a) there are four main areas with errors $\geqslant 10 \%$ near the asteroid surface. When a test particle is 2 km away from the mass center, the acceleration error is on average less than 4% except in the left head part. The dipole segment model is plotted with its right position. Such a result is similar to that reported by Colagrossi et al. (2015) to approximate the asteroid 4179 Toutatis by using two point masses. For case ${\boldsymbol{C}}$, considering all test points inside the outline boundary, the percent errors are very high with 50% for the two collinear points. Therefore, by using the simple model to approximate the exterior potential of an asteroid, the method with inner box removal improves the model accuracy.

Figure 8(e) presents the percent error of acceleration for case S with real-shape removal. Compared to Figure 8(a), no obvious improvement can be observed by decreasing the acceleration error. For instance, taking the 1% contour line of acceleration error as an example, eight lobes are obtained near HW1 where the farthest distance away from the origin is nearly the same for both cases. The relative acceleration errors near the four equilibria are at the same level for case B and case S. Such error for case S for point E₃ is a little better than that for case B.

Besides the approximation accuracy, another concern is the computation time. Calculating the acceleration value on the 100 × 100 point grid requires about 39.796 s for the polyhedral model on a personal laptop. Thus, the computational time to obtain the acceleration on a single field point is 1/10,000 of the 39.796 s. When the polyhedron is replaced by its corresponding dipole segment model, the required computing time is about 0.031 s. It indicates that the utilization of the dipole segment model is nearly 1300 times faster than the polyhedral model.

The potential errors between the two models in two cases are also reported in Figure 8. It is obvious that the case with the inner box removal in Figure 8(c) is better than that without, in Figure 8(d). From Figure 8(c), potential errors in the near vicinity of the asteroid surface at the neck regions are up to 5%. At 2 km from the asteroid surface, the potential error decreases to the level of 0.1% ∼ 1%. Figure 8(f) summarizes the potential error for case S, which is very similar to that given in Figure 8(c). Above all, the removal of inner box from the test zone is effective in the potential approximation. It is not necessary to consider the real-shape removal for preliminary studies to obtain the required parameters of the dipole segment model to complete the exterior potential approximation. Additionally, the dipole segment model can be used to roughly approximate the potential of the asteroid HW1 by considerably reducing the computational time with respect to the polyhedral model. Thus, only the case with inner box removal is considered in the following discussions.

To show the better performance of the new model, a numerical comparison is made among these three simplified models in the potential approximation of HW1. All parameters of HW1 are the same as given in Section 4.2. The method to determine parameters of the segment or the dipole is from Section 4.1 with inner box removal from the test zone. Numerical results are summarized in Table 4 where case ${\boldsymbol{A}}$ is the same as in Table 3 as a benchmark.

Table 4.  Four Equilibria of 1996 HW1 for the Traditional Segment and the Classical Dipole

	[μ1, μ2]	E1	E2	E3	E4
${\boldsymbol{A}}$	■, ■	[−3.212, −0.133]	[3.268, −0.084]	[0.150, −2.808]	[0.181, 2.826]
${\boldsymbol{D}}$	[0.0, 1.0]	[−3.250, 0.0]	[3.250, 0.0]	[0.0, −2.824]	[0.0, 2.824]
${\boldsymbol{E}}$	[0.408, 0.0]	[−3.193, 0.0]	[3.261, 0.0]	[0.175, −2.802]	[0.175, 2.802]
${\boldsymbol{A}}$. Polyhedral model (ε = ■, κ = ■, l = ■).
${\boldsymbol{D}}$. Traditional segment with inner box removal (ε = 3.182, κ = 0.777, l = 3.219 km).
${\boldsymbol{E}}$. Classical dipole without inner box removal (ε = 1.880, κ = 3.768, l = 1.902 km).

Download table as:  ASCIITypeset image

The traditional massive, straight segment corresponds to case ${\boldsymbol{D}}$ with μ₁ = 0.0 and μ₂ = 1.0. Only one free parameter ε is optimized, yielding a value of 3.182 and a characteristic length of 3.219 km. As for the classical dipole model (μ₂ = 0.0), two optimized parameters are obtained with μ₁ = 0.408 and ε = 1.880 in case ${\boldsymbol{E}}$. It indicates that the fixed distance between the two point masses is 1.902 km. Comparing case ${\boldsymbol{B}}$ with cases ${\boldsymbol{D}}$ and ${\boldsymbol{E}}$, all parameters of the dipole segment are located in a range between the parameters of the other two models. For example, the characteristic length of the dipole segment is 2.212 km, which is greater than the 1.902 km for the classical dipole but less than the 3.219 km for the segment model. In other words, optimized parameters of the segment and the dipole are the two boundaries for those of the dipole segment model. Thus, the dipole segment cannot theoretically be worse than the other two models.

With the same definition in Equation (32), the position errors of equilibrium points for both case ${\boldsymbol{D}}$ and case ${\boldsymbol{E}}$ are

$$\begin{eqnarray*}&&\mathrm{Case}\_{\boldsymbol{D}}:\qquad {[3.54 \% ,2.20 \% ,3.86 \% ,4.64 \% ]}^{{\rm{T}}}\\ &&\mathrm{Case}\_{\boldsymbol{E}}:\qquad {[3.44 \% ,2.16 \% ,0.66 \% ,0.63 \% ]}^{{\rm{T}}},\end{eqnarray*}$$

where case ${\boldsymbol{E}}$ is better than the segment of case ${\boldsymbol{D}}$. Taking case ${\boldsymbol{B}}$ into account regarding the position error, case ${\boldsymbol{B}}$ is better than case ${\boldsymbol{D}}$, particularly in terms of the triangular points. The mean error of case ${\boldsymbol{B}}$ is a little better than that of case ${\boldsymbol{E}}$, but the position error of E₄ is a little larger than the classical dipole.

Figure 9 reports the percent error of acceleration in the equatorial plane of HW1 for case ${\boldsymbol{D}}$ and case ${\boldsymbol{E}}$. Similarly, the equilibrium points within the polyhedral model are marked as "×," while "o" denotes an equilibrium location within the simplified model. The contour lines added to Figure 9 are in the same values as Figure 8(a), where the asteroid HW1 is represented by the vertices of its polyhedron. From Figure 9(a), the percent error is still 2% ∼ 4% at 3 km away from the mass center. Particularly, there are mainly six regions near the asteroid surface with an acceleration error of $\geqslant 10 \%$. For the dipole model in Figure 9(b), it should be better than the segment where the regions with an error of 10% are smaller.

Zoom In Zoom Out Reset image size

Taking the contour line of 1% error as an example, we see it is farther than that in Figure 8(a). Such a result indicates that the dipole segment can obtain a higher accuracy for the potential approximation of HW1 than the other two models. The right positions of the two simple models are also illustrated in Figure 9. The calculation time for the above two models is nearly the same and is approximately equal to 0.03 s. With nearly the same CPU time but higher accuracy for the potential approximation of elongated asteroids, the dipole segment is a better alternative simple model than the other two traditional ones.

To evaluate the performance of the dipole segment in the potential approximation of axisymmetrical elongated asteroids, four more asteroids are taken as targets: Gaspra, Bacchus, Itokawa, and 103P/Hartley-2 (a comet nucleus). The physical properties of these small bodies are summarized in Table 5, including the estimated bulk density, three-dimensional size, system mass, and rotating period. The overall dimensions of these asteroids range from 0.56 km (Itokawa) to 21.06 km (Gaspra), resulting in levels of the system mass from ${10}^{10}$ kg to ${10}^{15}$ kg. Their rotating periods are also different, spanning 7.042 hr (Gaspra) to 18.0 hr (Hartley-2), which could be used to represent different elongated asteroids.

Table 5.  Physical Properties for the Selected Elongated Asteroids

Asteroid	Density (g cm−3)	Dimension (km)	Mass (kg)	Period (hr)	Vertices and facets
Gaspraa	2.71	21.06 × 12.84 × 10.52	${2.2910}^{15}$	7.042	2522 and 5040
Bacchusa	2.00	1.15 × 0.54 × 0.53	$2.71\times {10}^{11}$	14.90	2048 and 4092
Itokawab	1.95	0.56 × 0.30 × 0.24	$3.46\times {10}^{10}$	12.132	25350 and 49152
Hartley-2a	0.34	2.52 × 0.94 × 0.78	$2.75\times {10}^{11}$	18.00	16022 and 32040

Notes.

^aWang et al. (2014). ^bGaskell et al. (2006).

Download table as:  ASCIITypeset image

The original polyhedral sources for the asteroids Gaspra, Bacchus, and Hartley-2 can be found in Wang et al. (2014). As it was visited by the Hayabusa spacecraft, high-fidelity polyhedral models are available for the asteroid Itokawa from Gaskell et al. (2006), where the number of facets ranges from 49,152 to the highest at 3,145,728. In this study, the polyhedron with 49,152 triangular facets (25,350 vertices) is adopted for Itokawa. By using the method specified in Section 4.1 with the inner box of the test zone, optimized parameters of the dipole segment for these four asteroids are reported in Table 6.

Let us first see the value of μ₂. Different from HW1 with μ₂ = 0.388, the smallest value of μ₂ in Table 6 is 0.5 for Gaspra. It indicates that for these four asteroids the segment is the dominator of the dipole segment model. In particular, the highest value of μ₂ is 0.799 for Bacchus. Such a result can be explained by considering their shapes. For example, the asteroid HW1 is a typical contact binary with a clear neck region similar to a dumbbell (recalling Figure 7). That is why the dipole is the dominator with μ₂ = 0.388. However, for the peanut-shaped asteroids Bacchus and Itokawa seen in Figure 10, their main shape feature is elongation, which can be well approximated by using a cylinder or a straight segment.

Zoom In Zoom Out Reset image size

By considering μ₁, the values for the first three asteroids in Table 6 are all at the level of 1%. In other words, m₁ is the dominator as the dipole part of the simple model. In particular, the mass m₁ for Gaspra is nearly the same as its segment part with μ₁ = 0.022 and μ₂ = 0.5. Such a result coincides with its matchstick shape with a very large head part. A different case is Hartley-2 with μ₁ = 0.191 and μ₂ = 0.57. It tells us that the mass of the segment is only a little greater than the connected dipole.

The mean percent errors of acceleration for these four asteroids are shown in Figure 10. The asteroids are represented with vertices of their polyhedron. Schematic maps of the dipole segment model are also given for each asteroid, where the approximation model for Itokawa cannot be seen because of the large number of vertices. Similar to Figure 9, the marker "o" denotes the equilibrium points of the dipole segment, and "×" is used for the polyhedron. As a general remark, we note that the equilibrium points nearly coincide with each other for the two models of these asteroids in the equatorial plane, which will be discussed later in detail.

The innermost contour line is relatively a large percent error, around 10%. But errors are measured very near the asteroid surface. For instance, there are only two small areas covered with an error of 10% near the large head part of Gaspra. As for Bacchus, the regions with $\geqslant 10 \%$ acceleration error are almost in the inner box and, therefore, excluded from the test zone. Another contour of interest is the line marking the 1% error, which can be taken as an acceptable approximation of the potential function, whereas a 0.1% error would correspond to a good approximation. In Figure 8, taking the mass center as the origin, the acceleration error is at the level of 1% on the circle, whose radius is set to be the longest axis of the asteroid. As a better case, in Bacchus, the acceleration error is nearly 0.1% when a test particle is 1.15 km away from its mass center.

In fact, fly-around orbits very close to an asteroid are usually unstable even for retrograde orbits. For example, Sun-terminator orbits near Itokawa designed for the Hayabusa spacecraft can remain stable for larger semimajor axes up to orbital radii of 1.5 km at perihelion over six months (Scheeres et al. 2006). However, orbits in similar shapes at semimajor axes lower than 1.0 km are unstable, because of both the irregular gravitational attraction and the solar radiation pressure force. Thus, Scheeres et al. (2006) pointed out that the former orbits are the only feasible ones for a spacecraft orbiting Itokawa. Seen in Figure 10(c), when a spacecraft is 1.5 km away from Itokawa, the acceleration error has been decreased to approximately 0.1%. Such a result is expected to be a good approximation for preliminary studies or qualitative analyses.

Another challenge is the position error of equilibrium points for these sample asteroids. Numerical results are summarized in Table 7 with two rows for each asteroid. The upper one lists the equilibrium locations within the polyhedral model, while the lower is for the dipole segment. The respective position errors defined in Equation (34) for the sample asteroids are

$$\begin{eqnarray*}\mathrm{Gaspra}: & {[0.56 \% ,0.20 \% ,2.42 \% ,2.01 \% ]}^{{\rm{T}}}\\ \mathrm{Bacchus}: & {[0.70 \% ,2.00 \% ,0.44 \% ,0.53 \% ]}^{{\rm{T}}}\\ \mathrm{Itokawa}: & {[0.36 \% ,0.36 \% ,1.82 \% ,1.61 \% ]}^{{\rm{T}}}\\ \mathrm{Hartley}-2: & {[1.47 \% ,0.25 \% ,1.02 \% ,0.87 \% ]}^{{\rm{T}}},\end{eqnarray*}$$

where the maximum value of 2.42% is for the point E₃ of Gaspra. All of those position errors are at the level of 1%, and the minimum value is 0.36%.

Taking the acceleration error in Figure 10 into account, we see the best case should be the approximation of Bacchus. All its position errors are less than 1% except the collinear equilibrium point E₂. Additionally, the acceleration errors for E₁ and E₄ are less than 0.1%, and the two others are less than 1%. For such cases, the dynamical characteristics near the equilibrium points of the dipole segment could be taken as good initial guesses for real asteroids. Note that all exterior equilibrium points of the four asteroids are unstable for approximated dipole segment models, which are the same as their corresponding polyhedral models. It is undoubted that the CPU time using the dipole segment model has been greatly reduced compared to its corresponding polyhedron. Such a dipole segment model may significantly enhance the computational efficiency for complex calculations regarding orbit propagation or the search for periodic orbits (Zeng & Alfriend 2017; Zeng & Liu 2017).

The dipole segment model is effective in terms of computational time with an acceptable approximation error for axisymmetrical elongated asteroids. Generally, the position errors of exterior equilibrium points are at the level of 1% between the dipole segment and the polyhedral model. The acceleration error in the equatorial plane is reduced to 1% when a test particle is sufficiently distant from the central asteroid, as measured relative to the asteroid longest axis. The potential error is expected to be even lower than the acceleration error regarding HW1, as demonstrated in Figure 8. To further improve the approximation accuracy, more simplified models with explicit potential forms may be proposed with additional segments or dipoles (Bartczak et al. 2006). Particularly, the method for the potential approximation should be an improvement with respect to Hirabayashi et al. (2010) and Zeng et al. (2015), which uses a test zone of the gravitational field as a benchmark by removing test points in the inner box contacting the asteroid surface.