Simulation of satellites and constellations for the assessment of collision avoidance operations

The recent rapid increase in the satellite population in low Earth orbit (LEO) causes the number of conjunctions between operational spacecraft to increase sharply as well. In order to be able to cope with the associated workload in the future, it is necessary to define precise and acceptable rules that determine which spacecraft has to evade and how the avoidance manoeuvres shall be performed. To enable the assessment of different conceivable rule sets, a holistic simulation framework has been developed within the “Rules4CREAM” (R4C) activity at the Technical University of Darmstadt. Based on current and expected trends in growth and distribution of satellites in LEO, long-term propagation in the order of years can be performed for arbitrary satellite populations with the R4C framework which enables the generation of realistic conjunction data. This paper presents the core modules of this simulation framework. The satellites are described by osculating Keplerian orbits and analytical models are implemented to represent the most dominant perturbations in LEO. Different modes of station keeping are depicted within the satellite operations module of the simulation. For constellations, a generic model has been developed that is able to model the simultaneous deployment of multiple constellation planes as well as other operational key aspects specific to constellations. The models are validated with other simulations and real-life data. Finally two methods for the all-on-all conjunction detection are examined and the number of detected conjunctions is compared.


Ω init
Initial RAAN to launch into Ω ref RAAN  Velocity of two satellites p p p , q q q Linearised orbit section of two satellites p 0 p 0 p 0 , q 0 q 0 q 0 Linearised position of two satellites at t 0 ū ū u , v v v Mean velocities of two satellites on linearised orbit section w w w Difference of linearised orbit sections of two satellites

Introduction
To avoid the generation of new space debris and therefore to allow safe and economic space operations in the future, the prevention of collisions in Earth's orbit is one of the key challenges. More than five thousand operational satellites are currently in orbit controlled by more than three hundred different operators [1]. Due to lower launch costs and new applications a record number of satellites has been launched in the past years, and the number of objects in low Earth orbit (LEO) will furthermore increase drastically in the near future, especially due to the deployment of multiple large satellite constellations. To enable safe, sustainable and economic operations, today's collision avoidance process (COLA) has to be evolved to take into account the increased number of potential collision events and for the fact that more events will involve two operational objects. Therefore, frameworks and systems are needed to enable coordinated and automated COLA processes, and unambiguous rules must be put into place that define which satellite has to perform a collision avoidance manoeuvre [2]. As part of the "Rules4CREAM" (R4C) activity, the concept of rule-based collision avoidance operations is simulated and assessed at the Technical University of Darmstadt in support of ESA's "Collision Risk Estimation and Automated Mitigation" (CREAM) activity [3]. The objectives of the R4C activity are the identification and modelling of a future low Earth orbit population, the development of different rulesets and their simulation-based application, a holistic analysis and assessment of benefits and disadvantages of individual rules and the derivation of requirements for future automated COLA systems.

Simulation of conjunction events
To enable the simulated application and the subsequent assessment of different rules and COLA operations, a representative set of conjunction events has to be generated by the simulation environment which is developed as part of the R4C activity. The simulation includes all operational satellites with a perigee below 2000 km, as events involving two manoeuvrable objects are the main focus of the research activity. Space debris objects can be included for additional studies. In addition to today's satellite population, potential future populations can be used to generate representative conjunction events. These have been generated by taking into account changing launch rates, satellite technologies and applications. The constellations used in the simulation are based upon extensive research of planned and proposed concepts. Perturbations and typical operations (depending on the satellite mission and type) are simulated over the course of the runtime and are the main focus of the simulation, as they are constantly changing the orbits of the objects and generate new and changing conjunctions. The goal of the simulation is not to generate and predict exact events between individual objects, but to generate a representative set of events, mapping a realistic approximation of the distribution of parameters. This includes for example the involved operators, the orbit types of involved satellites, mission types and the conjunction geometry. To ensure the generation of realistic results, the following features have been implemented: • Simulation of relevant perturbations on satellite orbits • Satellite and constellation operations including deployment, station-keeping and end of life operations • All-on-all conjunction detection with fixed screening thresholds • High performance long-term simulation with simulation time > 1 year in less than 24 h computing time

