Dynamical Stability of Polar Circumbinary Orbits and Planet-Formation in Planetary Disc of 99 Herculis

A possible polar-ring debris disc, the dynamics of which can be described by the outer hierarchical restricted three-body problem, has been detected in 99 Herculis. An empirical formula on the minimum radius beyond which test particles in polar orbits can keep stable within ${10^7}$ binary periods is provided through the numerical fitting, applying to the binary eccentricity $e_{1} \in \left[ {0,0.8} \right)$ and the mass ratio of binary $ \lambda \in \left[ {0.1,1} \right]$, where $ \lambda = m_0/m_1$ (${m_0}$ and ${m_1}$ represent the masses of the two binary stars). The polar planetary disc has the lowerest statistical accretion efficiency and moderate impact frequency of collisions among planetesimals (with a radius of 1-10km) compared to that in the circumbinary coplanar disc and the standard disc around the single host star. Colliding timescales in the circumbinary disk (both polar and coplanar configuration) are longer than $10^7$ yr exceeding the dissipation timescales of the gas disc. The stochastic simulations show that successive collisions cannot make planetesimal grow up which may explain the formation of the debris disc observed in 99 Herculis.


INTRODUCTION
Up until January 2023, approximately 737 planets have been discovered in binary systems 1 . These planets in binary systems can be classified into three types according to their orbital configuration, S-type (satellite-type), P-type (planet-type) and L-type (libration-type) (Dvorak 1986). Planets in S-type orbits, which are also called circumprimary planets, encircle one of the stellar binary components with the second star considered to be a perturber. Planets in P-type orbits, also known as circumbinary planets, encircle both members of a binary. The L-type planets librate around the Lagrangian equilibrium points L 4 or L 5 , which can stably exist if the binary mass satisfies µ < 0.04, where µ = m 0 /(m 0 + m 1 ) . The currently observed exoplanets in circumbinary systems include approximately 75 P-type planets, others are S-type planets. 2 Planetary formation in binary systems is different from that in single-star systems because of the disturbance from the movements of binary stars. The planetary disc will be truncated during several binary periods (Artymowicz & Lubow 1994), and forced into eccentric and processional motions (Larwood et al. 1996;Paardekooper et al. 2008;Fragner & Nelson 2010). The planets can also be forced into noncircular orbits by the disc eccentricity. The inward migration of planets will be stalled upon entrance to the tidally-truncated inner cavity in hydrodynamical simulations on the formation of P-type planets in Kepler 16, 34 and 35 systems (Pierens & Nelson 2013;Pelupessy & Zwart 2013).
2D locally isothermal hydrodynamical simulations of circumbinary discs with embedded planets were performed by Penzlin et al. (2021). The results strongly support the assumption that planets migrate to their present locations due to planet-disc interaction. During the migration, circumbinary planets could be captured into mean motion resonances, which may be associated with the final locations (Gianuzzi et al. 2023).
Some P-type planets were detected close to the edges of stable regions, where the disturbance is remarkably powerful to form planets, even considering the most favorable case of 100% efficient dust accretion (Moriwaki & Nakagawa 2004;Paardekooper et al. 2012;Meschiari 2012;Martin & Triaud 2014). Moreover, the self-gravity of the circumbinary disc can excite the eccentricities and prevent a full alignment of the planetesimal pericenters, thus resulting in sufficiently large impact velocities among planetesimals that damage the impacting planetesimals in the current location of circumbinary planets (Marzari et al. 2013). The accretion of planetesimals is possible for a 2 /a 1 > 20 in Kepler 16 (Paardekooper et al. 2012;Marzari et al. 2012;Meschiari 2012), where a 1 and a 2 are the semi-major axes of binary and planetesimals respectively. The gravity of an axisymmetric disc strongly suppresses the eccentricities of planetesimals beyond a 2 /a 1 ≈ 10-20, facilitating the easy growth of 1-10 2 km objects (Rafikov 2013a). The critical radial distance beyond which planetesimal accretion is possible increases with rising binary eccentricities and decreasing mass ratio based on the examination of the relative velocities among accreting planetesimals (Scholl et al. 2007).
In addition to the difficulties of circumbinary planets' formation at observed locations, the stable region around the binary where planets can survive for long times must be identified as a fundamental question of celestial mechanics. Holman & Wiegert (1999) simulated the empirical criteria for the largest and the smallest stable orbits of test particles in the orbital planes of S-type and P-type binaries within 10 4 T 1 (T 1 is the binary period) in the range 0.0 ≤ e 1 ≤ 0.7-0.8 and 0.1 ≤ µ ≤ 0.5.
Many factors could influence the stability of planets during and after the post-oligarchic evolution. Hong & van Putten (2019) extended the analysis of the chaotic region of coplanar P-type orbits by Dvorak (1986) to the un-restricted three-body problem and counter-rotating orbits. Mean motion resonances between P-type planets can interact with the binary via resonant and secular effects creating additional instabilities and driving chaos in multi-planet resonant systems (Sutherland & Kratter 2019). Thun & Kley (2018) found that massive planets can significantly alter the disc structure and remain on nearly circular orbits based on the planet-disc mass ratio, while low-mass planets are strongly influenced by the disc, with eccentricities excited to high values. If additional planets formed in the circumbinary disc, planet-planet scattering which takes place near the location of the currently discovered circumbinary planets, left a single planet with low eccentricity with 90% possibility (Gong 2016).
Most of the circumbinary planets were detected by transiting and eclipse time variation. Detected circumbinary planets were usually in the coplanar plane with binary orbits due to the restrictions of the two observation methods. Interestingly, several misaligned circumbinary planetary discs have been detected. The precessional circumbinary-ring model which is mildly misaligned with the binary orbital plane by 10 • ∼ 20 • has successfully interpreted the observations of KH 15D from 1995 to 2012 (Chiang & Murray-Clay 2004;Winn et al. 2004;Capelo et al. 2012). Lacour et al. (2016) found that the binary orbital plane of HD 142527 inclines 70 • relative to the outer circumbinary disc, which was considered as a transition disc. Through 3D hydrodynamical simulations on HD142527, Price et al. (2018) confirmed that all of the main observational features such as the spirals, shadows, and horseshoe can be explained by the interaction between the disc and the inner binary. However, there is no consensus on the inclination of the disc which was considered to account for the optical asymmetry of dust. Different methods find different solutions, which mainly range from 20 • to 28 • (Avenhaus et al. 2017;Hunziker et al. 2021).
The circumbinary debris disc in the 99 Herculis system was resolved, which may move in the plane perpendicular to the binary pericenter direction (Kennedy et al. 2012). The misalignment δi between the young circumbinary protoplanetary disc around GG Tau A binary system in the quadruple system GG Tau and the binary orbit is approximately 25 • ∼ 50 • (Cazzoletti et al. 2017;Aly et al. 2018).
Each stellar component of IRS43 has its circumstellar disc, and both are surrounded by a highlyinclined circumbinary disc (> 60 • ) (Brinch et al. 2016). An unusual gas-rich circumbinary disc in the young HD 98800 system is probably in the polar configuration based on the simulation of the disc dynamics, and the physical properties of such disc are similar to those around young single stars (Verrier & Evans 2008;Kennedy et al. 2019). The disc size (5-5.5 au) make it one of the smallest discs known (Ribas et al. 2018). Ziglin (1975) described the doubly averaged outer restricted elliptic three-body problem considering quadrupole approximation for the first time. HD 98800, which is a quadruple system, can be well approximated as a hierarchical triple-star system. High-inclination particles were found in long-term stable orbits inclined by 55 • ∼ 135 • to the inner binary. During the study on disc formation, Verrier & Evans (2009) found that inclination variations and nodal precession caused by the inner binary for certain initial longitudes can suppress Kozai cycles that would otherwise occur due to the outer star in the hierarchical three-body problem. Farago & Laskar (2010) provided a complete analytical description of the test particle in the secular and quadrupolar approximations of the outer hierarchical three-body problem using a vectorial formalism.
Given the stability of the inclined orbits, Pilat-Lohinger et al. (2003) recorded the escape time for inclined P-type orbits (0 • ≤ i ≤ 50 • ) in equal-mass binary systems with the binary eccentricity ranging from 0 to 0.5 within the integration time 5 × 10 4 T 1 , and distinguished different types of motions using the fast Lyapunov Indicator. Doolin & Blundell (2011) studied the stability of inclined orbits of circumbinary test particles and found the critical radius of stable planets in polar orbits is smaller than that in coplanar orbits. Considering the mass of the planet as well as the interaction between the planet and binary, Chen et al. (2019) extend the polar orbit to the generalized polar orbit by stationary inclinations, where the precession rates of the binary and planet are the same, and the relative inclination between the orbital plane of the planet and the binary are fixed. In contrast to that retrograde circulating orbits are usually the most stable around binary with small eccentricity, polar planets around highly eccentric are the most stable (Chen et al. 2020). N-body simulations were conducted to scan ∆e, ∆Ω 2 , and chaos indicators to study the global stability for different λ of the binary and initial Ω of the polar orbits by Cuello & Giuppone (2019).
However, an empirical formula for the critical radius of stable polar orbits has not been obtained, and the possibility of planet formation in circumbinary polar discs remains unknown. The paper aims to address the two issues. Firstly, the motions of test particles in polar orbits were briefly described considering the analysis of the elliptically restricted three-body problem in Section 2, which displays a libration mechanism in the longitude of the ascending node and the inclination relative to the plane of the binary. Then, an empirical formula of the stable boundary of circumbinary test particles in polar orbits with the longest integration time 10 7 T 1 was presented in Section 3. The binary eccentricity is in the range of [0, 0.8), and the mass ratio of binary star λ = [0.1, 1]. Lastly, statistical analyses and stochastic simulation of collisions among planetesimals in polar-ring discs were conducted in comparison with the results of the coplanar circumbinary disc and standard disc around the single star in Section 4.

