A Simple Analytical Formulation for Periodic Orbits in Binary Stars
aa r X i v : . [ a s t r o - ph ] O c t A Simple Analytical Formulation for Periodic Orbits in BinaryStars
Erick Nagel & Barbara Pichardo , ABSTRACT
An analytical approximation to periodic orbits in the circular restricted three-body problem is provided. The formulation given in this work is based in calcu-lations known from classical mechanics, but with the addition of the necessaryterms to give a fairly good approximation that we compare with simulations, re-sulting in a simple set of analytical expressions that solve periodic orbits on discsof binary systems without the need of solving the motion equations by numericalintegrations.
Subject headings: circumstellar matter, discs, periodic orbits – binary: stars
1. INTRODUCTION
The majority of low-mass main-sequence stars seem to be grouped in multiple systemspreferentially binary (Duquennoy & Mayor 1991, Fischer & Marcy 1992). These systemshave attracted more attention since the discovery that many T-Tauri and other pre-main-sequence binary stars possess both circumstellar and circumbinary discs from observationsof excess radiation at infrared to millimeter waves and direct images in radio (Rodr´ıguezet al.1998; for a review see Mathieu 1994, 2000). The improvement in the observationaltechniques allows, on the other hand, the testing of theories on these objects. Moreover,extrasolar planets have been found to orbit stars with a stellar companion, e.g., 16 Cygni B, τ Bootis, and 55 ρ Cancri (Butler et al.1997; Cochran et al. 1997). All these facts make thestudy of stellar discs in binary systems, as well as the possibility of stable orbits on themthat could be populated by gas or particles, a key element for better understanding stellarand planetary formation. Instituto de Astronom´ıa, Universidad Nacional Aut´onoma de M´exico, Apdo. Postal 70-264, 04510,M´exico, D. F., M´exico. Institute for Theoretical Physics, University of Zurich, Winterthurerstrasse 190, Zurich 8057, Switzer-land.
2. Methodology
In general, orbits around binary systems are calculated by approximate methods. Nu-merical calculations are the most extended. In this work we are going to use other knownmethod to search for an analytical approximation that will provide us of solutions for anyorbit in circumbinary or circumstellar discs of the circular problem. The method is basedon a perturbative analysis. In this case we choose some terms in the expanded equations ofmotion, which are relevant to the solution of the problem. The original equations are approx-imated with a set of equations that, in many cases, have an analytical solution. Comparingwith the numerical approach, this one has the advantage that the expressions are completelyanalytical and simple. We compare our calculations with numerical studies in order to pickas many terms in the expansion as we need to obtain the solution. The method can be usedas a way to find a quick answer and for studies in which an analytic expression makes theproblem manageable.
The analysis is restricted to the orbital plane, in polar coordinates the equations ofmotion are given by, ¨ r − r ˙ ψ = F r ,r ¨ ψ + 2 ˙ r ˙ ψ = F ψ , (1)where the origin of the system is located at the position of the star with mass M (hereafterwe call this star S ). It is important noticing that S can represent either the primary orthe secondary star. Here, r and ψ are the radial and angular coordinates respectively, andthe forces along this directions ( F r , F ψ ) have the form 4 – F r = − GM r + GM ∂ Φ p ∂r ,F ψ = GM r ∂ Φ p ∂ψ , with Φ p = 1( R − rR cos( ψ − Ψ) + r ) / − r cos( ψ − Ψ) R , (2)that represents the perturbing potential (divided by GM ) due to the star with mass M (hereafter, star S ). Φ p could be also expressed in terms of the Legendre Polynomials (Murray& Dermott, 1999). Here R and Ψ are the coordinates of S and r and ψ locate the particlethat is perturbed, P .The non-perturbed state corresponds to the case with Φ p = 0. Thus, the non-perturbedorbits are conics around S , which are perturbed by the presence of the other mass. Theset of equations (1) cannot be solved analytically, then we expand the perturbation given inequation (2) using r/R as a small parameter. This parameter at some points in an orbit thatbegins around a resonance is large, so this approach is not valid. This expansion improvesin precision the closer the particle is to S . Thus, equations (1) can be written as,¨ r − r ˙ ψ + GM r = GM rR (cid:20) (1 + 3 cos 2( ψ − Ψ)) + 34 rR (3 cos( ψ − Ψ) + 5 cos 3( ψ − Ψ))++ (cid:16) rR (cid:17) (cid:18)
98 + 52 cos 2( ψ − Ψ) + 358 cos 4( ψ − Ψ) (cid:19) + ... (cid:21) ,r ¨ ψ + 2 ˙ r ˙ ψ = − GM rR (cid:20) (3 sin 2( ψ − Ψ)) + 34 rR (sin( ψ − Ψ) + 5 sin 3( ψ − Ψ))++ 54 (cid:16) rR (cid:17) (cid:18) sin 2( ψ − Ψ) + 72 sin 4( ψ − Ψ) (cid:19) + ... (cid:21) , (3)which is equivalent to the expansion (6.22) in Murray & Dermott (1999).In this analysis the star S orbits S describing a circular orbit. We consider the massof the discs negligible compared with the mass of the binary system. Thus, the orbit of S is given by 5 – R = a, Ψ = r G ( M + M ) a ( t − t ) = Ω( t − t ) , (4)where a represents the fixed distance between the stars, and Ω is the angular velocity of S from the third law of Kepler. At t = t the star is in the reference line, which for a circularorbit can be any radial line that begins on S . The non-perturbed trajectory of the particle P , is a Keplerian orbit around S , in this manner the polar coordinates are expressed, r = a np ,ψ = s GM a np ( t − t ) = ω ( t − t ) , (5)where a np and ω have the same meaning of its counterpart in equation 4, but relating P and S . The perturbed position of the particle P can be written, r = a np (1 + r p ) , (6) ψ = ω ( t − t ) + ψ p , (7)where the terms ( a np r p ) and ψ p are the corrections in the trajectory due to the presence of S . We define now the parameters m o and τ as follows, m o = Ω ω − Ω , (8) τ = ( ω − Ω)( t − t ) , (9)here m o is a small parameter for r ≪ R . This condition is required for the expansion of theequations of motion (equations 3) to be useful. τ is the angle that locates the particle P ,measured from the line connecting the stars, for the non-perturbed orbit. The first step isto replace t for τ in equation (3), using equation (9). Equations (6 and 7) are substituted in 6 –equations (3) taking into account that r p and ψ p depend on τ . The resulting expressions arefinally divided by a np ( Ω m o ) . The mass is given in units of M + M , and M + M = 1, then,¨ r p − (1 + r p )(1 + m o + ˙ ψ p ) + (1 + m o ) (1 + r p ) = m o M r p )[(1 + 3 cos 2( τ + ψ p ))+ 34 M p m (1 + r p )(3 cos( τ + ψ p ) + 5 cos 3( τ + ψ p ))+ M p m (1 + r p ) ( 98 + 52 cos 2( τ + ψ p ) + 358 cos 4( τ + ψ p ))+ ... ] , (10)(1 + r p ) ¨ ψ p + 2 ˙ r p (1 + m o + ˙ ψ p ) = − m o M r p )[3 sin 2( τ + ψ p )+ 34 M p m (1 + r p )(sin( τ + ψ p ) + 5 sin 3( τ + ψ p ))+ 54 M p m (1 + r p ) (sin 2( τ + ψ p ) + 72 sin 4( τ + ψ p ))+ ... ] . (11)Equations (10 and 11) are dimensionless. They will be solved for r p and ψ p in termsof M ( M = 1 − M ), τ and m o . The value for M , characterizes the stellar masses ofthe system, τ defines the angular position of the point in the original non-perturbed orbit,and m o depends on the radial position of the non-perturbed trajectory a np , in units of theseparation between the stars, a . An analytic solution is found if an expansion for r p and ψ p in powers of m o is made as follows, r p = Σ ∞ j =2 r p j ( M , τ ) m jo , (12) ψ p = Σ ∞ j =2 ψ p j ( M , τ ) m jo . (13)Specifically, for periodic orbits, it is necessary to fulfill the next conditions, r p j ( τ + 2 π ) = r p j ( τ ) ,ψ p j ( τ + 2 π ) = ψ p j ( τ ) . (14) 7 –An useful assumption is that the angular position of the orbit, from the line that connectsthe two stars, is not perturbed by S , a fact necessarily true, which can be demonstratedusing symmetry arguments. We restrict the perturbed orbits to be symmetric with respectto the line joining the stars by, ψ p j (0) = 0 , (15)this condition is general enough since there are no restrictions provided by observationsabout the orientation of the symmetry axis in the plane of the discs with respect to the linejoining the stars.Equations (12 and 13) are substituted in equations (10 and 11) and expanded in powersof m o , ending with a polynomial in m o equal to zero. The only solution for such an equationis that every coefficient is zero. Each coefficient associated to a given power, is a secondorder differential equation. The solutions of the equations for m ko with k = 2 , ,
4, are solvedimposing the conditions given in equations (14 and 15) and are given by r p = − M HM cos τ − M cos 2 τ, (16) ψ p = − HM sin τ + 11 M τ, (17) r p = M HM − HM ) cos τ − M cos 2 τ − HM cos 3 τ, (18) ψ p = − ( 13516 HM + 2116 HM ) sin τ + 136 M sin 2 τ + 1532 HM sin 3 τ, (19) r p = 331288 M − M + 225512 H M − H M + q cos τ +( 23 M − M − H M − H M ) cos 2 τ +( 255256 HM − HM ) cos 3 τ − ( 38 M + 732 H M ) cos 4 τ, (20) ψ p = − (2 q + 51564 HM − HM ) sin τ + ( − M + 4118 M + 11251024 H M + 2532 H M )sin 2 τ +( − HM + 55128 HM ) sin 3 τ + ( 201256 M + 63256 H M ) sin 4 τ, (21)with 8 – q = 364 M H − M H + 121564 M H + 105128 M H . (22)It is worth noticing that for each order (in m k ) we obtain a system of two differentialequations that produce two integration constants that can be calculated by solving the setof equations of the next higher order, by imposing the conditions in equations 14 and 15.The parameter H comes from the approximation a np = m o H , where we will consider H constant, justified by the fact that H depends weakly on m o (Moulton 1963).Using equations (16, 17, 18, 19, 20 and 21), plus equations (12 and 13) in equations (6and 7) will let us find the perturbed trajectory of the particle P . Since a np is proportionalto m o , the radial corrections are of order m o and the angular corrections are of order m o .Thus, equations (6 and 7) represent coordinates depending on time, which means velocityand accelerations can be obtained by direct differentiation. It is worth noticing that theorbits are given in the reference system where the stars are fixed in the X axis. An exampleof the orbits and rotation curves produced by this model are given in Section 4. To calculate the total radial size of discs we will assume we are in the regime of lowpressure gas. The permitted orbits for gas to settle down there, are all those that do notintersect themselves or any other (Paczy´nski 1977, PSA05). In this paper we will use the sameapproximation to find discs radii except that the intersections are found analytically (otheranalytical approximation to the radial extension for circumstellar discs for any eccentricitycan be found in PSA05). First, we choose a given M , to fix the binary system. A physicalintersection occurs when r − r = 0 , (23) ψ − ψ = 0 , (24)where r i is the radius of the orbit given by equation (6) and ψ i spans the 2 π -angular-rangeof the same orbit. The subindex i takes the values i = 1 , a np and τ that simultaneously satisfyequations (23 and 24). Note that equation (23) can be expanded in a series of Taylor in thevariable a np . Due to the assumption that we are interested in infinitesimally close orbits, alinear approximation is good enough. In this manner, equation (23) can be transformed to 9 – drda np = 0 , (25)which is a function of M , a np and τ . Also, using equation (7), the equation (24) can bewritten as ψ p, − ψ p, = 0 . (26)From equations (17, 19 and 21) we directly identify that all the terms are proportional tosin( kτ ), which results proportional to sin τ . Thus, sin τ can be factorized in equation (26)giving the solutions τ = 0 and τ = π . Equation (26) has other solutions but none allows tofind a root for equation (25). Taking ( τ = 0, τ = π ) and using any of these options in equation(25), we end up with an equation for a np . The solution of equation (25) using τ = 0 , π givestwo values of a np , which correspond to the intersections. We take the intersections at thesmallest radius (that represent the maximum radii where gas particles would be able to settledown), thus the second intersection at larger radius are not useful in this analysis. In bothcases the intersecting orbits are tangent to each other, and contiguous orbits between themintersect at an angle. Because we do not consider this kind of orbits as part of the disc, thedisc naturally ends up in the orbits with smaller a np .The minimum value for a np comes from τ = π , this is, the intersection appears atthe opposite side of the star S . In Table 1, the solution of equation (25) for differentvalues of M , with τ = π , is given in the column named a np,τ = π . The following two columns( r ( τ = 0), r ( τ = π )) give the radial positions of the innermost intersecting orbit. The averageof these values is given in the column < r > , the next column gives the same average butfrom PSA05, and the last column is the approximation to the Roche Lobe radius by Eggleton(1983) given by, R i a ≈ R i (Egg) a = 0 . q / i . q / i + ln (1 + q / i ) , where (27) q = 1 − M M , (28) q = M − M , (29)(30)where i refers to either star S i with i = 1 ,
2. 10 –From Table 1, < r > and < r >
P SA (for particular masses of the stars, M and M = 1 − M ) coincides within a 20% of error for M ≤ .
5, decreasing to less than 10% for M > .
5. As expected the last non-intersecting orbit is contained inside the Roche lobe, aswe can see in the Figure 1 for M = 0 .
2. A good approximation for the Roche lobe is givenusing the approximation of Eggleton (1983).
3. Orbits in Circumbinary Discs
The method described in Section 2.1 is followed in this case with few differences. Firstof all, the origin is now on the center of mass of the system. The perturbing potential Φ p ,takes the form,Φ p = GM r h − R r cos( ψ − ψ )+( R r ) i − / + GM r h
1+ 2 R r cos( ψ − ψ )+( R r ) i − / , (31)where R = M R , R = M R are the distance to the star S and S from the origin,respectively. It is important to mention that S is located to the left of the center of mass,and that it can represent either the secondary or the primary star. a is the distance betweenthe stars. r and ψ − ψ are the coordinates for the particle P that trace the circumbinaryorbit, where the angle is measured from the radial vector which aims to S . In this case,equation (1) can be used with the new definitions as follows, F r = ∂ Φ p ∂r , (32)Table 1: Sizes of the circumstellar discs. M a np,τ = π r ( τ = 0) r ( τ = π ) < r > < r > P SA R RL F ψ = 1 r ∂ Φ p ∂ψ . (33)Thus, an expansion of equation (1) for R r ≪ R r ≪ r − r ˙ ψ + G ( M + M ) r = − GM p a r n
34 [1 + 3 cos 2( ψ − ψ )] + R − R r [3 cos( ψ − ψ ) + 5 cos 3( ψ − ψ )] + ... o ,r ¨ ψ + 2 ˙ r ˙ ψ = − GM p a r n
32 sin 2( ψ − ψ ) +3( R − R )8 r [sin( ψ − ψ ) + 5 sin 3( ψ − ψ )] + ... o , (34)where M p = M M = M (1 − M ).The third term on the left side of the first equation represents the contribution to theforce, in the case that all the stellar mass is concentrated in the origin. The right side ofboth equations are the perturbing force due to the mass distribution between the stellarcomponents. If the right side of equation (34) is equal to zero, we find the non-perturbedorbit in the form, r = a np , (35) ψ = s G ( M + M ) a np ( t − t ) = ω ( t − t ) . (36)Note that the definition of ω differs from the circumstellar case in the mathematicalexpression, however the meaning is the same, non-perturbed angular velocity of the particle P . Again, the perturbed position of P is given by equations (6 and 7). Equation (9) changesto τ = (Ω − ω )( t − t ) , (37)where Ω = q G ( M + M ) a is the angular velocity of the binary system, then τ represents theangle between S and P .An useful parameter is given by 12 – µ = ω Ω = ( aa np ) / , (38)where µ is the parameter used to expand equation (34), then µ is small as we are only takinginto account circumbinary orbits far away from the stars. Changing t → τ (equation 37)and expanding equation (34) in a power series of µ / , we can write the equations of motionas ¨ r p − (1 + r p )( µ − µ + ˙ ψ p ) + µ (1 − µ ) r p ) = − M p (1 − µ ) µ / (1 + r p ) n
34 [1 + 3 cos 2( τ + ψ p )] +(2 M − µ r p ) [3 cos( τ + ψ p ) + 5 cos 3( τ + ψ p )] + ... o , (39)(1 + r p ) ¨ ψ p + 2 ˙ r p ( µ − µ + ˙ ψ p ) = − M p (1 − µ ) µ / (1 + r p ) n
32 sin 2( τ + ψ p ) +38 (2 M − µ (1 + r p ) [sin( τ + ψ p ) + 5 sin 3( τ + ψ p )] + ... o . (40)In the same way as in Section 2.1, differential equations are extracted at different ordersif we expand r p and ψ p as follows, r p = Σ ∞ i =4 r p i ( M , τ ) µ i/ , (41) ψ p = Σ ∞ i =4 ψ p i ( M , τ ) µ i/ . (42)The perturbed trajectories must be closed, then equations (14) are imposed. Also, thefact that the orbit is symmetric with respect to the line joining the two stars, allows to findthe following condition, ˙ r pi ( τ = 0) = 0 . (43) 13 –Substitution of equations (41 and 42) in equations (39 and 40) gives a couple of poly-nomials in µ i/ . A solution can be found if every coefficient of the µ i/ -terms is zero. Usingconditions given in equations (43, 14) allows directly to find the set of solutions.Moulton (1963) describes the method and gives explicitly the expressions for r p i and ψ p i for i ≤
15. The integration constants with the full expressions are provided in the nextequations r p = 14 M p ,r p = 0 ,r p = 0 ,r p = 0 ,r p = 364 M p (5(1 − M ) − M p + 5 M ) ,r p = 0 ,r p = 916 M p cos 2 τ,r p = 0 ,r p = 1768 M p (175(1 − M ) − − M ) M + 711(1 − M ) M − − M ) M + 175 M ) − M p − M )(3 cos τ + 59 cos 3 τ ) ,r p = 34 M p cos 2 τ,r p = − M p + M p
64 (25(1 − M ) − M p + 25 M ) cos 2 τ +1751024 M p ((1 − M ) − M p + M ) cos 4 τ,r p = − M p (1 − M )(cos τ + 527 cos 3 τ ) , (44)and for the angle we have, ψ p = 0 ,ψ p = 0 ,ψ p = 0 ,ψ p = 0 ,
14 – ψ p = 0 ,ψ p = 0 ,ψ p = 38 M p sin 2 τ,ψ p = 0 ,ψ p = − M p (1 − M )(sin τ + 59 sin 3 τ ) ,ψ p = 316 M p sin 2 τ,ψ p = 532 M p ((1 − M ) − M p + M ) sin 2 τ + 35256 M p ((1 − M ) − M p + M ) sin 4 τ,ψ p = 14 M p (1 − M )(9 sin τ − τ ) , (45)where the integration constant for the radial equation 39 is calculated using the solution forsix orders ahead. While for the angular equation 40 the integration constant is obtainedfrom the solution for ten orders ahead.It is worth to mention that the approximation to circumbinary periodic orbits given byMoulton is far from the solution as we have compared with the numerical solution, thus wehave developed all the necessary terms in the expansion of the potential to reach a fairlygood approximation to the numerical solution of the problem.In this way, an analytical expression is found for circumbinary orbits with high accuracyusing the expressions 44, 45, and the equations (41, 42, 38, 6 and 7). Gap ) Our purpose in this section is to provide a good estimation for the radius of the innerboundary of a circumbinary disc. Then we are looking for the closest stable orbits to thebinary system. In this case, it is required to calculate orbits decreasing in size until con-secutive orbits intersect each other. The procedure we follow was to find a few new termsin the approximate solution for the disturbed orbit and calculate the larger inner orbit thatintersects a contiguous one with the method described in Section 2.2. In this case, as in thecircumstellar disc calculation (see Section 2.2), the angular correction ψ p i is proportional tosin τ , the intersections occur in the same manner at τ = 0 and τ = π , i.e., in the line joiningthe stars. The larger the number of terms considered for the approximation, the closer theorbit is to the orbit calculated in PSA05. However, the difference in sizes is large in spiteof taking into account the terms up to µ / = ( a np a ) − / that gives high precision in the 15 –circumstellar discs case. Table 2 gives for a set of values for M , the non-perturbed radius a np , the physical radius of the orbit in the intersection with the line connecting the stars, r ( τ = 0) and r ( τ = π ), and the analogous values in PSA05. Here, for each M there is anintersection at τ = 0 or at τ = π . Note that, for example, the system with M = 0 . M = 0 .
8, only rotated an angle π , then, the intersection at a np = 1 . τ = 0 in the former one, corresponds to a np = 1 . τ = π in the latter. Thus, Table2 gives radii for M ≤ . In the numerical approach, the solution for a closed orbit is searched by successiveiterations. This means that an unstable orbit (surrounded by orbits far from the solution) isvery hard to find. Unlike numerical calculations, analytically, one can calculate either stableand unstable orbits without any possible identification. One has to apply another criteria tolook for the long-lived orbits as the ones in real discs, this criteria is taken from a stabilityanalysis.In the case of circumstellar discs, the criteria to end up a disc was the intersection oforbits. In the circumbinary discs case, this is not applicable since in general orbits startbeing unstable before any intersection of orbits.Message (1959) defines an orbit as stable or unstable using the equation that describesnormal displacements from a periodic orbit in the three-body restricted problem. Thesedisplacements are solved with terms proportional to exp( icτ ), where c is given by c − (4 + λ − + λ ) c + 2( λ − λ − ) c + λ − λ − ν ν − = 0 , (46)whereTable 2: Size of the central gap in the circumbinary disc M a np,τ =0 r ( τ = 0) r ( τ = π ) r ( τ = 0) P SA r ( τ = π ) P SA λ ± = θ p − − ( β + β ± ) θ p θ p − − β ± β ± θ p θ p − , (47) ν ± = − β θ p ± , (48) β k = 1 θ p − ( k + c ) , (49)and θ p , θ p ± are the main coefficients of the Fourier expansion of Θ( τ ), where Θ( τ ) is thefunction involved in the normal equation displacement, d qdτ + Θ( τ ) q = 0 , (50)which is given in Message (1959). The function Θ( τ ) depends on the shape of the orbit,i.e. depends on τ .If for a specific orbit, c has an imaginary term then the orbit is unstable.Calculation of c is made for analytical orbits given by equations (6 and 7) with r p and ψ p estimated with equations (41, 42) and the coefficients given by equations 44 and 45. Orbitssmaller and larger that the intersecting-orbits given in Table 2 are considered. The valuestaken for the parameter a np are given in Table 3.The parameter a = a np,P SA corresponds to the analytical orbit closer to the innermostorbit of the circumbinary disc, found in PSA05, for several stellar mass values, M . Thus, ifthis set of orbits is (un)stable we expect that the numerical orbits also are (un)stable. Wesolve equation (46) for the orbits in Table 3 and the results are given in Table 4. There, thelargest value of the imaginary part of c , looking at all the modes found, are given by M ,ranging from 0 . . M ax | Im ( c ) | is to decrease when the orbit is moving awayfrom the binary system, as expected. Note that M ax | Im ( c ( a np,P SA )) | is quite small, conse-quently, this value can be taken as a threshold for instability and the orbit in PSA05 willTable 3: Values of the parameter a np for the orbits studied in the instability analysis. M a np, a np,Inters a np, a np, a np,P SA a np,
4. Test: Comparison with Numerical Results
In this section we present a simple comparison between the results presented in thiswork and precise numerical results performed with the usual methods to integrate the motionequations in a circular binary potential. We show for this purpose two tests, the first is adirect qualitative comparison between numerical and the analytical approximation to somerepresentative discs. The second including velocities, are the correspondent rotation curves.
The equations of motion of a test particle are solved for the circular restricted threebody problem using an Adams integrator (from the NAG FORTRAN library). Cartesiancoordinates are employed, with origin at the center of mass of the binary, in an inertialreference frame. Here we look for the families of circumstellar and circumbinary orbits ata given stellar phase of the stars (equivalent to look for the families in the non-inertialframe of reference). We calculate the Jacobi energy of test particles, in the non-inertialframe of reference co-moving with the stars, as a diagnostic for the quality of the numericalintegration, conserved within one part in 10 per binary period.Table 4: Maximum imaginary Part of c for the Orbits Studied in the Instability Analysis M ax | Im ( c ) | \ M c ( a np, ) 2 . × − . × − . × − . × − . × − c ( a np,Inters ) 6 . × − . × − . × − . × − . × − c ( a np, ) 2 . × − . × − . × − . × − . × − c ( a np, ) 4 . × − . × − . × − . × − . × − c ( a np,P SA ) 2 . × − . × − . × − . × − c ( a np, ) 0 0 0 0 0 18 – In the case of binary systems is well known the shape and extension of discs are differentfrom the single stars case, where the periodic orbits are circles with a radial unlimited (bythe potential of the single star) extension. In binary systems the stellar companion exertsforces able to open gaps, limiting the extension of the discs and producing at the same timedeformations in the shape of the discs making them visually different from single star cases.In this section we show a comparison between the discs obtained by the analytical approxi-mation provided in this work, the numerical solution and the keplerian approximation.To construct the analytical circumstellar discs we show in this section, we have employedthe equations 6 and 7 that can be written in the inertial reference system, r = a np (1 + r p ) , (51) ψ = τ + ψ p , (52)with τ in the interval [0 , π ].In the Figure 1, we show the comparison for circumstellar discs between the approxima-tion presented in this work for periodic orbits in a binary system (left panels) and numericalcalculations (right panels). The value for M is indicated on the top of the figure. Thedarker curve represents the Roche lobe.Likewise, in Figure 2 we show a qualitative comparison of a circumbinary disc betweenthe analytical approximation (left panels) and numerical calculations (right panels). Thedistance of the stars with respect the center of mass of the binary is shown.The discs present a slight deformation of the orbits depending on the radius: the outer-most orbits are farer from being circles. The orbits are, a very good approximation, ellipseswith an eccentricity depending on the radius to the central star, in the cincumstellar discs,and to the center of mass in the circumbinary discs case. Thus, we have calculated the m=2Fourier mode for all the orbits given by A k = 1 N N X i =0 s ( φ i ) cos ( kφ i ) ,B k = 1 N N X i =0 s ( φ i ) sin ( kφ i ) , (53) 19 – aa a a Fig. 1.— Circumstellar discs comparison between orbits with the approximation presentedin this work and numerical results described in subsection 4.1 for the case M = 0 .
2. Leftpanel shows the analytical approximation provided in this work. Right panel shows thenumerical approximation. The darker curve surrounding the discs is the Roche Lobe. Theaxes are in units of the distance between the stars, a
20 – aa a a M =0.2 M =0.2
Fig. 2.— Circumbinary discs comparison between orbits with the approximation presentedin this work and numerical results described in subsection 4.1 for the case M = 0 .
2. Leftpanel shows the analytical approximation provided in this work. Right panel shows thenumerical approximation. The axes are in units of the distance between the stars, a . Thedistance of the primary and secondary star with respect to the center of mass of the systemis indicated by the black asterisks 21 –where the index i refers to the number of evenly distributed (by interpolation) in angle pointsin a given loop, and N is the total number of points in the loop, s ( φ i ) is the distance tothe point i from the center of mass of the binary. In this manner, the average inner radiusof a circumbinary disk is given in units of the semimajor axis a by the coefficient A K with k = 0. The ellipticity ( ell ) by p A k + B k with k = 2 in the same units and transformedfirst to the ratio b i /a i , where a i , b i are the semimajor and semiminor axis of the orbit i , andfinally transformed to eccentricity ( ecc ) by, ecc = √ − ell (54)which represents a more sensitive geometrical characteristic of the orbits than ellipticity.We have not considered in this comparison higher Fourier modes since they are negligiblecompared to the m=2 mode.In the Figure 3 we show the eccentricity vs the average radius for a) the orbits producedby the analytical approximation provided in this work (continuous lines in the 9 frames),b) the full numerical solution (open triangles), and for reference c) the keplerian ( ecc = 0)solution (dashed lines). In this figure we show three different masses M = 0 . , . , . We also show here, for a quantitative comparison involving positions and velocities, therotation curves for binary systems with M = 0 . , . , . v c ) at two differentregions of the discs. For the primary discs we compared the rotation curves obtained alongthe x axis, on the left part of the discs (where the radial velocities are negative and progradewith the rotation of the binary). For the secondary we chose to compare the rotation curvesalong the x axis but on the right part of the disc (where the radial velocities are positive andprograde with the rotation of the binary). Once the velocities are calculated, we transformedthe results to the inertial frame (in the circumstellar discs) by simply adding (or subtractingdepending on the given disc) the velocity of the correspondent star, and adding or subtracting 22 –Fig. 3.— Comparison between the orbital eccentricity vs radius produced by the analyticalapproximation provided in this work (continuous lines), the full numerical solution (opentriangles), and for reference the keplerian ( ecc = 0) solution (dashed straight lines). Threeexamples of different masses, M = 0 . , . , .
4, are indicated by the upper labels at left,middle and right panels respectively. The upper frames are referred to the circumprimarydiscs, the middle to the circumsecondary discs, and the lower are the circumbinary discs 23 –also a factor Ω r , where Ω is the angular velocity of the star and r is the radius of a givenpoint in he rotation curve. The last because of the chosen reference system where this worksolves the equations, which anchorages the discs to rotate with the system.In the Figure 4 and 5 we show the outer (top panels) and inner (bottom panels) regions ofthe primary and secondary rotation curve discs, respectively. By “outer” or “inner” regionswe mean, in a reference system located in the primary (or in the secondary) star, measuringand angle from the line that joins the stars, the inner region corresponds to angle zero fromthis line, and the outer regions correspond to and angle of 180 o . For both figures we plotthe three different cases, M = 0 . , . , .
4. Three different types of lines indicate, a)thekeplerian rotation curve on the top (filled circles); b) the analytical approximation providedin this work (continuous line); c) the numerical solution (open circles). The velocity andradius are given in code units (in such a way that G = 1, M + M = 1, a = 1 and Ω = 1).The rotation curves are in approximately 70% of the total radius of the correspondentdisc, keplerian curves. For the last 30% the keplerian discs velocity are systematically overthe numerical result that solves with high precision the restricted three body problem. Theanalytical solution we provide here is very close to the numerical solution as expected.In the Figure 6 we show the rotation curves as in Figures 4 and 5, but for the circumbi-nary disc. Although in this case both approximations (the one given in this work, and thenumerical one) are very close to a keplerian rotation curve, it is not exactly keplerian. Thevelocities in the analytical approximation are sistematically under the keplerian curve andthey are practically the same as the ones provided by the numerical solution, until the end(the beginning of the gap) where the numerical solution goes slightly over the analyticalsolution. We present only one case of mass ratio because other mass ratio would give verysimilar results.In the case of the full (numerical) solution, it is worth to notice we obtain in some casesdiscotinuities on the rotation curves and orbits mainly due to resonances. This is the case,for instance, in the Figure 6, where we appreciate a discontinuity very close to the resonance3:1 ( r ∼ − . a in the figure). The discontinuity is also noticeable in the Fig. 2 at the sameradius. Due to the perturative fit performed in this work, finding this kind of discontinuitiesdue to resonances, that result in general very narrow in radius but where the physics ishighly non-linear is out of the reach of our approach. 24 –Fig. 4.— Comparison between the keplerian rotation curve (filled circles), the analyticalapproximation provided in this work (continuous line), and the numerical solution (opencircles), for the outer region (top frame), and inner region (bottom frame) of the primarydisc, in the different cases M = 0 . , . , .
4. 25 –Fig. 5.— Same as Figure 4 but for the secondary star 26 –Fig. 6.— Comparison between the keplerian rotation curve (filled circles), the analyticalapproximation provided in this work (continuous line), and the numerical solution (opencircles), for the left region (top panel), and right region (bottom panel) of the circumbinarydisc, for M = 0 .
2. 27 –
5. CONCLUSIONS
We have constructed an analytical set of equations based on perturbative analysis thatresult simple and precise to approximate the solution for the circular restricted three bodyproblem. We choose some terms in the expanded equations of motion, which are relevant tothe solution of the problem. The original equations are approximated with a set of equationsthat, in many cases, have an analytical solution. We show the goodness of our analyticalapproximation by direct comparison of geometry and rotation curves with the full numericalsolution.The set of equations we prsent can provide any periodic orbit for a binary system, eithercircumstellar or circumbinary, and not only the outer edge radius of circumstellar discs orthe inner edge radius (gap) of the circumbinary disc, without the need to solve the problemnumerically.The relations provided are simple and straightforward in such a way that our approxi-mation can be used for any application where initial conditions of periodic orbits or completeperiodic orbits are needed, or for direct study of the three body restricted circular problem.For instance, in order to introduce on a hydrodynamical or particle code, our initial condi-tions would result in much more stable discs than with keplerian initial conditions, or thanby constructing the discs directly from hydrodynamical simulations of accretion to binariesthat will require long times to obtain stable discs. Since our approximation is completelyanalytical it will have the obvious additional advantage of being computationally, extremelycheap, and easier to implement than the numerical solution.The periodic orbits respond uniquely to the binary potential and do not consider otherphysical factors such as gas pressure, or viscosity that work to build the fine details ofdiscs structure. They are however, the backbone of any potential and their shape andbehavior give the general discs phase space structure. In this manner, apart of all thepossible hydrodynamical or particle discs applications, we can directly use them to studyfrom the general discs geometry to rotation curves, or rarification and compression zones byorbital crowding, etc.
Acknowledgments
We acknowledge Linda Sparke, Antonio Peimbert and Jorge Cant´o for enlighteningdiscussions. B.P. thanks project CONACyT, Mexico, through grant 50720. 28 –
REFERENCES
Bate, M. R., 1997, MNRAS, 285, 16Bate, M. R. & Bonnell, I. A. 1997, MNRAS, 285, 33Bonnell, I. & Bastien, P. 1992, IAU Colloquium 135, 32, 206Butler, R. P., Marcy, G. W., Williams, E., Hauser, H. & Shirts, P. 1997, ApJ, 474, 115Cochran, W. D., Hatzes, A. P., Butler, R. P. & Marcy, G. W. 1997, ApJ, 483, 457Duquennoy, A. & Mayor, M. 1991, A&A, 248, 485Eggleton, P. P. 1983, ApJ, 268, 368Fischer, D. A. & Marcy, G. W. 1992, ApJ, 396, 178Holman, M. J. & Wiegert P. A. 1999, AJ, 117, 621Lubow, S. H. & Shu, F. H. 1975, ApJ, 198, 383Murray, C.D. & Dermott, S.F. 1999,
Solar System Dynamics , Cambridge University Press,p. 229.Mathieu, R. D. 1994, ARA&A, 32, 465Mathieu, R. D., Ghez, A. M., Jensen, E. L. N. & Simon, M. 2000, in Protostar and PlanetsIV, ed. V. Mannings, A. P. Boss & S. S. Russell (Tucson: Univ. Arizona Press) 731Message, P. J. 1959, AJ, 64, 226Moulton, F.R.,
Periodic Orbits , Washington: Carnegie institution of Washington, 1963c1920.Papaloizou, J. & Pringle, J. E. 1977, 181, 441Paczy´nski, B. 1977, ApJ, 216, 822Pichardo, B., Sparke, L. S. & Aguilar, L. A. 2005, MNRAS, 359, 521Rodr´ıguez, L. F., D’alessio, P., Wilner, D. J., Jo, P. T. P., Torrelles, J. M., Curiel, S., G´omez,Y., Lizano, S., Pedlars, A., Cant´o, J. & Raga, A. C. 1998, Nature, 395, 355Rudak, B., Paczy´nski, B. 1981, Acta Astron, 31, 13 29 –Wiegert, P. A., Holman, M. J. 1997, AJ, 113, 1445