Simulation framework
When the simulation is started, a database containing all individual satellites and all constellations is loaded. This satellite and constellation related information includes among others the initial orbits, constellation parameters like the spacing between the individual planes, information about operators, the size and shape of the satellites, operational parameters like station-keeping modes and the age as well as the projected lifetime of the objects. Based upon these information the initial state of the population is generated. During the runtime, the orbits of all satellites are assumed to be constant for a certain period of time, typically 1 h to 24 h of simulation time. This simplification was chosen to increase performance, which is particularly necessary when several tens of thousands of objects of a possible future LEO population have to be considered, taking into account multiple constellations. During this time frame the position of each object on its orbit is propagated and the all-on-all conjunction detection is performed. Subsequently, the effects due to perturbations and satellite operations, that were relevant during the simulated time frame, are calculated for each individual object. They are combined and the orbits of all objects are changed simultaneously according to the calculated influence on the orbital parameters. An overview of the framework is shown in Fig. 1.
In the following chapters the individual sub-modules for the simulation of perturbations, satellite and constellation operations as well as conjunction detection are presented.

Single satellite model
The single satellite model is used to represent individual satellites, that are not part of constellations and models the perturbations acting on them, as well as the operations, focusing on station keeping. The structure of the model with its submodules is shown in Fig. 2.

Perturbations
The implemented perturbation models cover the most dominant perturbations in LEO. The implemented models are summarised in Table 1. Due to the large step size in the order of one hour to one day, only long periodic and secular perturbation effects are modelled.
The combined perturbations are verified and validated by comparing them with two-line element (TLE) data and ESA's Orbital SpaceCraft Active Removal (OSCAR) tool.

Operations
The effects of the residual atmosphere are highly noticeable in LEO and lead to a decreasing altitude of the satellites. Most satellites in LEO are therefore equipped with a station keeping system. In order to get an accurate and representative simulation, the station keeping has to be incorporated.
The operations module models the station keeping and potential end of life operations. For the station keeping, the three cases which are most prominent in LEO satellite operations are implemented: • Maintaining the ground track. • Maintaining the local altitude. • Keeping the orbit Sun-synchronous and maintaining the local time of ascending node (LTAN).
The purpose of the station keeping is not to optimise the manoeuvres and model the perfect station keeping, but to represent the general station keeping behaviour. Therefore, the implemented methods are highly simplified and focused on reliability and a short computation time.

Maintaining the ground track
Maintaining a ground track is especially important for Earth observation satellites. To mimic a realistic behaviour in the simulation this task is spilt into monitoring and controlling the ground track at the equator and at high latitudes. The ground track at the equator is mainly defined by the nodal    [4] period. The atmospheric drag decreases the semi-major axis (SMA) and therewith the nodal period. This causes the ground track of the satellite to drift east. In order to correct that drift, the SMA needs to be raised above the nominal SMA to cause a drift to the west. The simulation allows to define different bands of acceptable deviation from the nominal ground track for each satellite. To use the full allowed band, the manoeuvre has to be triggered when the ground track reaches the eastern boundary and has to be sized to stop the westward drift at the western limit. Vallado [4] derived a method to calculate the needed change in SMA relative to the nominal SMA. The needed change in SMA relative to the current SMA can be calculated by: Figure 3 shows the implemented station keeping of the ground track by highlighting the SMA correction manoeuvres and showing the deviation of the ground track at the equator. The deviation of the ground track at the highest latitudes is directly related to the inclination. If the ground track deviation exceeds the limit, an inclination correction manoeuvre is performed.

Maintaining the local altitude
The simulation implements two different kinds of altitude keeping. The simple altitude keeping maintains the SMA (4 cos 2 i − 1) ,