THREE-BODY PROBLEM
The complete Hamiltonian of the hierarchical three-body system can be described in Jacobian coordinates. The hierarchical three-body secular approximation (Harrington 1968(Harrington , 1969 can be obtained by adopting Delaunay's canonical elements, considering the quadrupole moment, and averaging Hamiltonian over the short timescales by von Ziepel transformation (Kozai 1962;Harrington 1968). In the outer-restricted three-body problem, the outer body is assumed to be a massless test particle revolving around the binary pair, thus the outer body does not affect the inner orbit. Let the longitude and pericentre of the inner orbit satisfy g 1 + h 1 = π without losing generality. This condition, combined with the nodal difference between the two orbits in the invariable plane reference system meets h 1 − h 2 = π, leads to some orbital elements of the inner binary can be eliminated from the general quadrupole Hamiltonian (Ford et al. 2000;Naoz et al. 2013). Omitting constant terms, a " simplified Hamiltonian" is obtained as follows, where, e 1 is the eccentricity of the inner binary orbit, and Ω 2 is the longitude of outer orbit ascending node. i 2 is the inclination of the outer orbit relative to the orbital plane of the binary.
The range of H is [0, 1 + 4e 2 1 ]. Separatrix between circulation and libration of Ω 2 is decided by the following, The ascending node of the test particle librates near 90 • or −90 • when 1 + 4e 2 1 > H > H c .
Simultaneously, the inclination oscillates. The maximum and minimum of inclination appear when A possible polar-ring debris disc was detected in 99 Herculis by the Herschel Key Program DEBRIS (Dust Emission via a Bias-free Reconnaissance in the Infrared/Submillimeter). This program characterized extrasolar analogs to the asteroid and Kuiper belts of the Solar System, which are collectively called "debris discs" (Pilbratt et al. 2010). 99 Herculis is a binary system, comprising an The theoretical analysis of the outer-restricted hierarchical problem is provided in this section to describe the circumbinary polar orbit in quadrupole approximation briefly. However, the longterm stability of polar orbits remains unknown. Only the quadrupole term is considered for the theoretical analysis based on the expansions of complete Hamiltonian in the ratio of the semi-major axis of the inner orbit to that of the outer one. The orbits may show different movements in the octupole Hamiltonian. For example, the eccentricity of near-polar orbits may be significantly excited (Li et al. 2014). Therefore, numeric calculations are necessary to investigate the stability of polar orbits with different orbital elements. Next, the stable region of test particles in polar orbits of the outer-restricted hierarchical three-body problem will be presented. Doolin & Blundell (2011) seem in conflict, but they are not. The critical inclination for libration, Pilat-Lohinger et al. (2003) did not fix the longitude of ascending nodes of planetary orbits specifically according to the separatrix between circulation and libration. They chose the Ω 2 arbitrarily, the critical inclination turning into the libration zone can be quite larger than i c , as revealed in Figure 1. The critical inclination turning into libration for the arbitrary node can be obtained by the following, Doolin & Blundell (2011) presented a density plot of the stability measure across entire parameter spaces with the longest integration time 5 × 10 4 T 1 . Firstly, the peninsulas of instability in the libration region which appear symmetrically on either side of i = π/2 for the binary eccentricity e 1 ≥ 0.3 converge upon each other as e 1 → 0.6. This result implies that high binary eccentricity e 1 ≥ 0.6 could be more favourable to the stability of orbits in the libration region compared with 0.3 ≤ e 1 < 0.6. The innermost semi-major axis of stable orbits in the libration region is smaller than that of stable prograde orbits in the circulation region. This finding means that the critical semi-major axis of stable planets in polar orbits is smaller than that of coplanar prograde orbits, thus, raising the possibility of planet formation in the polar planetary discs. An important point to note is that the innermost stable orbits in the libration region appear when the inclination of orbit in the vicinity of 90 • according to the stability maps in Doolin & Blundell (2011) and Cuello & Giuppone (2019). Then, we will find out the innermost stable orbits for different binary parameters, the mass ratio of binary star λ and the binary eccentricity e 1 .
In order to obtain the innermost stable polar orbits, a large number of numerical simulations need to be carried out. Without loss of generality, the semi-major axis of the binary orbit is set as 1 AU, and the binary eccentricity e 1 in the range [0, 0.8] with scanned interval ∆e 1 = 0.05. Simultaneously, the mass ratio of binary star λ = [0.1, 1], where m 0 = 1M ⊙ and m 1 varied with an interval 0.05M ⊙ .
The semi-major axes of test particles range from 2 AU to 6 AU, with an interval of 0.01 AU. For a specific semi-major axis, eight inclinations were chosen from arcsin Eight mean anomalies distributed evenly in the range [0 • , 360 • ] for each inclination. The initial longitude of ascending nodes h 2 was set as 90 • in all cases. The orbits are circular initially, and the arguments of the pericenter are chosen arbitrarily. After 10 7 T 1 , the initial semi-major axes of the innermost orbits that remained stable are regarded as the critical stable boundary.
The following empirical formulas are computed using multivariable linear regression analysis of the minimum stable semi-major axis a c of test particles within 10 7 T 1 with different binary mass ratio λ and eccentricity, for 0.1 ≤ λ ≤ 1, 0.25 ≤ e 1 ≤ 0.6 with the coefficient of determination R 2 ≈ 0.8664. Several fitted curves were plotted, and the corresponding raw data are presented in Figure 2.
Apparently, the tendencies of minimum stable orbits with λ are different for several eccentricity intervals in Figure 2. The values of minimum stable semi-major axes of test particles moving in the polar orbits around the binary with e 1 = 0.2 are directly plotted in Figure 2 directly. Because the results with e 1 = 0.2 additional to the fitting processes of Equation 4 and 5 did not yield good fitting effects, which can be seen in For fixed binary eccentricity, the binary with comparable masses, λ = 1, is favourable to the stability of polar orbits. These orbits around the binary with mild eccentricity e 1 < 0.2 can stably exist in the location closer to the binary compared with those around moderate and highly eccentric  binaries. The stable boundary of polar orbits revolving the binary with high eccentricitye 1 > 0.6 change slightly for different binary mass ratios. The numerical simulation on the stability of polar orbits revealed that finding planets in polar orbits at a close location in the binary systems with comparable masses and mild eccentricities e 1 < 0.2 is more hopeful. In two other cases, the coplanar disc of the same binary stars (namely coplanar case) and the general disc around a single star (namely standard case) were simulated for comparison. "General disc" is based on the model of "minimum mass Solar Nebula", but with different scaling factors of solid surface density and solid enhancement beyond the ice line (f d and f i ce). The host in the standard case is a single star with a mass of 1.4 M ⊙ (the total mass of 99 Herculis).
Considering that the precession range of the orbital inclination in the polar disc is 10 • for the initial inclination distributed in the range 85 • ∼ 95 • , in order to make the three planetary discs comparable, the orbital inclination ranges of the planetesimals in the coplanar disc and the standard disc were set 0 • ∼ 5 • and initial longitudes of ascending nodes are chosen randomly.
In the numerical simulation, the orbits of planetesimals in planetary disks need to be as close as possible to the state of natural dynamic evolution, so as to reflect the real collision results. If we start from circular orbits initially, the planetesimals in polar disc need about 6×10 4 yr (about 10 3 T 1 , where T 1 is the binary period) to get a natural state from fixed initial longitudes of ascending node at 90 • , and need 9 × 10 5 yr (about 1.5 × 10 4 T 1 ) in the coplanar circumbinary disc from initial circular orbits.
Computing time costs are expensive for 30000 planetesimals. So, we make efforts to let planetesimals move initially naturally to the greatest extent on the bases of dynamical properties.
In the polar protoplanetary disc, initial inclination i is chosen in the range 85 • ∼ 95 • randomly.
In eccentric coplanar binary disc, considering the secular perturbation of binary, the initial eccentricities of planetesimals are chosen in the range [0, e pump ] arbitrarily (Moriwaki & Nakagawa 2004), The eccentricities of planetesimals are chosen from 0 ∼ 0.05 randomly in the general disc around the single star.
The inflated radii of planetesimals were adopted as dr = 2r The inflated radius is artificially set to control the collision speed accuracy. In pure Kepler motion of two-body problem, the velocity of the celestial body moving around the host satisfies v = GM r 1/2 . The differential of the formula can be obtained, dr = 2r 3 2 (GM ) −1 2 dν , which means that the velocity difference between two particles in Kepler motion at adjacent distance dr is dν . Therefore, the error between the relative velocity recorded here and the relative velocity of collision can be guaranteed by limiting the distance, dr, between two adjacent planetesimals before the collision. This error is proportional to the distance, that is to say, the relative velocity recorded when the distance between planetesimals is smaller will be closer to the relative velocity at the time of collision. In our code, we judge whether the two planetesimals are within the critical distance, dr. On the one hand, the collision requires the two planetesimals to be in close contact. On the other hand, we can control the error range between the relative speed calculated within the critical distance and the realistic collision speed.
During our simulations, the gas drag need not be considered. Because we found gas damping has little effect on the orbits of the planetesimals and collision results in our studies. On the one hand, for the outer part of the disc around 99 Her which locates far away from the binary, the local gas density is too small to dampen the planetesimals. The gas damping force is proportional to the gas density and the relative velocity between the gas and the planetesimal. According to the general model of gas disc (Sano et al. 2000;Umebayashi et al. 2013), at midplane, the gas drag force at 66 au is 2.8 × 10 −5 times than that at 1 au in "minimum mass Solar Nebula". On the other hand, for the inner part of the disc, the dynamics of planetesimals are strongly excited by the binary, leading to the damping effect of the gas disc working inefficiently. Rafikov (2013a) found that the gas drag does not resolve the fragmentation barrier issue in Kepler circumbinary systems. Because fast relative precession of planetesimal and binary orbits results in inefficient planetesimal apsidal alignment. Rafikov (2013b) also demonstrated that gas drag regulates eccentricity behaviour only for bodies with radii less than 1 km, which is below the adopted planetesimal size (1-10 km) in our article.

Outcomes of Collision
During the realistic collisions among planetesimals, the perturbation of gravitational force between the two planetesimals entering into the range of close encounters will gradually increase which can't be ignored. In our numerical simulations, gravitational interactions among planetesimals are neglected to save computation costs. We can use the physical parameters and the states of motion between the colliding planetesimals during entering into a close encounter to obtain the outcomes after collision approximately. The outcomes of collisions depend on the target mass M targ , the projectile mass M p , the target radius R targ , the projectile radius R p , impact velocity V i as well as impact factor b, which are defined as b = sin θ, θ is the acute angle between the line of impact velocity and the barycenter line. The impact velocity vector ⃗ V i is the relative velocity between the target and the projectile.
During our process to calculate the outcome of the collision between two bodies, we set the heavier one M targ , and the lighter one M p . Figure 4 summarized the regimes and the calculations on the mass of the largest remnant after each collision briefly. The detailed process and scaling laws used to calculate maps of collision outcomes can be found in Leinhardt & Stewart (2011);Stewart & Leinhardt (2012).
We have tried three methods to deal with the results of collisions. The first way, we assumed that the two planetesimals move in a straight line at a constant speed within the critical distance. A collision occurs when the distance between the planetesimals is less than the sum of the physical radii of the planetesimals. The statistical results show that almost no collisions have occurred, because this processing method ignores gravitational interactions between the planetesimals. The second way, all of the collisions are regarded as head-on collisions and the diversity of actual collision results is ignored. The third way considers the collision parameters and the calculation of the collision results is based on the coordinates and velocities of the planetesimals after they enter into the critical distance.
This processing method will increase the proportion of hit-and-run collisions, which remains the two planetesimals undamaged. However, we can eliminate the artificially enlarged effect of the hit-andrun regime by sufficient collisions. Considering that the growth of planetesimals requires multiple collisions, we made stochastic simulations to find out the final mass of the largest remnant for a specific planetesimal after successive collisions. As long as the cumulative number of collisions is sufficient, the artificially enlarged effect of hit-and-run will be eliminated, and will not affect the final mass of the largest remnant. This is also the reason why we made statistical simulations in our article.

Efficiency of accretion
The efficiency of accretion ξ is defined as follows, where, M t is the mass of the target, and M p is the mass of projectile. M lr is the largest mass of the remnant part after a collision between two planetesimals. The detailed calculations of M lr can be obtained in Section 4.2. Clearly, ξ = 1 means that the two planetesimals have merged together, thus demonstrating perfect accretion. 0 < ξ < 1 means collisions lead to partial accretion, and ξ < 0 implies target erosion.
The ratios of different efficiencies of accretion in the three discs are shown in Figure 4 and Table   2. There is one thing worth noticing, the peaks appearing in the vicinity of ξ = 0 include the overestimated part of hit-and-run. Hit-and-run can remain the planetesimals almost the same as before the colliding, and lead to ξ ≈ 0. The perfect accretion (ξ = 1) in the polar-ring case is 0%, while that in the coplanar and standard case is approximately 0.1261% and 0.04117%. Figure 4 shows that the accretion efficiencies among planetesimals that can cause accretion including part and perfect accretion in the polar protoplanetary disc are lower than either the coplanar binary disc or the single star disc. Meanwhile, the erosion collision (ξ < 0) occurrence in polar-ring cases reveals  Figure 4. The collision outcomes among planetesimals in the polar protoplanetary disc of 99 Herculis system, the coplanar disc of the same binary stars (coplanar case), and the general disc around single stars with the mass of 99 Herculis system (standard case) for comparison. The coplanar case has the same initial conditions as the polar case except for inclination randomly chosen from 0 • -5 • . The standard case has the same initial conditions as the coplanar case except the host is a single star with mass 1.4 M ⊙ , the total mass of 99 Herculis. slightly higher frequencies (50.81%) than those of the two other cases (49.15% and 44.17%). Because of the expensive computer costs, we only simulated one group for each case, so there are no error bars. If enough numerical simulations are conducted, the efficiencies of accretion in the polar and coplanar circumbinary disc may be very close, presumably within error bars. In general, collisions among the planetesimals in the polar protoplanetary disc are not favourable for accretion compared with the single-star disc.
According to the statistical data of collision results, we can basically deduce that there is a high probability that planetesimals cannot grow through collisions. In view of the following two reasons, we continued to carry out stochastic simulation: on the one hand, the proportion of hit-and-run results was overestimated in the process of collision result processing, and it was necessary to eliminate its influence through enough collisions; on the other hand, stochastic simulation can give out the maximum mass remained quantitatively after a series of collisions.
Before stochastic simulations of successive collisions based on the ratios of efficiencies of accretion, collisional timescales require research, which determines the number of collisions that may occur within the age of 99 Herculis.

Collisional timescale
During the integration time 10 5 yr, the numbers of collisions in the polar, coplanar and standard case are 7915,11894 and 2429 respectively as shown in Table 2. Collisions among planetesimals occur most frequently in the coplanar case for the strong dynamical disturbing from the binary. The number of collisions in the polar case is moderate for mild gravitational disturbances from the binary in the perpendicular plane which affect the inclination more than eccentricity. The average numbers of collisions in the polar and coplanar circumbinary discs are respectively three and five times as many as those in the standard disc surrounding the single star.
In single-star systems such as our solar system, the objects of the Kuiper belt are distributed at 30-50 AU from the Sun. The chance of collision among planetesimals will be remote considering the long period of objects and the shortage of disturbance. Generally, the collisional timescale of planetesimals with typical radius R p at semi-major axis a in the standard disc around the single star with a mass M A can be obtained by where, for "minimum mass Solar Nebula", f d = 1, and f ice = 4.2 beyond the ice line, which lies in about 2.7 AU for solar-type stars. According to , collisional timescales of planetesimals in coplanar circumbinary disc T C col and in polar circumbinary disc T P col can be estimated by The impact rates n S imp , n C imp and n P imp can be read out from our numerical simulations of three discs containing the same total numbers and physical radius distribution of planetesimals.
Although the impact rates of circumbinary planetesimals are several times than that in the standard disc, as shown in Figure 5, the collisional timescales T P col of planetesimals in the polar circumbinary disc, and T C col in the coplanar circumbinary disc with a moderate f d = 10 are longer than 10 7 yr which is beyond the observed dissipation timescale (< 6M yr) of gas disc (Haisch Jr et al. 2001;Wyatt 2008). Ribas et al. (2015) confirmed that the dissipation timescale of the gas disc is directly related to the stellar mass. The gas disc around high-mass stars (> 2M ⊙ ) dissipated up to two times earlier than low-mass ones. That means, gas giants are hardly formed in such distant locations through core accretion for the collisional timescales are too long to absorb enough gas. Next, stochastic simulations were conducted to find out whether solid cores of protoplanets can form.

Stochastic simulations
However, Batygin & Brown (2016) show that a possible distant giant planet with more than 10 Earth mass in the solar system can explain the cluster phenomenon of distant Kuiper Belt objects in the argument of perihelion and physical space. Kenyon & Bromley (2015) showed that super-Earth mass planets can form at 125-250 au around solar-type stars from swarms of 1 cm-10 m planetesimals within 1-3 Gyr in the annuli with a mass of approximately 15M ⊕ considering collisional damping.
However, this simulation did not include the gas drag, which is important for the vertical distribution and radial motion of dust and large bodies. The gas drag eventually becomes negligible once bodies of planetesimal size have formed.
After obtaining the probability distribution law which has been shown in Figure  1. The discrete distribution function of the accretion efficiency, F (ξ) = ξn<ξ p n , can be obtained by the data of Figure 4, for the distribution law as below.
4. For a ε, an accretion efficiency ξ = (ξ n−1 + ξ n )/2 will be chosed. The mass of the largest remnant part after a collision will be calculated by Equation 9. The remnant continues to participate in the next collision.
5. If ξ < −M p /M t leads to M lr < 0, which is a non-physical result, we need to regenerate a random number ε.
6. Sampling with probability proportional to size 100N c times, then the final mass of planetesimal after enough collisions will be obtained.
We carry out 10000 runs of probability proportional to size sampling for the collisions occurring in each kind of protoplanetary disc. The distribution of the final largest mass of the remnant part M lr /M ti are listed in Table 3. M ti is the initial mass of planetesimal before collisions. It is clear that the planetesimal can hardly grow bigger. In most cases, the final mass after 100N c successive collisions is around the original mass at the beginning of collisions. It may be a possible mechanism to produce dust continuously, making a debris disc observed in 99 Herculis. An older debris disc (> 100M yr) requires replenishment of dust through mutual collisions among a population of greaterthan-kilometer-sized planetesimals. Some of the youngest(about 10-Myr-old) debris discs may be the remnant of protoplanetary discs (Wyatt & Dent 2002). 99 Herculis is a main-sequence star with an age 6-10Gyr. So, the debris disc observed asks for a dust-generating process whether it is a long-lived debris disc surviving for Gyr timescales or a transient ring generated from a recent collision. Our Herculis will not make planetesimals grow, but produce dust steadily.
Since the molecular gas is collocated with the dust in the debris disk, a new semi-analytical equivalent of the numerical model proposed by Kral et al. (2016Kral et al. ( , 2017, which assumes CO is produced from volatile-rich solid bodies. P-type exoplanets in circumbinary coplanar orbits have been detected by Kepler (Martin & Triaud 2014). A majority of these planets are located at the boundary of stable zones over long timescales (Holman & Wiegert 1999). The stability limit is due to overlapping firstand second-order mean motion resonances with the binary, and is mainly influenced by the overlaps of three-body mean motion resonances in massive multi-planet systems (Wang et al. 2019). However, the formation of planets in the circumbinary coplanar disc is possible for a 2 > 20a 1 (Paardekooper et al. 2012;Marzari et al. 2012;Meschiari 2012), or a ≈ 10 − 20a 1 considering the gravity of an axisymmetric disc which can strongly suppress the eccentricities of planetesimals beyond and facilitate the easy growth of planetesimals (Rafikov 2013a). Moreover, the critical radial distance beyond which planetesimal accretion is possible increases with the rising binary eccentricities (Scholl et al. 2007). Thus, planet can hardly form around 65 − 130 au (3.94 − 7.88a 1 ) in the circumbinary coplanar Ying Wang et al.
protoplanetary disc with a 1 = 16.5 au and e 1 = 0.76 of 99 Herculis. By comparison, the accretion efficiency of planetesimals in the circumbinary polar protoplanetary disc is lower than that in the coplanar one according to our simulation. So, we can infer the formation of planets in the inner region or around the location where the current debris disc exist is difficult in the system of 99 Herculis.

Discussions
The innermost stable orbits simulated by Cuello & Giuppone (2019) and Chen et al. (2020) locate much inner (about 2.5a 1 ) for their integration time is 10 5 T 1 and 510 4 T 1 , while 10 7 T 1 in our article.
Cuello & Giuppone (2019) studied the evolution of misaligned circumbinary discs through hydrodynamical simulations. Viscous torques exerted by the binary make retrograde configurations easier to become polar than prograde circumbinary discs. Childs & Martin (2021) show that about five circumbinary planets form in polar and coplanar orbits in the vicinity of 5.4a 1 . However, it simulates the late stage of the formation of planets with purely gravitational interaction for the Moon-sized planetesimals and Mars-sized embryos. That means, every impact can be regarded as a perfectly inelastic collision, which leads to perfect accretion.
While the physical size of planetesimals in our simulations distribute in 1-10km. Various outcomes of collision including perfect accretion, partial accretion, hit-and-run, graze-and-merge and catastrophic disruption, super-catastrophic disruption, have been considered.
A possible third component which is about 2.4 times as faint as 99 Herculis B was reported three times (Kennedy et al. 2012). According to the three positions in sky coordinates, its possible stable orbits may be in Kozai cycles. The third component might be involved in the formation of the polar disk of 99 Herculis, which is interesting to further study. Lepp et al. (2023) find out that polar circumtriple orbits only exist within a critical radius, outside which circulating orbits precess about the binary angular momentum vector. Through smoothed particle hydrodynamics simulations by Ceppi et al. (2023), the wide range of disc inclinations in hierarchical systems with more than two stars may result from the secular oscillation of their orbital parameters. Chen et al. (2022) investigated the orbital dynamics of circumbinary planetary systems with two planets in polar orbits around the binary star. Under binary-planet and planet-planet gravitational interaction, the tilt angles of planets oscillate complicatedly.
The study on planetary formation in this paper is under the framework of the core accretion model. Some giant planets such as the four giants in HR 8799 system were detected far away from the central star through the method of imaging. If gravity instability gives rise to collapses in the solid component of the disc material, then giant planets can form in the outer regions of the protoplanetary disc. Planets formed beyond 100 au in solar-like gas discs through disc fragmentation can migrate inward, and produce giant protoplanets at a distance of a few tens of AU from the protostar via high-resolution numerical hydrodynamics simulations (Vorobyov & Elbakyan 2018).
Because typical collisional timescales in the polar disc are at least one order of magnitude higher than the dissipation timescale of the gas disk, gas giants are practically impossible to form by the core accretion model. In this model, the formation of planets mainly depends on the collisions among the planetesimals at the early stage. However, we couldn't exclude the possibilities of other models such as gravitational instability, pebble or dust-accretion producing planets or planet embryos.

Conclusions
The motion of planetesimals in polar orbits within the libration region of the outer-restricted threebody problem is simulated in this paper to study the stability of circumbinary polar orbits and the planetary formation in the circumbinary polar disc of 99 Herculis. Firstly, the empirical formulas of the stable critical semi-major axes of the polar orbits applied to 0.1 ≤ λ ≤ 1, 0.0 ≤ e 1 ≤ 0.15 and 0.65 ≤ e 1 ≤ 0.8 in Equation 4, as well as 0.1 ≤ λ ≤ 1, 0.3 ≤ e 1 ≤ 0.6 in Equation 5, are presented.
Secondly, the collision outcomes, colliding timescales and stochastic simulation of successive collisions among the planetesimals in the polar protoplanetary disc of the 99 Herculis system are statistically analyzed. The statistical results show that the collisions of planetesimals (with the physical radius 1-10 km and density 3g · cm −3 ) in the polar disc are the most unhelpful to the accretion and grow-up of planetesimals compared with the coplanar case and the standard case. Typical collisional timescales in the polar disc are at least one order of magnitude higher than the dissipation timescale of the gas disk. Furthermore, collisions of planetesimals in the polar protoplanetary disc of 99 Herculis will not make planetesimals grow, but produce dust steadily, which may explain the formation of the detected debris disc around 120 AU. Thirdly, considering the various outcomes of collisions among the planetesimals (1-10km), the performances of planetesimal growth by collisions in reference groups including the coplanar case and standard case are similar. The main differences between the three cases lie in the impact rate and the collision timescales.