⏟⏞⏞⏞⏞⏟⏞⏞⏞⏞⏟
Δ a based on Vallado algorithm while keeping the eccentricity small. As Earth observation satellites need to maintain their local altitude, an advanced altitude keeping is implemented which regulates the SMA, eccentricity and argument of perigee. The station keeping is triggered if the SMA exceeds a threshold. The SMA is over corrected to the nominal SMA plus the current deviation. This change determines the size of the manoeuvre. The calculation of the manoeuvre size is derived from the Lagrange equations [4] with the assumption of small velocity changes: In order to keep the local altitude, the eccentricity and argument of perigee need to remain nearly constant. To naturally reduce the fluctuations, the station keeping keeps them close to the frozen values. The location of the manoeuvre is chosen so that the eccentricity is corrected back to the frozen value. There are two locations f 1 and f 2 which can achieve the necessary change in eccentricity Δe . One leads to a decrease, one to an increase of the argument of perigee. The location that brings the argument of perigee closer to the frozen value is chosen. The calculation of the maneuver location is also derived from the Lagrange equations [4]:  Figure 4 shows the implemented station keeping. The undermost graph show the deviation of the radius at points with a fixed argument of latitude. With the assumption of a rotation symmetric Earth and a nearly constant inclination this gives a good impression on how the altitude above a certain point on Earth changes.

Maintaining the LTAN
Sun-synchronous orbits maintain their orientation towards the Sun and keep their local time of the ascending node (LTAN) constant. A deviation in LTAN occurs when the nodal rate of the orbit, mainly caused by the J 2 effect, deviates from the average rate of the Sun's motion. The (9) implemented station keeping model tracks the current nodal rate, compares it to the rate of the Sun's motion and integrates the deviation over time. This yields the deviation from the LTAN. Once the deviation from the LTAN exceeds the limit, an inclination correction manoeuvre is triggered. If the altitude of the orbit remains constant due to other methods of station keeping, the nominal inclination at which the orbit is Sun-synchronous remains the same. If the altitude decreases, the nominal inclination at which the orbit is Sun-synchronous changes. The correction manoeuvre overcorrects the inclination to the nominal inclination plus the current deviation from the inclination. This causes the LTAN to drift in the opposite direction.

Constellation model
The characteristic of a satellite constellation is the interdependence of the individual satellites. In order to model this behaviour, constellations are modelled as a whole and not as individual satellites. The constellations are split into subconstellations that have the same SMA, eccentricity and inclination. A subconstellation consists of multiple orbital planes. Since the satellites within one plane follow the same orbit and are equally spaced, only one satellite is modelled to represent the entire plane. This increases the performance and decreases the runtime. The individual satellites are then generated within the conjunction detection.
To be able to model an arbitrary constellation, the model is designed to be as generic as possible. The constellation design chosen for the model is the Walker constellation. The orbital planes have the same SMA, eccentricity and inclination. The planes are identified by the right ascension of the ascending node (RAAN) relative to a reference plane. The Walker constellation design assumes circular orbits and evenly spaced constellation planes. It is equivalent to Streetof-Coverage constellations with equally spaced planes or to Flower constellations with circular orbits. Figure 5 shows how the constellation is implemented in the simulation and details the submodules.

Deployment and orbit raising
Constellation satellites are usually deployed at an altitude lower than their operational altitude and use their electric propulsion system to spiral up to their operational orbit. During this process they cut through a wide altitude band and can therefore lead to conjunctions with several different satellites, which is why this phase is of particular interest for this simulation. Since there is very limited technical information about the deployment process used by the different operators, TLE data of constellation deployment and orbit raising are analysed.

Analysis of TLEs
Starlink and OneWeb are two constellations that have already launched and populated several orbital planes, which is why they are chosen for the TLE analysis. They both deploy the satellites at an altitude lower than the operational and use the electric propulsion system to raise the orbit. For the initial population of the orbital planes, they follow different approaches. OneWeb typically populates one orbital plane per launch, which is why they have a continuous orbit raising. Starlink on the other hand populates three orbital planes with one launch. These three planes have a spacing in RAAN towards each other which is achieved by incorporating drift phases into the orbit raising. Figure 6 shows the SMA during this phase.
After deployment all Starlink satellites raise their orbit for approximately 10 days. After this initial phase, the satellites of the three orbital planes split up into three groups. The first keeps raising its orbit until it reaches the operational altitude. The others remain at their altitude and stay in a drift phase. Due to the SMA difference with respect to the operational SMA and the SMA of the first group, they experience a different nodal drift rate due to the J 2 effect. This causes the second and third plane to drift away from the first plane. Once the desired spacing is achieved, the second plane begins raising its orbit. The third plane remains in the drift phase until the required spacing towards the second plane is reached.

Design of a generic deployment model
The implemented model is based on the analysis of the TLE data and is designed as generic as possible to be able The orbital planes are characterised by their RAAN relative to the reference plane which is continuously drifting due to the J 2 effect. The magnitude of the drift is dependent on the inclination and SMA. All orbital planes have the same inclination and also the same SMA when they are in their operational orbit, so the RAAN relative to the reference plane remains constant. As discussed in the previous chapter the nodal drift during the orbit raising differs from the drift at the operational SMA. For the altitude range of the orbit raising the course of the nodal rate due to the J 2 effect over the time with constantly increasing SMA can be assumed as linear. Integrating the nodal rate over time yields the change in RAAN: The initial RAAN at which the plane has to be deployed is calculated by integrating the above expression over the orbit raise duration:  The change in SMA and RAAN during the orbit raising in one time step is calculated by: In the model, it is distinguished between launches that populate one plane per launch and launches with multiple planes per launch. If only one plane is populated, the orbit is continuously raised to the operational altitude. In the case of multiple planes per launch, the first plane raises its orbit continuously as well, while the other planes perform drift phases to change their RAAN relative to the first plane.
To keep the model as generic and as accurate as possible, the satellites are subjected to the complete perturbation model and the change in RAAN is monitored. The drift phase ends when the change in RAAN is equivalent to the difference in the final RAAN between the first and current plane. Once this point is reached, the orbit is raised.
The use of this approach enables the modelling of constellations with arbitrary numbers of planes, differences in RAAN between the planes and orbit raise speeds.

Validation of the model with Starlink TLE data
To validate the designed deployment and orbit raising model the simulated deployment is compared to the TLE data of the Starlink launch L1 of November 11th 2019. The comparison is shown in Fig. 7.
The aim of the simulation is to depict the general behaviour during the deployment and orbit raising in a realistic way, but not to represent it in every little detail to reduce the computation time. The model follows the real data closely and is sufficiently accurate for the desired use case.

Operations
The operation of constellations is comprised of the synchronisation of the orbital planes and the station keeping.

Orbital plane synchronisation
To avoid conjunctions between the satellites on different orbital planes the mean anomaly of each satellite has to be set to: where M ref is the mean anomaly of the reference plane's first satellite, F the phasing parameter and Ω rel the RAAN of the plane relative to the reference plane [7]. The different orbital planes experience almost the same perturbations, but the small differences lead to a shift in argument of perigee and mean anomaly. This leads to an invalid phasing and conjunctions between satellites of the same constellation.
To keep the synchronisation, the reference plane is simulated under the influence of the entire perturbations. The argument of perigee of the individual planes is then set to the same as the reference plane and the mean anomaly for correct phasing is calculated by Equ. 14 in every time step.

Station-keeping
Besides the synchronisation of the planes, there are two other main requirements for station keeping. The first is keeping the overall structure of the constellation by maintaining the desired spacing in RAAN. The second is keeping the desired altitude. To keep the simulation simple and runtime efficient, highly simplified station keeping models are implemented. The spacing in RAAN is done equivalently to the argument of perigee. The RAAN of the reference plane is propagated under the influence of the entire perturbation model. The RAAN of the individual planes is then set to the RAAN of the reference plane plus the desired spacing. Analysis of the Starlink TLE data showed that, with some exceptions, the satellites maintain their altitude within a band of ±120 m. This requires a high number of station keeping manoeuvres. Since Starlink uses electric propulsion systems, the manoeuvres can not be assumed to be impulsive. Due to the small variations in altitude and the difficulty of modelling the continuous manoeuvres the implemented model neglects the perturbing effects on the SMA and keeps it fixed for constellation satellites with electric propulsion.

Spare satellites
The high number of satellites in constellations make it inevitable that some of the satellites will fail. In order to assure continuous service, spare satellites are needed. The simulation is able to model in plane spare satellites, overpopulated planes and satellites in a parking orbit. The satellites in parking orbits have the same relative RAAN and follow the same deployment. They will be placed at a different altitude which can be specified.

End of life and replacement
Once the satellites of a constellation reach their intended end of life, an active de-orbit maneuver is simulated. In the simulation, this process is modelled by applying a constant decrease of the SMA on the satellites orbit, which acts in addition to the natural perturbations. Once the satellites reach an altitude lower than a user specified minimum altitude, they are removed from the simulation. The lifetime of the constellation satellites is usually pretty short, often in the order of five to eight years. This requires a frequent replacement of the satellites. If the constellation is set up by populating one plane per launch, the replacement is scheduled so that the new satellites reach the operational altitude when the old ones start to deorbit. In the case of multiple planes per launch the planes reach the operational altitude after different periods of time. This would lead to operational gaps or to overpopulated orbits. The simulation implements two different models. The first replaces the planes with multiple planes per launch as well. In this scenario the new satellites are launched when the old satellites re-enter. The second scenario uses individual launches for each plane to avoid the operational gap.

Conjunction detection
To perform all-on-all conjunction detection between the active satellite population a stochastic and a deterministic mode have been implemented. While the stochastic method has a high calculation speed, it can only reflect the quantitative and qualitative character of the conjunctions. The deterministic method, on the other hand, allows to determine the actual locations and times of closest approach. By implementing two methods, the obtained results can be compared and thus be validated.

Stochastic
The conjunction detection strategy for the stochastic mode has been derived from the orbit trace method. This method allows to calculate a conjunction probability for each combination of objects on Keplerian orbits using the equation where r 1 and r 2 are the bounding spheres of the two satellites and d min is the minimum orbit intersection distance (MOID) of the two orbits. Furthermore, v 2 is the velocity of object 2, 1 is the angle between the velocity vector v 2 of object 2 and the relative velocity vector v 1 − v 2 , and T 1 , T 2 are the orbital periods [8].

Selection of method
The orbit trace method has been chosen as it supports large time steps that are necessary to conduct the long-term investigations of this research in adequate computation time. Additionally, the method integrates well with the chosen propagation approach, as the conjunction probabilities are calculated based on Keplerian orbits. Constellation satellites can moreover be easily incorporated in this method as outlined in Sect. 4.1.3.

Procedure
To discard object pairs that cannot produce a conjunction efficiently before the application of the orbit trace method, several additional calculations are performed prior. The logic is shown in Fig. 8.
The apogee/perigee filter allows to sort out object pairs whose ranges of radius values of their respective orbits do not overlap [9].
If the subsequent constellation filter is toggled, all conjunctions between satellites belonging to the same constellation are excluded, assuming that the constellation operators are capable of synchronising the orbits of their satellites in order to prevent collisions. The orbit path filter discards candidate pairs based on the distance between their orbits [9]. For the object pairs that pass through the filter, it moreover delivers the distance d min between the orbits that is needed for calculating the conjunction probability in Equ. 16.
However, as the orbit path filter exhibits a singularity for coplanar orbits, the respective object pairs are identified by the coplanarity filter.
The distance between these objects with coplanar orbits is then calculated in the coplanar orbit intersection filter by solving a quadratic Equ. [10].
All the object pairs whose orbits have been determined to come close enough to produce a conjunction are then forwarded to the orbit trace method and the conjunctions are triggered based on the calculated conjunction probabilities of Equ. 16.

Incorporation of constellation satellites
To incorporate constellation satellites into the conjunction detection module, the conjunction probability is modelled as a Bernoulli process. For this purpose, the probability p of a conjunction to occur between the two orbits is calculated with Equ. 16, assuming there is only one satellite on each orbit. Let there be i satellites on the first orbit and j satellites on the second orbit. The amount of combinations between the two satellite groups then amounts to Assuming that all of the n conjunction combinations are equally likely to occur and independent, the probability P for k conjunctions to happen is Of interest is the probability P(k > 0) that at least one combination between the satellites i and j has a conjunction. The complementary probability P(k = 0) can be used to calculate The Bernoulli probability P(k > 0) is depicted in Fig. 9 for a probability of p = 1 ⋅ 10 −3 for a single event. As can be seen, the conjunction probability increases with the number of satellites on any of the two orbits and approaches 1 for a large number of satellites.

Deterministic
The deterministic conjunction detection strategy is based on the geometric filter chain as proposed by Woodburn and Dichmann [11]. In this procedure, candidate pairs are systematically excluded on the basis of their orbital geometry before determining the time and location of the closest approach for the remaining pairs.

Selection of method
In contrast to the (smart) sieve methods [12,13], which require time steps in the order of 5-15 min [13,14], the  Fig. 8 Flow chart of the implemented stochastic conjunction detection algorithm Fig. 9 Bernoulli probability P(k > 0) for at least one conjunction to occur for a probability of p = 1 ⋅ 10 −3 of a single event and for variable numbers of satellites i on the first orbit and j on the second orbit proposed geometric filter chain is more suitable for the objective of the R4C simulation as it is designed to examine large time periods in the order of years. The geometric filter chain integrates well with the chosen propagation approach, since it operates on Keplerian orbits. Assuming that the objects under investigation move along fixed Keplerian orbits during each time step, time steps of the order of hours can be chosen.

Procedure
As shown in Fig. 10, the initial logic that is applied for the deterministic method is the same as for the stochastic method. However, the condition of the coplanarity filter for considering two orbits coplanar is extended to also include orbit pairs that are at most separated by the chosen conjunction screening threshold at the most distant orbit points. The location, where conjunctions are possible can not be further reduced based on the orbit geometry for these near coplanar orbits. For this reason, they are directly forwarded to a sampling algorithm. This algorithm searches for the conjunctions between all of the near-coplanar object pairs by sampling the positions of the objects over the duration of the time step.
For the non-coplanar object pairs, in contrast, the location where conjunctions are generally possible can be narrowed down. A conjunction is only possible to occur within a true anomaly window around the line of mutual nodes of the two orbits based on the maximum distance from the relative node that a satellite could be and still be within the chosen screening threshold with respect to the other orbit. The subsequent time filter determines the time intervals of each object when it is within the respective true anomaly window. By searching for the overlaps between the time intervals of common residence of any object pair, the search intervals for the conjunctions can be further reduced. For any time window longer than 10 min, sampling is applied to search for conjunctions. In the conducted investigations more than 75 % of the candidate pairs have been on non-coplanar orbits with time windows shorter than 10 min. To reduce the computational load the orbits have been linearized within these time windows. This allows to approximate the conjunction range by solving a simple equation instead of performing computationally expensive iterations or sampling. The conducted investigations showed that all conjunctions of candidate pairs with a time window shorter than 10 min which have been detected by performing fine sampling of the objects could also be detected by applying the linearisation. The mean error remained in the order of centimeters for all conducted investigations.

Linearisation
To apply the linearisation, the Keplerian ellipses are approximated by line segments within the time windows. The velocity at the arithmetic mean of the temporal start and end point is used as the constant velocity for the interval With this approximation the positions of the satellites at time t have been modelled with the parametric equations where p 0 p 0 p 0 and q 0 q 0 q 0 are the respective positions of the satellites at the start of the interval t = t 0 and ū ū u and v v v are their respective mean velocities (see Fig. 11) [15]. The distance between the satellites at time t can thus be described by   [15]. Equation 22 has a minimum when D(t) = d(t) 2 has a minimum. For this reason, it is sufficient to calculate in order to save calculation time [15]. This function has a minimum when The time of closest approach (TCA) can be obtained by solving equation 24 for t: The conjunction range at TCA can finally be obtained by calculating

Incorporation of constellation satellites
For the non-coplanar combinations, the orbit path filter and the true anomaly window method can be applied to (23) constellation satellites in the same way as for single satellites by using one representative satellite of each constellation plane as reference satellite, since both algorithms operate on the orbit geometry and not on the individual satellites. The other satellites of a constellation plane are first considered when the combinations are processed by the time filter. The algorithm replicates the time windows according to the number of satellites on the respective orbit. This process is depicted in Fig. 12 for the case of a single satellite and a constellation plane with four satellites being processed. The assumption is made that the satellites are equally spaced along the constellation plane within the time domain, so that the time windows can be distributed equally. For the coplanar combinations, sampling is performed for all possible combinations of constellation satellites.

Comparison of methods
In order to ascertain that the deterministic and the stochastic methods produce consistent results, a set of simulation runs is carried out and the obtained conjunction data is compared. The chosen settings for all subsequent simulation runs are listed in Table 2.

Influence of time step
In the first investigation the influence of the time step length on the quantity of the detected conjunctions is examined. The three time step lengths 6 h, 3 h and 1 h are studied for both the stochastic and the deterministic method by running 16 Monte Carlo runs respectively, where random initial mean anomalies are chosen for the satellites for each run. The distribution of the number of detected conjunction events is depicted in Fig. 13.
The median of the results of the stochastic conjunction detection method (marked orange) consistently amounts to around 187,600 for all time steps and agree well with the results of the deterministic method (marked blue). The stochastic results slightly overshoot the amount of conjunctions  in comparison with those of the deterministic method for a time step of 6 h and 3 h and slightly undershoot it for a time step of 1 h. The relative deviation of the average amount of detected conjunctions between the deterministic and the stochastic method is below 2 % for all time steps.
The greater expansion of the kernel density of the deterministic method in Fig. 13 shows that the conjunction number has a greater variance for the deterministic method than for the stochastic mode. This is most likely attributed to the influence of the randomised initial mean anomaly of the satellites. If a satellite combination on synchronised orbits starts with anomalies that cause them to have a conjunction, the conjunction will be detected by the deterministic method repeatedly for every orbit revolution. With the stochastic method, on the other hand, this effect is averaged out by the orbit trace method, where the anomaly is not taken into account. In the same way, the downwards deviations can be reasoned. This observation shows the importance of the anomalies in the initial database for the deterministic method.

Influence of screening threshold
To study the influence of a variation of the screening threshold on the results, another set of Monte Carlo runs is performed, where the time step length is fixed to 3 h and the screening threshold is varied (0.1 km, 1.0 km, 5.0 km, 10.0 km; compare Table 2).
The results shown in Fig. 14 exhibit similar average amounts of detected conjunctions for each threshold. The quartiles of the results for the two conjunction detection show overlapping ranges for all thresholds. As in the previous investigation, it can be observed that the deterministic results scatter more than the results of the stochastic method.

Distribution of conjunction range
The number of detected conjunctions is not the only figure of merit that has to be considered in order to evaluate the quality of the results. In order to produce realistic conjunction data, the simulation shall also yield a realistic distribution of the conjunction range corresponding to the conjunction events. This distribution is shown for the stochastic and deterministic method in Fig. 15. To consider the influence of the time step length, the simulation has been conducted for a time step of 6 h, 3 h and 1 h.
It can be seen, that the density of the deterministic method exhibits few conjunctions with small conjunction ranges and more conjunctions the larger the conjunction range becomes. In contrast, the distribution of the stochastic method deviates very strongly, peaking at low conjunction ranges. The reason for this deviation is how the conjunction range is determined for the deterministic methods and the stochastic method. In the deterministic method, the conjunction range is calculated by applying the linearisation as explained in Sect. 4.2.2.
In the stochastic approach, there is no way of determining the conjunction range, as the anomalies of the satellites are not considered by the orbit trace method. For this reason, the distance between the orbits d min , previously determined by the orbit path filter and the coplanar orbit intersection filter, is used to approximate the conjunction range and the exact points along the line of mutual nodes are assumed to be the points of closest approach. As can be seen in Fig. 15, this does not represent the expected distribution as predicted by the deterministic method.

Influence of constellation satellites
All of the previous investigations have been conducted using a satellite population database which contained both the operative single satellites and operative constellation satellites. This causes all satellites to be treated as single satellites by the simulation framework internally, whereby the developed constellation model is not used. For validating the constellation model, the constellation satellites are filtered out of the single satellite population database, leaving 1,100 operative single satellites. Subsequently, artificially created constellations are added to a separate constellation database. To ensure that the created generic constellation satellites cause a significant number of conjunctions with the single satellite population, the height distribution of their orbits is analysed. For this purpose, the distribution of the SMA of all single satellites is shown in Fig. 16. An accumulation of the SMA in the interval from 6800 to 7040 km can be observed.
Five artificial constellations are created that are equally distributed with respect to their SMA within the interval from 6800 to 7040 km in 60 km increments (red lines in Figure 16). Each of these constellations consists of 8 orbit planes, equally distributed with respect to their RAAN, where all orbits have an inclination of 52 • and an eccentricity of 0.0. The setup is depicted in Figure 17 for the first (lowest altitude) and fifth (highest altitude) generated constellation. The other three constellations show the same setup, but are not shown for reasons of clarity of the illustration. As can be seen, the orbital planes of the individual constellations are stacked above each other forming a fence in the critical SMA range. An increased amount of conjunctions is thus expected to occur along these planes.
Four sets of simulations are conducted that differ in the number of satellites placed on the constellation orbit planes. 1, 2, 10 and 100 satellites are chosen, which lead to a total number of constellation satellites of 40, 80, 400 and 4,000 respectively. Each of the four sets of simulations with the varying satellites on the constellation planes is conducted 16 times to ensure statistical significance of the results. The number of detected conjunctions is shown in Fig. 18. For one constellation satellite per orbit plane, the same characteristics can be observed that have been discussed in Sect. 4.3.1, i.e. the stochastic method is slightly overestimating the number of conjunctions showing a smaller variance of the results than the deterministic method. Increasing the number of constellation satellites per plane to 2 yields an increase of about 1000 average detected conjunction for both methods. Two aspects can be observed in the other two plots, which show the detected conjunctions for 10 and 100 constellation satellites per plane respectively. Firstly, the conjunction numbers that are determined for the two methods slowly diverge as the number of satellites per plane increases. The stochastic method overestimates the conjunctions for 1 and 2 satellites per plane by around 2 % with respect to the deterministic result. For 10 satellites per plane this relative difference amounts to 4 % and for 100 satellites per plane the relative deviation reaches 8 %. Secondly, the variance of the stochastic method is increasing for increasing numbers of constellation satellites, reaching a similar variance as the one of the deterministic method for 100 satellites per plane. The results in the first row of Fig. 18 show that modelling the conjunction probabilities as a Bernoulli process is able to correctly integrate the conjunction probabilities for single satellites, since it provides a continuous transition from the single satellite run (upper left diagram) to constellations with few satellites on their orbit planes (upper right diagram). For constellations with up to 10 satellites on each plane the results show good agreement with predicted conjunction numbers in the same order. However, an increasing overestimation of the number of conjunctions can be observed for the stochastic method for increasing numbers of constellation satellites per orbital plane. This could indicate that the Bernoulli process assumption of stochastically independent events might not be justified for modelling the conjunction probability of very large constellations. Other stochastic models shall therefore be analysed in future studies to see if they can provide a better approximation also for larger constellations.

Choice of simulation parameters
The propagation module and the conjunction detection module of the presented simulation framework have been designed and validated for time steps between 1 h and 12 h to allow for long-term simulations while still reflecting the important characteristics of satellite operations such as station keeping in the created results. A time step of 6 h has shown a good trade-off between accuracy and run time and is the recommended standard time step length. Furthermore, it is recommended to activate all implemented perturbation models for the best accuracy of the simulation results created. It has been found that the perturbation models do not significantly deteriorate the simulation run time in comparison with the conjunction detection module which accounts for the largest share of the calculation time. While in theory, arbitrarily high screening thresholds can be chosen for the simulation runs, it is recommended to choose screening thresholds below 20 km. For larger screening thresholds the computation time increases sharply for the deterministic conjunction detection method, because many satellite combinations pass through the filter chain. The stochastic method is capable of handling larger screening thresholds, albeit the accuracy of the predicted amount of conjunctions has not yet been tested for very large values. The deterministic method is the recommended standard conjunction detection method as it enables the creation of realistic simulated conjunction events with tangible satellite positions at TCA, which can be used for further processing.

Conclusion
A large-scale orbit simulation has been presented in this paper that was developed with the objective of generating representative list of conjunctions for arbitrary object  populations in LEO to allow the assessment of different rule sets for collision avoidance operations. The most dominant long periodic and secular perturbations in LEO are considered for the propagation of the satellite population and validated by comparison with historic TLE data and ESA's OSCAR tool. In order to realistically depict the dynamics of the LEO environment of active satellites, station keeping, end of life operations and a generic deployment model for constellations have been developed. The deployment process has been validated with TLE data of the Starlink launch L1. For the detection of conjunctions, a deterministic and a stochastic method have been implemented that are both able to process constellation satellites, leveraging the movement of satellites along the same orbit to increase performance. A comparison of the detected number of conjunction events for different simulation settings has shown that both methods produce similar results. While the stochastic method requires simpler computational operations to be performed, it has the disadvantage that the distribution of the conjunction range cannot be represented correctly.