Nadine Friesen, Tim Sander, Karl Nachtigall, Nils Nießen 11institutetext:
Institute of Transport Science, RWTH Aachen, Aachen, Germany
11email<EMAIL_ADDRESS>
WWW home page: https://www.via.rwth-aachen.de/ 22institutetext: Chair of
Traffic Flow Sciences, TU Dresden, Dresden, Germany 33institutetext: Lehr- und
Forschungsgebiet Kombinatorische Optimierung, RWTH Aachen Aachen, Germany

# The complexity of the Timetable-Based Railway Network Design Problem

Nadine Friesen 11 Tim Sander 22 Christina Büsing 33 Karl Nachtigall 22 Nils
Nießen 11

###### Abstract

Because of the long planning periods and their long life cycle, railway
infrastructure has to be outlined long ahead. At the present, the
infrastructure is designed while only little about the intended operation is
known. Hence, the timetable and the operation are adjusted to the
infrastructure. Since space, time and money for extension measures of railway
infrastructure are limited, each modification has to be done carefully and
long lasting and should be appropriate for the future unknown demand.

To take this into account, we present the robust network design problem for
railway infrastructure under capacity constraints and uncertain timetables.
Here, we plan the required expansion measures for an uncertain long-term
timetable. We show that this problem is NP-hard even when restricted to
bipartite graphs and very simple timetables and present easier solvable
special cases.

This problem corresponds to the fixed-charge network design problem where the
expansion costs are minimized such that the timetable is conductible. We model
this problem by an integer linear program using time expanded networks.

To incorporate the uncertainty of the future timetable, we use a scenario-
based approach. We define scenarios with individual departure and arrival
times and optional trains. The network is then optimized such that a given
percentage of the scenarios can be operated while minimizing the expansion
costs and potential penalty costs for not scheduled optional trains.

###### keywords:

Network Design, Timetabling, Robust Optimization, Railway Planning, Railway
Network Design, Strategic Timetabling

##  1 Introduction

Traditionally, the planning process for public transportation and railways is
executed sequentially: infrastructure planning, line planning and timetable
planning. Network design is included in the infrastructure while timetabling
is usually executed later in the planning horizon. Because of that, the
infrastructure is often assumed to be fixed during the tactical timetabling.
This leads to strict constraints for the timetable and reduces the options to
adjust the timetable to the transportation demand which possibly changes
between the infrastructure planning and the actual operation. Therefore more
and more railway companies in Europe shift from this traditional approach to a
timetable-based one by including a long-term timetable into infrastructure
planning. This long-term timetable is then used as input for further planning
steps, including network design. This approach is used in several western
European countries, including Switzerland, Germany and the Netherlands.
However, it is rarely covered in the scientific literature. The ideas and the
conception of strategic timetables in Germany have been described by Weigand
et Heppe in [11]. Only recently, Polinder et al proposed a mathematical
program for strategic timetabling in [19].

Network design problems describe the decision which nodes or lines of a
network should be expanded or build to meet a given demand. They occur in many
different contexts for example in transportation, communication or
electricity. The classic network design problem is described in [14] by
Leblanc. The problem is modelled as a nonlinear mixed-integer problem. An
overview over models and algorithms for transportation network design is given
in [16]. A more recent review by Chen et al. can be found in [1]. The fixed
charge capacited network design problem with multi-commodity flow is closely
related to the railway network design problem. If an arc is used, a fixed
charge is applied, independently of the amount of flow that is transported
over this arc. A survey of this problem is given in [7]. They also provide
possible relaxations of the problem. The capacited network design problem for
a multi-commodity flow describes a network design problem where the capacity
of each arc is limited and the flow starts from different sources and should
arrive at different sinks. In [3] this problem is described and strategies to
solve this problem with Benders decomposition are analyzed.

In the context of network design for railway infrastructure in the traditional
planning approach, vague demand prognoses are assumed. Based on these, the
necessary capacity expansions and new constructions are determined. The
network design problem in railways without a timetable or temporal constraints
has been studied for example in [24] where Spönemann provides a MIP
formulation for the problem as well as some computational results. In [13]
Kuby et al. give an example of the railway network design where extensions of
the network are determined by using a mixed-integer linear program. In a
second step, they define different stages of the extensions through a
heuristic backwards time sequencing procedure.

Most network design problems do not consider a temporal aspect and assume a
static demand that does not vary over time. This can be assumed if a worst-
case demand over all time steps for the network can be given. However, if such
a worst-case cannot be determined, the timetable needs to be included in the
network. As the problem becomes computationally intensive if we additionally
include a temporal dimension, there exist only very few approaches that
consider a temporal aspect. In [10] Guihaire et al. give a survey over transit
network design. They include the planning of the transit routes network, the
frequencies and departure times which covers problem related to both network
design and timetabling. Zhao et al. combine in [28] route network design,
vehicle headway and timetable assignment and provide a metaheuristic search
scheme consisting of simulated annealing, tabu and greedy search methods to
solve the problem.

Timetables in the railway network design problem have been studied in [21].
Schöbel at al. provide a formalization of the network design problem under
timetable constraints as well as sketch some possible algorithms and
heuristics. No implementations or computational results are given. In [22],
the approach is extend with some considerations to reduce the computation time
of the problem. In [8], a neighborhood search algorithm is given where as in
[9] a meta-heuristic approach is developed.

Due to the long planning horizons for railway infrastructure, the strategic
timetable is subject to uncertainty. One method to handle this uncertainty is
robust optimization. In the context of railways, robust optimization has been
extensively studied. A survey about robust optimization in railway planning
can be found in [15]. However, most research considering robustness in
railways focuses on robust timetables and not on robust infrastructure.
Results about robust infrastructure in the case of failure which occur in the
operation can be found in [5].

Robustness in network design problems can be found in different contexts:
Ukkusuri et al. [26] consider the network design problem under demand
uncertainties for transportation networks. The uncertainty in the demand is
modelled through random variables with known distributions for the entries of
the origin-destination-matrices. Pishvaee et al. ([18]) describe a
probabilistic programming approach for the supply network design problem under
uncertainties. In [1] Chen and al. present a bi-objective-reliable network
design problem model that optimizes the reliability of the capacity and the
travel time under demand uncertainty. A simulation-based multi-objective
genetic algorithm solution procedure is developed. In [17] Mudchanatongsuk et
al. present a robust optimization-based formulation for the network design
problem under transportation cost and demand uncertainty. An approximation to
this is shown to be done efficiently for a network with single origin and
destination per commodity and general uncertainty in transportation costs and
demand that are independent of each other. For a network with path
constraints, an efficient column generation procedure to solve the linear
programming relaxation is given. These approaches do not take a temporal
aspect into account.

Up to the authors’ knowledge, little research considering the network design
problem under both uncertain demand and timetable constraints has been
conducted yet. In [20], we integrate the timetable conditions by modelling the
arrival and departure times as variables and provide computational results for
a scenario-based robust approach. In this paper, we aim to provide another
model for the same railway network design problem under timetable constraints.
A comparison of both models can be found in [4] where we describe the
theoretical differences between the models and compare their computational
results for a test instance of a network around Dresden.

The model presented here is based on a time expanded network. This considers
the influence of timetables as well as the influence of uncertainty on the
railway network design problem. We propose an optimization model that allows
for the expansion cost-optimal calculation of railway networks while the
demand is given as a set of trains with a respective starting and end node as
well as earliest departure and latest arrival times. We prove that this
problem is NP-hard to solve and present some graph classes for which the
problem easier to solve. Additionally, we present a formalization of the
optimization problem for the deterministic and the robust case.

This paper is structured as follows: First, we define the problem and
introduce the notation in section 2. In section 3 we describe our mathematical
model for the network design problem which we adapt for the robust approach in
section 3.1. Then, we provide the proof of the NP-hardness of the problem in 4
as well as some easier special cases. Some computational results are given in
section 3.2. We conclude this paper in section 5 with some final remarks and a
short outlook onto further research.

##  2 Problem definition and notation

In this paper, we expand the classic network design problem by incorporating
timetable constraints and adapting railway specific capacity measures. We aim
to find an expansion cost minimal railway network to operate a given long-term
timetable on such that the capacity constraints are respected.

The problem is closely related the well-known fixed-charge network design
problem as described in [25]. Instead of costs for each unit of flow as in the
more common minimum cost flow problem, we consider costs that occur if an arc
is used at all, independently of the amount of flow on this arc. Here, these
fixed charges correspond to the expansion costs for arcs. If the capacity of
an arc is not sufficient, the arc is expanded and the expansion costs incur.
We consider a setting in which it is a binary decision to expand an arc, so
the arc is either expanded fully or not at all.

We consider a railway network consisting of stations represented by nodes and
lines represented by arcs $(i,j)$. Each arc has a travelling time $tt_{ij}$
which is the number of time steps, a vehicle needs to pass that arc.

Since the capacity of railway lines is limited, it is necessary to include a
capacity measure to estimate the number of tracks needed for the operation of
the input timetable. Railway line capacity depends on several different
factors, most notably the type and amount of trains running on a line, the
number of tracks and the train control system in place. Here, we assume the
capacity to be given for each line. The capacity of each arc is given by
$c_{ij}$ which is the maximum number of vehicles that can start in node $i$ to
node $j$ in a predefined time interval $\widetilde{t}$. The capacity of a line
can be extended by $\widetilde{c_{ij}}$ for the cost of $k_{ij}$. The lines
that can be built, expanded or already exist are predefined as well as the
corresponding costs. This is due to the fact that not all lines could be built
due to financial and spacial restrictions.

To integrate the temporal aspect we use a time-expanded network. That is, for
each node $i$ we have as many copies $i_{t}$ as we have time steps $t\in T$
Where $T$ is the set of all time steps. The edges of the time expanded network
rely on the edges of the original network. That is for each edge $ij$ in the
original network we define edges $i_{t}j_{t^{\prime}}$ in the time-expanded
network where $t^{\prime}-t$ is the travelling time for edge $ij$. To describe
these edges we use a adjacency matrix $A$ where $a_{i_{t}j_{t^{\prime}}}=1$ if
the edge $ij$ exists and the travelling time between $i$ and $j$ is
$t^{\prime}-t$ and $t,t^{\prime}\in T$. As we only model a fixed travel time
for all vehicles in this paper, we omit $t^{\prime}$ for simplicity in the
following.

Further, we have a long-term timetable as an input which is described through
a set of trains $V$ where each train $v\in V$ has a departure and an arrival
node $d_{v}$ and $a_{v}$ with an earliest departure time $t_{d_{v}}$ and a
latest arrival $t_{a_{v}}$ time at these nodes. These trains then have to be
routed such that the, potentially expanded, capacity is not exceeded.

For each edge and each pair of vehicles, there exists a minimum headway time
$M_{v_{1}v_{2}ij}$ which states the minimal time train $v_{2}$ can leave onto
line $ij$ after $v_{1}$ left station $i$.

We define binary variables $b_{ij}$ that state if the edge $ij$ is expanded.
The binary variables $x_{v,i,j,t}$ describe if the edge $ij\in E$ is used by
vehicle $v\in V$ which leaves node $i\in N$ at time point $t\in T$.

Additionally, we introduce VIA-nodes which are stations a train has to pass.
For VIA-Nodes, we define a pair $(v,n)$ where train $v\in V$ needs to pass
node $n\in N$ that is $\exists i\in N,t\in T:x_{i,n,v,t}=1$. Furthermore,
connections are specified through a station and for the connection a pair of
trains. In these stations, both trains have to stop and allow the transition
of passengers from the first train to the second one. Connections are given as
a tuple $(n,v_{1}.v_{2})$ where train $v_{1}$ should arrive at node $n$ before
train $v_{2}$ to ensure a connection.

##  3 Mixed Integer Linear Programming Formulation

In this section, we want to formalize the network design problem for railways
under timetable constraints. We model this optimization problem through a
time-expanded network. The objective of the network design problem for railway
infrastructure is to minimize the expansion costs of a railway network while
the constrains ensure that the timetable is conductible.

The optimization model can informally be described as follows:

$\displaystyle\underset{b_{ij}\text{ }binary,ij\in E}{\text{minimize}}$
$\displaystyle\sum_{ij\in E}ExpansionCost_{i,j}\cdot b_{ij}$ such that
(Expanded) capacity is not exceeded Earliest departure times are respected
Latest arrival times are respected Minimum headway times are respected Flows
are conserved Connections are respected VIA-Nodes are respected

With the notations from Section 2 we obtain the following mathematical
formalization.

#### Objective function

$\displaystyle\min\sum_{ij\in E}k_{i,j}\cdot b_{i,j}$ (1)

The objective function is to minimize the expansion costs. For that, the
variable $b_{ij}$ for each arc ${ij}$ becomes true if the capacity of this arc
is expanded. The costs for the expansion of that track $ij$ are given by
$k_{ij}$. To obtain the costs for the whole network, we build the sum of these
costs over all arcs.

#### Capacity

$\displaystyle\sum_{v\in V,t\in[t_{0},t_{0}+\widetilde{t}]}x_{i,j,t,v}\leq
c_{ij}+b_{ij}\cdot\widetilde{c_{ij}}$ $\displaystyle\forall ij\in E,t_{0}\in
T$

We include the capacity $c_{ij}$ as the number of trains on a line which
cannot be exceeded per a given time interval $\widetilde{t}$. The binary
variables $x_{i,j,v,t}$ are true if train $v$ leaves node $i$ at time $t$ to
get to node $j$. Hence, the sum of these $x_{i,j,v,t}$ needs to be less or
equal to the capacity $c_{ij}$ for a predefined time interval
$\widetilde{t}\in\\{1,\max\\{T\\}\\}$. This needs to be respected for each
time interval $[t_{0},t_{0}+\widetilde{t}]$ for all time points
$t_{0}\in\\{0,\dots,\max\\{T\\}-\widetilde{t}\\}$. This capacity $c_{ij}$ can
be extended by $\widetilde{c_{ij}}$ to allow more trains per time window. If
the line is extended, $b_{ij}$ has to be true and therefore the costs of the
extension have to be respected in the objective function.

#### Departure and arrival

$\displaystyle\sum_{t<t_{d_{v}}}x_{d_{v},j,t,v}=0$ $\displaystyle\forall v\in
V,\forall j\in N$ $\displaystyle\sum_{t\geq
t_{d_{v}}}a_{d_{v_{t}},j_{t^{\prime}}}\cdot x_{d_{v},j,t,v}\geq 1$
$\displaystyle\forall v\in V,\forall j\in N$
$\displaystyle\sum_{t>t_{a_{v}}}x_{i,a_{v},t,v}=0$ $\displaystyle\forall v\in
V,\forall i\in N$ $\displaystyle\sum_{t\leq
t_{a_{v}}}a_{i_{t^{\prime}},a_{v_{t}}}\cdot x_{i,a_{v},t-tt(i,a_{v}),v}\geq 1$
$\displaystyle\forall v\in V,\forall i\in N$

These constraints ensure that no vehicle $v$ leaves before their earliest
departure time $t_{d_{v}}$ which is no $x_{i,j,t,v}$ is $1$ for $t<t_{d_{v}}$.
Furthermore, they make sure that the departure from the starting node $d_{v}$
only occurs along an existing line. This is not ensured by the flow
constraints as the starting and end nodes are not included in the flow
constraints we use in this model. These constraints are implemented for the
arrival node $a_{v}$ and latest arrival time $t_{a_{v}}$ accordingly.

#### Minimum Headway time

$\displaystyle
x_{i,j,t_{1},v_{1}}\cdot(M_{ijv_{1}v_{2}}-(t_{2}-t_{1}))-(1-x_{i,j,t_{2},v_{2}})\cdot
max\\{0;M_{ijv_{1}v_{2}}-(t_{2}-t_{1})\\}\leq 0$ $\displaystyle\forall ij\in
E,t_{1}<t_{2}\in T,v_{1}\neq v_{2}\in V$ (4)

Furthermore, we implement minimum headway times. These are defined as the
amount of time that has to pass between the departure of two trains following
each other on the same track. They depend on the speed and the acceleration of
the trains, the blocking distance and the distance to the next station which
permits a change of the train sequence. Therefore, the minimum headway time is
individually determined for each line $ij$ and each pair of trains
$(v_{1},v_{2})$.

The most intuitive way to implement the minimum headway time is to multiply
the difference between the departure times and the minimum headway time with
the binary variables $x_{i,j,t_{1},v_{1}}$ and $x_{i,j,t_{2},v_{2}}$ for all
$ij\in E,t_{1}<t_{2}\in T,v_{1}\neq v_{2}\in V$ which indicate that the trains
are driving on that line at the indicated time. By that, we obtain the
constraint

$x_{i,j,t_{1},v_{1}}\cdot(M_{i,j,v_{1},v_{2}}-(t_{2}-t_{1}))\cdot
x_{i,j,t_{2},v_{2}}\leq 0$ (5)

To avoid this quadratic constraint, we use the linear constraint given in 3.
If both trains $v_{1}$ and $v_{2}$ use the line $(i,j)$, we obtain
$(M_{ijv_{1}v_{2}}-(t_{2}-t_{1}))\leq 0$ which is exactly the constraint we
are aiming for. If $x_{i,j,v_{1},t_{1}}=1$ and $x_{i,j,t_{2},v_{2}}=0$, we
obtain
$(M_{ijv_{1}v_{2}}-(t_{2}-t_{1}))-max\\{0;(M_{ijv_{1}v_{2}}-(t_{2}-t_{1}))\\}\leq
0$ which is true trivially. If $x_{i,j,v_{1},t_{1}}=0$, equation 3 becomes
$(1-x_{i,j,t_{2},v_{2}})\cdot-max\\{0;M_{ijv_{1}v_{2}}-(t_{2}-t_{1})\\}\leq 0$
which is $0\leq 0$ for $x_{i,j,t_{2},v_{2}}=1$ and
$-max\\{0;M_{ijv_{1}v_{2}}-(t_{2}-t_{1})\\}\leq 0$ for
$x_{i,j,t_{2},v_{2}}=0$. Both are easy to verify.

#### Flow constraints

$\displaystyle\sum_{i\in N,t\in T,j\neq a_{v}}a_{i_{t},j_{t^{\prime}}}\cdot
x_{i,j,t,v}-\sum_{i\in N,t^{\prime\prime}\in T,j\neq
d_{v}}a_{j_{t^{\prime}},i_{t^{\prime\prime}}}\cdot x_{j,i,t^{\prime},v}=0$
$\displaystyle\forall v\in V,j\in N,t^{\prime}\in T$ (6)

The train paths have to be consistent. That is that no trains appear or
disappear somewhere else than their origin or destination node. This is made
sure through the inclusion of flow conservation constraints. All trains that
enter a time-space-node have to leave that node through an existing line
except for their arrival node. By including the parameter
$a_{i_{t},j_{t^{\prime}}}$ which is true if the connection is possible we make
sure that only existing lines are used.

#### Connections

Connections are given as a tuple of two trains $v_{1}$ and $v_{2}$ and a
station $i$ where the connection should occur. By the constraints

$\displaystyle\sum_{n=1}^{|N|}a_{n_{t},j_{t}^{\prime}}\cdot
x_{n,i,t-tt_{n,i},v_{1}}-\sum_{n=1}^{|N|}a_{i_{t},n_{t^{\prime}}}\cdot
x_{i,n,t,v_{2}}\geq 0$ $\displaystyle\forall t\in T,(v_{1},v_{2},i)\in
Connections$ (7)

we ensure that $v_{1}$ arrives before or at the same time as $v_{2}$.
Furthermore, we need to ensure that train $v_{2}$ leaves station $i$. This is
done by

$\displaystyle\sum_{t\in T}\sum_{n=1}^{|N|}x_{i,n,t,v_{2}}\geq 1$
$\displaystyle\forall(v_{1},v_{2},i)\in Connections$ (8)

This already implies by equation 7 the arrival of train $v_{1}$ at the station
as elsewise equation (7) could become negative.

#### VIA-Nodes

$\displaystyle\sum_{t\in T}\sum_{j=1}^{|N|}x_{i,j,t,v}\geq 1$
$\displaystyle\forall(v,i)\in VIA-Nodes$ (9)

We can define a VIA-Node $i$ for a train $v$ as a tuple $(v,i)\in VIA-Nodes$.
This ensures that the vehicle $v$ passes node $i$ at some time point.

###  3.1 Robust Extension of the Model

Since extension measures have a long planning and construction horizon, the
long-term timetable, we use to determine the necessary extensions, is not
final. The demand or the traffic policy in the future could differ. To
incorporate this uncertainty of the timetable, we developed a robust version
of the model.

We model the uncertainty through two different extensions of the deterministic
model:

  * •

Optional trains and

  * •

Different timetables as scenarios.

First, we introduce optional trains, which can but do not have to be routed.
If they are not routed, this leads to a penalty which is implemented in the
objective function. The objective function is therefore changed to

$\min\sum_{ij\in E}k_{i,j}\cdot b_{i,j}-\sum_{v\in
V_{opt}}x_{d_{v},j,t,v}\cdot k_{v}$ (10)

where $k_{v}$ is the penalty of the optional train $v$ and $V_{opt}\subseteq
V$ is the set of optional trains. The departure and arrival are not necessary
for the optional trains. Therefore, the constraints that force the train to
leave and arrive are omitted for the optional trains $v\in V_{opt}$.

Arrival and Departure (11) $\displaystyle\sum_{t<t_{d_{v}}}x_{d_{v},j,t,v}=0$
$\displaystyle\forall v\in V$ $\displaystyle\sum_{t\geq
t_{d_{v}}}a_{d_{v_{t}},j_{t^{\prime}}}\cdot x_{d_{v},j,t,v}\geq 1$
$\displaystyle\forall v\in V\backslash V_{opt}$
$\displaystyle\sum_{t>t_{a_{v}}}x_{i,a_{v},t,v}=0$ $\displaystyle\forall v\in
V$ $\displaystyle\sum_{t\leq t_{a_{v}}}a_{i_{t^{\prime}},a_{v_{t}}}\cdot
x_{i,a_{v},t-tt(i,a_{v}),v}\geq 1$ $\displaystyle\forall v\in V\backslash
V_{opt}$

As all $x_{i,j,v,t}$ on the time-space-path of the train equal
$x_{d_{v},j,v,t}$ for some $t$ due to the flow constraints, the path is
completely true or all variables for this train are $0$. Therefore, the
capacity constraints and minimum headway time constraints do not change due to
the introduction of optional trains. Connection constraints and VIA-node
constraints are not implemented for optional trains.

Secondly, we integrated $m$ different scenarios into the optimization. Each
scenario has its own set of trains $V_{n}$ where each train in
$V=\bigcup_{n\in\\{1,\dots,m\\}}V_{n}$ has to be routed. Therefore the
departure and arrival time constraints (3), the flow constraints (3), the
connection constraints (7) and VIA-node constraints (9) stay the same and have
to be fulfilled for all trains. The constraints (13) and (12) describe
interactions between the trains and are used for every $V_{n}$ independently:

Capacity (12) $\displaystyle\sum_{v\in
V_{n},t\in(t_{0},t_{0}+\widetilde{t})}x_{i,j,t,v}\leq
c_{ij}+b_{ij}\cdot\widetilde{c_{ij}}$ $\displaystyle\forall ij\in E,t_{0}\in
T,t_{0}<\max\\{T\\}-\widetilde{t},n\in\\{1,\dots,m\\}$ Minimum Headway times
(13) $\displaystyle
x_{i,j,t_{1},v_{1}}\cdot(M_{ijv_{1}v_{2}}-(t_{2}-t_{1}))-(1-x_{i,j,t_{2},v_{2}})\cdot
max\\{0;M_{ijv_{1}v_{2}}-(t_{2}-t_{1})\\}\leq 0$ $\displaystyle\forall ij\in
E,t_{1}<t_{2}\in T,v_{1}\neq v_{2}\in V_{n},n\in\\{1,\dots,m\\}$

With these adaptions, we can model a uncertainty set of timetables in which
different scenarios with different set of trains can occur. The optimization
model can decide if the network can be extended to be able to schedule
optional trains. This happens through the comparison of the penalty for not
scheduling the train and the needed expansion cost for scheduling it. This
uncertainty set allows us to represent a uncertain long-term timetable for
operation in a few decades. Based on this, we can decide which infrastructure
measure can be taken to obtain a robust infrastructure.

###  3.2 Computational Results

The following test instances are run on an Intel i7-10700 CPU 2.90GHz machine
with 16 GB RAM using Gurobi 9.1.0 as optimizer. All the instances are solved
to optimality. The runtime is taken as the average out of 4 resolutions of the
same instance. The test instances are a small extract of one hour of the
German strategic timetable D̈eutschlandtaktïn the area of Dresden, see [23].

In the deterministic case, we obtain the running times in table 1:

Timesteps | Trains | Nodes | Arcs | Runtime [s] | Constraints | Variables
---|---|---|---|---|---|---
60 | 20 | 30 | 100 | 0,7 | 19875 | 1661
60 | 22 | 30 | 100 | 1,73 | 26621 | 2313
60 | 24 | 30 | 100 | 1,81 | 32265 | 2965
60 | 26 | 30 | 100 | 3,77 | 44879 | 4019
60 | 28 | 30 | 100 | 16,5 | 56507 | 5073
Table 1: Runtime of the optimization model for the deterministic case

Already in these small examples, we see the exponentially rising runtime and
number of constraints.

For the case with more scenarios, the runtime heavily depends on the partition
of the vehicles onto the scenarios as shown in table 2. If the vehicles are
approximately evenly distributed over the scenarios, the runtime decreases
compared to the deterministic case. This happens as the constraints for the
minimum headway time (which are the majority of the constraints) only occur
for vehicles that are in the same scenario. If a higher percentage of the
vehicles occur in both scenarios, the number of constraints and the runtime
rise.

Trains | in Sc 1 | in Sc 2 | in Sc 3 | in Sc 4 | Nodes | Arcs | Runtime[s] | Constr. | Variables
---|---|---|---|---|---|---|---|---|---
28 | 7 | 7 | 7 | 7 | 30 | 100 | 1,19 | 45502 | 5076
28 | 9 | 8 | 7 | 2 | 30 | 100 | 1,59 | 45284 | 5076
28 | 10 | 9 | 7 | 1 | 30 | 100 | 1,15 | 45612 | 5076
28 | 12 | 10 | 4 | 0 | 30 | 100 | 0,54 | 45612 | 5076
Table 2: Runtime of the optimization model for four scenarios and 60 timesteps

##  4 Complexity

In this section, we are going to analyze the complexity of the network design
problem for railways under timetable constraints. We will show, that the
problem is already NP-hard for very simple graphs and timetable structures.
Nevertheless, we will give some examples of specific graph classes in which
the problem is (pseudo-)polynomial solvable. We will focus on given, fixed
departure times. Otherwise, we obtain, even if our network consists of only
one edge, an interval scheduling problem. This is shown to be NP-hard in [2]
if there exists three or more departure options for any train.

The proof in this section is adapted from [12] who studied the budget-
constrained minimum cost flow problem. The network design problem for railways
we study in here is analogous to the dual problem of the budget-constrained
minimum cost flow problem. In this problem, a minimum cost flow is searched
for. The costs occur per unit of flow on each arc. Furthermore, a fixed-charge
is applied if the capacity of an arc is not sufficient. This fixed charge is
budgeted.

For this proof, we reduce the ExactCoverBy3Sets problem onto the railway
network design problem. In this proof, we omit VIA-nodes and connections for
simplicity.

######  Definition 4.1.

ExactCoverBy3Sets: Given a set $X$ with $3q$ elements and a set $C$ of
3-element subsets $C_{i}$ of $X$. Does there exist a subset
$\widetilde{C}\subset C$ where each element of $X$ occurs exactly once?

######  Theorem 4.2.

The railway network design problem is strongly NP-hard to solve even on
bipartite graphs and if all trains have the same, fixed departure time.

######  Proof 4.3.

We reduce the ExactCoverBy3Sets-problem to our problem. This problem is shown
to be strongly NP-hard in [6].

If we are given an instance of ExactCoverBy3Sets, we design a network
consisting of a source and a sink vertex, a vertex $v^{\prime}_{i}$ for each
subset $C_{i}$ and a vertex $v_{j}$ for each element $x_{j}\in X$. We then
construct arcs from the source $s$ to each node $v^{\prime}_{i}$ and from each
node $v_{i}$ to the sink $t$. Furthermore, we connect each node
$v^{\prime}_{i}$ with the nodes $v_{j}$ if $x_{j}\in C_{i}$. All of these arcs
have a capacity of 0. Let the expandable capacity
$\widetilde{c_{s,v^{\prime}_{i}}}$ of the arcs between the source and each
$v^{\prime}_{i}$ be $3$ and the expandable capacity of the arcs between each
$v_{i}$ and $t$ and the expandable capacity of the arcs between
$v^{\prime}_{i}$ and $v_{j}$ equal 1. Let all arcs have a travelling time of
1. The expansion cost for the arcs from the source to each $v^{\prime}_{i}$
are 3 while all other arcs have expansion costs of 0. Alternatively, this can
be modelled by defining no expandable capacity and having already a capacity
of 1 in these arcs.

For a small example given by $X=\\{x_{1},\dots,x_{6}\\}$ and
$C=\big{\\{}\\{x_{1},x_{3},x_{4}\\},\\{x_{1},x_{4},x_{5}\\},\\{x_{2},x_{5},x_{6}\\}\big{\\}}$,
the resulting network is shown in figure 1.

Figure 1: Resulting graph modelling an Exact3Cover as a railway network design
problem

In this example we can easily verify that the sets corresponding to
$v^{\prime}_{1}$ and $v_{3}^{\prime}$ would lead to an Exact3Cover.

We define a timetable with $X$ trains which all have the same earliest
departure time $0$ at their departure node $s$ and latest arrival time $3$ at
their arrival node $t$. This results in a given, specific departure time which
is the same for all trains (and as the minimal travel time is $3$, the
departure has to be at the earliest possible departure time). All trains start
at node $s$ and end at node $t$.

We consider the decision problem if it is possible to construct a routing with
a maximum cost of $3q$ such that $3q$ trains can be operated. We show that
this is true if and only if the given instance of ExactCoverBy3Sets is a yes-
instance.

Let there exist a schedule for $3q$ trains with cost $\leq 3q$. In this case
at most $q$ lines from the source to the nodes $v^{\prime}_{i}$ are expanded
and can therefore be used to route the trains. Since the capacity of the arcs
between $s$ and $v^{\prime}_{i}$ is 3 each and $3q$ trains are scheduled, we
obtain that at least $q$ arcs have to be expanded. Therefore, 3 trains have to
use each of the expanded arcs between the source and the $v^{\prime}_{i}$
since there exists exactly as many options to leave $s$ as there exists
trains. Each trains then has to use exactly one of the edges between
$v^{\prime}_{i}$ and $v_{j}$. Since only $q$ nodes $v^{\prime}_{i}$ are used,
all arcs between these $v^{\prime}_{i}$ and their corresponding three
successors $v_{j}$ have to be used to route the trains. This leads to a subset
$\widetilde{C}=\\{C_{i}|(s,v^{\prime}_{i})\text{ is expanded}\\}$ with $q$
sets such that all elements of $X$ are covered.

Let us now assume that there exists an ExactCoverBy3Sets given as
$\widetilde{C}\subset C$. With the same correspondence as before, we can
expand the arcs $(s,v^{\prime}_{i})$ for each $C_{i}\in\widetilde{C}$. This
leads to expansion costs of $3q$. Then three trains can be send through each
of these arcs. As there exists an exact cover, these trains can then proceed
to the three nodes $v_{j}$ following $v^{\prime}_{i}$ and each $v_{j}$ is only
visited by one train for which the capacity of 1 between $v_{j}$ and $t$ is
sufficient. This leads to a routing for all $3q$ trains with a cost of $3q$.
By this, we show that ExactCoverBy3Sets can be reduced polynomially onto the
network design problem for railways. If we have a train schedule, we can
easily verify if it is feasible. Therefore, the problem is in NP and we have
shown that it is strongly NP-hard.

As we used an extremely simple timetable and a svery simple costs structure to
prove the NP-hardness, a further simplification of the timetable does not seem
helpful to obtain easier solvable special cases. Therefore, we focus in the
following on specific graph classes to find (pseudo-)polynomial solvable
subproblems.

###  4.1 Complexity on Arborescences

In this subsection, we show that the network design problem for railway is
polynomially solvable on arborescences. An arborescences is an abstraction of
a tree onto digraphs. One vertex is defined as a root and all arcs point away
from the root. There exists no cycles. This leads to a unique path between
each pair of nodes.

######  Theorem 4.4.

The railway network design problem is polynomially solvable on an
arborescence, if all trains leave at a fixed departure time.

######  Proof 4.5.

If we are given a set of trains with an origin and a departure, the route for
each train is unique and therefore fixed. Hence, we can count the number of
vehicles on each arc (either per time step or directly if we set the root
vertex as departure for all trains). We can now verify if the routing is
feasible and if the (expandable) capacity on each arc is respected. It follows
directly which arcs need to be expanded. It can be easily seen that all these
steps can be executed in polynomial time in the input size.

However, as the trains potentially use the same lines, we again obtain an
interval scheduling problem if we omit the fixed departure time.

###  4.2 Complexity on Series-parallel graphs

We focus in this section on so called series-parallel graphs, for which we
obtain pseudo-polynomial solvability. We neglect here the fact that trains use
the capacity only during their respective travelling time on the line as we
would again obtain an NP-hard interval scheduling problem even on one arc.
(This corresponds to choosing the time interval for the capacity constraints
$\widetilde{t}=\max{t}$.)

######  Definition 4.6.

A series-parallel graph is a graph $G=(V,E)$ with a source $s$ and a sink $t$
which can be recursively constructed. Each series-parallel graph is build from
two series-parallel graphs through composition as shown in figure 2. The
simplest series-parallel graph is defined as one edge between $s$ and $t$.
This simple graph can then be expanded through series and parallel
composition. For a series composition of two series-parallel graphs $G$ and
$G^{\prime}$ the sink $t$ of $G$ is contracted with the source $s^{\prime}$ of
$G^{\prime}$ to obtain a new graph with source $s$ an sink $t^{\prime}$. For a
parallel composition the sources $s$ and $s^{\prime}$ are contracted to obtain
a new source $s$ as well as the sinks $t$ and $t^{\prime}$ are contracted to a
new sink $t$.

Figure 2: Series and parallel composition for a series-parallel graph

######  Theorem 4.7.

The network design problem for railway infrastructure under timetable
constraints is pseudo-polynomial solvable in

$O(|V|\cdot max_{v}(tt_{v})\cdot(1\cdot m+max_{v}(tt_{v})\cdot n+|V|\cdot m)$
(14)

for problems with a uniform given departure time.

######  Proof 4.8.

By $K_{G}(v,tt)$ we denote the minimal expansion cost on a series-parallel
graph for a number of $v$ vehicles and a maximum travel time for all vehicle
of $tt$ time steps.

We start with the case that the network $G=(N,A)$ only consists of one edge.
If the extendable capacity is smaller than the number of trains that want to
pass the edge, there exists no feasible solution and the minimum costs are
infinite. If the capacity without any extension is already greater or equal
than the number of trains, no extensions are necessary, so no costs occur.
Finally, if the extended capacity is greater or equal to the number of trains
while the existing capacity is smaller, the minimal costs are exactly the
extension costs for this edge.

Therefore, we obtain for the expansion costs for the edge $ij$

$K_{G}(v,tt)=\begin{cases}0&\text{if }v\leq c_{ij}\\\ k_{ij}&\text{if
}c_{ij}<v\leq c_{ij}+\widetilde{c}_{ij}\\\ \infty&\text{if }v\geq
c_{ij}+\widetilde{c}_{ij}\end{cases}$ (15)

#### Series Composition

For a series composition $G=G_{1}\circ G_{2}$ we show that

$K_{G}(v,tt)=\min_{tt_{1}\in(0,tt)}K_{G_{1}}(v,tt_{1})+K_{G_{2}}(v,tt-
tt_{1}).$ (16)

The minimal expansion costs for graph $G$ denoted by $K_{G}(v,tt)$ are clearly
smaller or equal to the sum of the costs $K_{G_{1}}(v,tt_{1})+K_{G_{2}}(v,tt-
tt_{1})$ for al $tt_{1}\in(0,tt)$ as this is a specific division of $G$ into
subgraphs for which the expansion scheme is minimal. This is not necessarily
the minimal expansion scheme for $G$.

Let us assume there exists a flow $x$ such that

$K_{G}(v,tt)<\min_{tt_{1}\in(0,tt)}K_{G_{1}}(v,tt_{1})+K_{G_{2}}(v,tt-
tt_{1}).$ (17)

Let $x_{1}$ denote the flow of $G$ restricted on $G_{1}$ and $x_{2}$ the flow
restricted on $G_{2}$. Then either $x_{1}$ has to have expansion costs lower
than $K_{G_{1}}(v,tt_{1})$ or $x_{2}$ has to have expansion costs lower than
$K_{G_{2}}(v,tt-tt_{1})$ for some $tt_{1}\in(0,tt)$. This contradicts the
definition of $K_{G_{1}}$ as minimal, as the flow $x$ has to pass $G_{1}$
completely and then pass $G_{2}$ completely by the definition of a series-
parallel graph.

#### Parallel Composition

Let $G=G_{1}|G_{2}$ be the parallel composition of two series-parallel graphs
$G_{1}$ and $G_{2}$. Then we obtain for the expansion cost of $G$

$K_{G}(v,tt)=\min_{v_{1}\leq v}K_{G_{1}}(v_{1},tt)+K_{G_{2}}(v-v_{1},tt)$ (18)

This holds, since we look for the optimal partition of the vehicles onto
$G_{1}$ and $G_{2}$ with respect to the expansion costs.

For each series-parallel graph, a tree can be constructed whose nodes refer to
the series and parallel compositions and whose leaves represent each single
edge. In [27] Valdes et al. show that such a tree for $G$ can be computed in
$O(m+n)$ where $m$ is the number of arcs and $n$ the number of vertices in
$G$. This tree has $O(m)$ inner nodes for parallel composition and $O(n)$
inner nodes for series composition as well as $m$ leaves which correspond to
single arcs. We compute $K_{G}(v,tt)$ for each node in the composition tree
for $tt\leq\max_{v}(tt_{v})$ and $v\leq|V|$. If the values are computed from
the leaves up, it can be assumed that the values for all subgraphs are known.
Therefore, the value $K_{G}$ can be computed in $O(1)$ for a single edge, in
$O(max_{v}(tt_{v}))$ for a series composition as we minimize over the
travelling times and in $O(V)$ for a parallel composition as we minimize over
the partition of the vehicles onto the parallel graphs. Then we obtain the
minimal expansion cost with $K_{G}(|V|,max_{v}(tt_{v}))$. Therefore, we obtain
a running time of

$O(|V|\cdot max_{v}(tt_{v})\cdot(1\cdot m+max_{v}(tt_{v})\cdot n+|V|\cdot m).$
(19)

##  5 Conclusion and Outlook

In this article, we analysed the network design problem under timetable
constraints. We have shown that the railway network design problem under
timetable constraints is NP-hard to solve on bipartite graphs and have shown
some easier solvable cases. Furthermore, we provided a formalization of the
macroscopic railway network design problem under timetable constraints for
both the deterministic and the robust problem formulation using time expanded
networks. The robust formulation can be used as a long-term timetable. This
still allows some flexibility to adapt the timetable to changing traffic
demand or political conditions during the long planning process of railway
infrastructure. Based on this timetable, the infrastructure planning can be
conducted more demand-oriented and economical.

Our next research goal is to find heuristics and approximations to reduce the
running time of the implementation of this NP-hard optimization problem.

#### Acknowledgements

This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) – 2236/1 and NI 1597/4-1.

## References

  * [1] Anthony Chen, Zhong Zhou, Piya Chootinan, Seungkyu Ryu, Chao Yang, and SC Wong.  Transport network design problem under uncertainty: a review and new developments.  Transport Reviews, 31(6):743–768, 2011.
  * [2] Julia Chuzhoy, Rafail Ostrovsky, and Yuval Rabani.  Approximation algorithms for the job interval selection problem and related scheduling problems.  Mathematics of Operations Research, 31(4):730–738, 2006.
  * [3] Alysson M Costa, Jean-François Cordeau, and Bernard Gendron.  Benders, metric and cutset inequalities for multicommodity capacitated network design.  Computational Optimization and Applications, 42(3):371–392, 2009\.
  * [4] Nadine Friesen, Tim Sander, Karl Nachtigall, and Nils Nießen.  Modelling time in the timetable-based railway network design problem.  Available at SSRN 4470191.
  * [5] Bosco García-Archilla, Antonio J Lozano, Juan A Mesa, and Federico Perea.  Grasp algorithms for the robust railway network design problem.  Journal of Heuristics, 19(2):399–422, 2013.
  * [6] Michael R Garey and David S Johnson.  Computers and intractability.  A Guide to the, 1979.
  * [7] Bernard Gendron, Teodor Gabriel Crainic, and Antonio Frangioni.  Multicommodity capacitated network design.  In Telecommunications network planning, pages 1–19. Springer, 1999\.
  * [8] Igor Grujičić, Günther Raidl, and Andreas Schöbel.  Variable neighborhood search for integrated timetable based design of railway infrastructure.  Electronic Notes in Discrete Mathematics, 47:141–148, 2015.
  * [9] Igor Grujičić, Günther Raidl, Andreas Schöbel, and Gerhard Besau.  A metaheuristic approach for integrated timetable based design of railway infrastructure.  Road and Rail Infrastructure III, 2017.
  * [10] Valérie Guihaire and Jin-Kao Hao.  Transit network design and scheduling: A global review.  Transportation Research Part A: Policy and Practice, 42(10):1251–1273, 2008.
  * [11] Andreas Heppe and Werner Weigand.  Spurplangestaltung und betriebliche infrastrukturplanung.  In Handbuch Eisenbahninfrastruktur, pages 467–523. Springer, 2019\.
  * [12] Michael Holzhauser, Sven O Krumke, and Clemens Thielen.  Budget-constrained minimum cost flows.  Journal of Combinatorial Optimization, 31(4):1720–1745, 2016.
  * [13] Michael Kuby, Zhongyi Xu, and Xiaodong Xie.  Railway network design with multiple project stages and time sequencing.  Journal of Geographical Systems, 3(1):25–47, 2001.
  * [14] Larry J Leblanc.  An algorithm for the discrete network design problem.  Transportation Science, 9(3):183–199, 1975.
  * [15] Richard M Lusby, Jesper Larsen, and Simon Bull.  A survey on robustness in railway planning.  European Journal of Operational Research, 266(1):1–15, 2018.
  * [16] Thomas L Magnanti and Richard T Wong.  Network design and transportation planning: Models and algorithms.  Transportation science, 18(1):1–55, 1984.
  * [17] Supakorn Mudchanatongsuk, Fernando Ordóñez, and Jie Liu.  Robust solutions for network design under transportation cost and demand uncertainty.  Journal of the Operational Research Society, 59(5):652–662, 2008\.
  * [18] Mir Saman Pishvaee, Jafar Razmi, and S Ali Torabi.  Robust possibilistic programming for socially responsible supply chain network design: A new approach.  Fuzzy sets and systems, 206:1–20, 2012.
  * [19] Gert-Jaap Polinder, Marie Schmidt, and Dennis Huisman.  Timetabling for strategic passenger railway planning.  Transportation Research Part B: Methodological, 146:111–135, 2021\.
  * [20] Tim Sander, Nadine Friesen, Karl Nachtigall, and Nils Nießen.  Robust railway network design based on strategic timetables.  submitted.
  * [21] Andreas Schöbel and Gerhard Besau.  Timetable based design of railway infrastructure.  In EURO-ZEL 20th International Symposium, 2012.
  * [22] Andreas Schöbel, GR Raidl, Igor Grujicic, Gerhard Besau, and Gottfried Schuster.  An optimization model for integrated timetable based design of railway infrastructure.  In Proceedings of the 5th International Seminar on Railway Operations Modelling and Analysis–RailCopenhagen, pages 765–774, 2013.
  * [23] SMA und Partner AG.  Zielfahrplan Deutschlandtakt: Dritter Gutachterentwurf Sachsen/Sachsen-Anhalt/Thüringen, 2020.
  * [24] Jacob Christian Spönemann.  Network design for railway infrastructure by means of linear programming.  PhD thesis, Aachen, Techn. Hochsch., Diss., 2013, 2013.
  * [25] Tommy Thomadsen and Thomas Stidsen.  The generalized fixed-charge network design problem.  Computers & Operations Research, 34(4):997–1007, 2007.
  * [26] Satish V Ukkusuri, Tom V Mathew, and S Travis Waller.  Robust transportation network design under demand uncertainty.  Computer-Aided Civil and Infrastructure Engineering, 22(1):6–18, 2007.
  * [27] Jacobo Valdes, Robert E Tarjan, and Eugene L Lawler.  The recognition of series parallel digraphs.  SIAM Journal on Computing, 11(2):298–313, 1982.
  * [28] Fang Zhao and Xiaogang Zeng.  Optimization of transit route network, vehicle headways and timetables for large-scale transit networks.  European Journal of Operational Research, 186(2):841–855, 2008\.


# Fractal Geometry of Angular Momentum Evolution in Near-Keplerian Systems

M. Atakan Gürkan
Astronomical Institute “Anton Pannekoek”, University of Amsterdam, Kruislaan
403, 1098 SJ Amsterdam, The Netherlands
Leiden University, Leiden Observatory, P.O. Box 9513, 2300 RA Leiden, The
Netherlands E-mail<EMAIL_ADDRESS>
(Accepted 201x XXX xx. Received 201x XXX xx; in original form 201x XXX xx)

###### Abstract

In this paper, we propose a method to study the nature of resonant relaxation
in near-Keplerian systems. Our technique is based on measuring the fractal
dimension of the angular momentum trails and we use it to analyze the outcome
of $N$-body simulations. With our method, we can reliably determine the
timescale for resonant relaxation, as well as the rate of change of angular
momentum in this regime. We find that growth of angular momentum is more rapid
than random walk, but slower than linear growth. We also determine the
presence of long term correlations, arising from the bounds on angular
momentum growth. We develop a toy model that reproduces all essential
properties of angular momentum evolution.

###### keywords:

methods: statistical — methods: N-body simulations — stellar dynamics —
Galaxy: centre

††pagerange: Fractal Geometry of Angular Momentum Evolution in Near-Keplerian
Systems–References††pubyear: 201x

##  1 Introduction

Dynamical systems where the gravitational potential is dominated by a single
mass are called near-Keplerian (Tremaine, 2005), since the orbits of smaller
masses in such systems are very close to Keplerian conic sections. A
characteristic feature of these systems is the existence of various distinct
timescales over which the dynamical variables change. For positions and
velocities, this is the dynamical time $t_{\mathrm{dyn}}\sim(a^{3}/GM)^{1/2}$,
where $M$ is the mass of the central object and $a$ is a given orbit’s
semimajor axis. The bound orbits (ellipses) precess over the precession time
$t_{\mathrm{prec}}\sim(M/Nm_{\star})t_{\mathrm{dyn}}$, where $m_{\star}$ is
the mass of smaller objects and $N$ is the number of these objects within $a$.
The timescale for the evolution of energy is the relaxation time,
$t_{\mathrm{rlx}}\sim(M^{2}/m_{\star}^{2}N\ln\Lambda)t_{\rm dyn}$. These
timescales form an hierarchy $t_{\mathrm{rlx}}\gg t_{\mathrm{prec}}\gg
t_{\mathrm{dyn}}$ and lead to three regimes for angular momentum evolution
(Tremaine, 1999).

For short timescales ($t\lesssim t_{\mathrm{dyn}}$), the changes in angular
momentum has no correlation since the torques felt at different parts of the
orbit vary rapidly. For intermediate timescales ($t_{\mathrm{dyn}}\lesssim
t\lesssim t_{\mathrm{prec}}$), where we are effectively averaging over orbits,
the changes in angular momentum are correlated, since the configuration of the
orbits change slowly. For long timescales ($t\gtrsim t_{\mathrm{prec}}$), the
correlation is lost again, since the orbits are randomized as they precess in
different directions at different rates. The presence of an intermediate
regime with enhanced angular momentum evolution was first recognized in the
literature by Rauch & Tremaine (1996). They called this process resonant
relaxation, since it is the result of a near resonance between angular and
radial frequencies of the orbit. It plays a central role in various scenarios
that are proposed to take place in the vicinity of supermassive black holes
(see, e.g., Magorrian & Tremaine, 1999; Alexander, 2008).

There are two major uncertainties regarding resonant relaxation. The first one
is the boundaries, especially the upper limit, for the intermediate regime.
The timescale over which the torques are coherent is evidently related to the
precession time $t_{\mathrm{prec}}$, but the detailed nature of this
relationship is unknown. Order of magnitude estimates for the timescales are
not sufficient, since in some cases the lifetimes of the stars in the systems
under question are comparable to the timescales of the dynamical processes
(for an example at the Galactic centre, see Madigan et al., 2009).

The nature of evolution of angular momentum over intermediate timescales is
also uncertain. In this regime, since the torques are correlated, the
evolution of angular momentum is not going to be diffusive ($\Delta
L\propto\Delta t^{1/2}$). Rauch & Tremaine (1996) seem to suggest that the
evolution in this regime is ballistic ($\Delta L\propto\Delta t$, see their
Eq.6 and Fig.1) and this is also adopted by Eilon et al. (2009); but their
simulations contain too few stars to draw conclusions in this regard. Indeed,
in this work we find that the nature of the angular momentum evolution in the
intermediate regime is between diffusive and ballistic.

We expose this behaviour and determine the timescales by carrying out
simplified simulations of near-Keplerian systems and a fractal analysis of the
evolution of angular momentum. In section 2, we give a brief review of
fractals and present a new method for determining fractal dimension, suitable
for trails obtained in numerical simulations. In section 3 we describe our
simulations and the fractal analysis of the angular momentum evolution. In
section 4, we demonstrate that such an evolution can be mimicked by a simple
random walk that retains a limited term memory. We discuss the limitations and
implications of our findings in section 5.

##  2 Fractal Dimension

Fractals (Mandelbrot, 1982) are geometric objects that exhibit a number of
properties that make them suitable for modelling physical processes and
natural structures. The defining property of fractals is that their Hausdorff-
Besicovitch dimension (hereafter fractal dimension, $D$) exceeds their
topological dimension. A number of methods for calculating the fractal
dimension is given by Mandelbrot (1982), here we give a brief sketch and
develop a new one that is suitable for our purposes.

A fractal is a self similar object; that is, part of it exhibits similar
properties to the whole, sometimes only in a statistical sense. If a self
similar object is made up of $n$ copies of itself, each of which is smaller
(in length) by $1/m$, then the object has fractal dimension $D=\ln n/\ln m$.
It is trivial to see that this leads to correct numbers for self similar
Euclidean objects. For example, any line segment can be thought to be composed
of $n=3$ identical copies of itself, each of which is smaller by $1/3$
($m=3$), giving $D=1$; or a rectangular prism can be cut into $n=8$ identical
copies of itself, each of which is smaller (in length) by $1/2$ ($m=2$),
giving $D=3$. We can use this technique to calculate the dimension of a well
known fractal, the Koch curve (Fig. 1). This curve consists of $n=4$ identical
copies of itself111Koch curve can also be seen as consisting of two copies of
itself, each of which is smaller by $1/\sqrt{3}$., each of which is smaller by
$1/3$, leading to a fractal dimension $D=\ln 4/\ln 3\sim 1.262$.

Figure 1: The Koch curve, obtained by repeatedly transforming a line segment.
The result of the first transformation is shown as a broken line in red, to
demonstrate the self similar structure more clearly. Each part of the Koch
curve resting on the straight parts of this broken line is a smaller but an
otherwise identical copy of itself.

Another way to obtain the fractal dimension is to measure the length of a
curve in increasing detail. Since as we use smaller and smaller rulers, we
will be resolving more and more detail; the measured length of the curve $L$,
is a function of the ruler length $\varepsilon$ (Mandelbrot, 1967):

$L(\varepsilon)\propto\varepsilon^{1-D}\,.$ (1)

As an example, let us assume that we measure the length of the Koch curve by a
given ruler and obtain $L_{0}$. When we reduce the ruler length by $1/3$, we
shall be traversing the smaller copies in exactly the same manner we traversed
the whole curve. Since there are four smaller copies, each of which will be
measured to have a length $1/3$ times the whole curve, we have
$L^{\prime}=4/3\times L$. To state in more general terms, if we modify our
ruler length by $\varepsilon^{\prime}\leftarrow(1/m)\varepsilon$, the measured
length changes by $L^{\prime}\leftarrow(n/m)L$. It is easy to see that this
scaling is satisfied when $L\propto\varepsilon^{1-(\ln n/\ln m)}$. This method
of determining the (fractal) dimension of a curve readily generalizes to real
life curves, which are self similar only in a statistical sense (Mandelbrot,
1967).

For our purposes, it is more convenient to interpret a curve as the motion
trail of an object. As we check the position of the object on this trail at
decreasing time intervals, we will be resolving more details and the total
trail length we calculate is going to increase. In other words, the trail
length is going to be a function of the sampling interval $\chi$. For example,
if we sample the motion of an object on the Koch curve $4$ times, we will be
measuring the trail shown as the broken line in Figure 1, leading to an
increase in length by $4/3$, with respect to going from the beginning to the
end in one step. In more general terms, when the sampling interval is modified
by $\chi^{\prime}\leftarrow(1/n)\chi$, the measured length becomes
$L^{\prime}\leftarrow(n/m)L$. In this approach, the fractal dimension is
determined through the relation

$L(\chi)\propto\chi^{\frac{1}{D}-1}\,.$ (2)

This technique of measuring fractal dimension of a trail can be easily applied
to the results from dynamical simulations.

##  3 Simulations

###  3.1 Simulation Method and Parameters

The system we simulated has three components. At the centre lies a stationary
supermassive black hole (SMBH) of mass $M_{\rm SMBH}=4\times 10^{6}M_{\sun}$.
Around the SMBH, we have $N_{\rm field}=1200$ field stars each with mass
$m_{\rm field}=2.5M_{\sun}$, semi-major axes distributed in a powerlaw cusp
$M_{\rm cusp}(r)\propto r^{3/2}$ from $a_{\rm min}=0.0001$ to $a_{\rm
max}=0.03\,{\rm pc}$, with eccentricities between $e_{\rm min}=0$ and $e_{\rm
max}=0.95$, distributed following the distribution function $g(e)\propto e$.
The final component is the massless test stars all with semimajor axis $a_{\rm
test}=0.01\,{\rm pc}$, and eccentricities $e_{\rm test}=0.2$, $0.75$, and
$0.85$ (12 stars each). Other orbital elements (inclination, longitude of the
ascending node, argument of pericentre and mean anomaly) of all stars are
picked at random. The exact values chosen do not matter too much, since over
the course of the simulation the eccentricities are randomized. We chose these
parameters to have a system that somewhat resembles the environment of S-stars
observed at our Galactic centre. There are large uncertainties regarding the
star distributions at this environment (Merritt, 2010, and references
therein), and $t_{\rm dyn}\ll t_{\rm prec}\ll t_{\rm rlx}$ hierarchy may not
even exist there. However, this condition would hold for a region around a
SMBH that developed a Bahcall-Wolf cusp (Bahcall & Wolf, 1976, 1977), so we
expect our method would be applicable to such systems.

The details of the code that is used for the simulations is going to be
explained in detail elsewhere, here we point out only the key features. Both
the field stars and the test stars feel the potential of the SMBH, including
the general relativistic (GR) correction that leads to the prograde precession
of the orbits222We use the treatment of Saha & Tremaine (1992, their eq. 30).
The test stars also feel the individual potential of the field stars, which
leads to retrograde precession of their orbits and changes in angular momentum
and energy. The field stars do not feel the individual potential of each
other, but instead see the potential of a smooth cusp that is consistent with
their distribution. This approximation decreases the time required for force
computation significantly and was already employed by Rauch & Tremaine (1996)
in their $N$-body simulations. For the interactions between the test stars and
field stars we use a softening kernel ($\mathrm{K}_{2}$) of Dehnen (2001),
with a softening length $10^{-4}$ pc. The units we adopted are
$G=M_{\mathrm{field}}=3000M_{\sun}=1\,\mathrm{pc}=1$.

The extended mass distribution in the background cluster precesses the orbits
in a retrograde fashion, and the GR effects lead to prograde precession. By
our choice of parameters, these two effects cancel each other for a star with
$a\sim 0.01\,\mathrm{pc}$ and $e\sim 0.6$. Mass precession becomes more
effective as semi-major axes get larger and eccentricities get smaller, while
the GR precession has the opposite behaviour.

We carry out the integration of the orbits using a high order Runge-Kutta-
Nyström method, discovered by Blanes & Moan (2000, their SRKN${}^{b}_{11}$).
We split the Hamiltonian into Keplerian and perturbation parts (Kinoshita et
al., 1991) to increase efficiency and avoid spurious precession. This scheme
advances the Keplerian orbital elements correctly, except for the truncation
and the roundoff errors (the largest accumulated error per star amounts to
$\sim 10^{-5}$ over the course of the whole simulation). In particular, the
evolution of the Runge-Lenz vector does not exhibit a linear drift as in the
case of potential energy-kinetic energy splitting (Hut et al., 1995). Our
treatment of the GR perturbation leads to a small error in mean motion, but
since we are interested in changes that take place over many orbits, this
error is not important. We use shared adaptive timesteps, but time-symmetrize
the integration with the method of Hut et al. (1995). Our simulations last 30
code units which corresponds to a few precession times of the slowest
precessing test stars.

Figure 2: The evolution of angular momentum (black) and energy (red) for a
test star with initial eccentricity $e_{0}=0.75$. For comparison, a Brownian
motion (blue) and a fractional Brownian motion (yellow; $c=18$, $n=4000$) data
as described in sec. 4 are also plotted. All curves are sampled at the same
abscissae, and the ordinates are adjusted to have unit variance.

We record the energy and angular momentum of each test star throughout the
simulation. In Figure 2 we show the evolution of these quantities for a star
with initial eccentricity $e_{0}=0.75$. Even by eye, it is possible to tell
that these quantities show a different behaviour. For comparison, in this
figure we also plot the curves generated by Brownian and fractional Brownian
motion, as described in Section 4.

###  3.2 Analysis of Simulation Results

Figure 3: Measured trail lengths as a function of duration of the sampling
intervals for angular momentum trails and the twice-broken power-law fits made
to them. Sets are artificially offset from each other to avoid confusion. The
gray lines in the background are to lead the eye and have slopes $-1/2$.

We measure the length of the angular momentum trail of each test star by
sampling it at intervals ranging from the full length of the simulation down
to a fraction of the orbital period. The dependence of the total measured
trail length on the sampling interval is shown in Figure 3 for a few stars
with different initial eccentricities. Starting from long sampling intervals
and moving towards shorter ones, a number of features can be observed on this
figure:

  * •

For long timescales, the curves do not have the slope $-1/2$ (corresponding to
dimension $D=2$, the value for Gaussian random walks (Mandelbrot, 1982)), but
are somewhat steeper. This is to be expected since the energy of an orbit
changes through relaxation and this process is much slower; hence the angular
momentum evolution is bounded unlike a true random walk, and cannot be
described by simple diffusion, even for long timescales.

  * •

There is a marked transition to a more coherent motion as indicated by a
decrease in the slope. The point of this transition can be determined with
reasonable accuracy for a given star, but it is not common for all stars.

  * •

Even though the slope decreases, it never becomes zero; hence the evolution of
angular momentum is never ballistic ($\Delta L\not\propto\Delta t$).

  * •

This slope can also be determined with reasonable accuracy for a given star,
but varies from star to star.

  * •

Looking at the eccentricity evolution of the stars reveals that the transition
point is later for the stars that moved into $e\sim 0.6$ region where
precession is slow. This is in harmony with the expectation that for more
slowly precessing stars, the coherent torques last longer.

  * •

The slope increases again for short timescales, but much before the period of
the stars is reached. This randomization of the torques is a result of the
stochastic nature of the processes that develop the torques and dominates the
shorter timescale randomization that would result from orbital motion, at
least down to the timescales we resolve.

Our results verify the presence of an intermediate regime, where the angular
momentum evolution is enhanced. The evolution of angular momentum in this
regime is not as rapid as ballistic growth $\Delta L\propto\Delta t$, but more
rapid than diffusive growth $\Delta L\propto\Delta t^{1/2}$. This manner of
evolution can be seen as the generalization of Brownian motion called
fractional Brownian motion. Mandelbrot & van Ness (1968) describes various
properties of this motion, along with applications. Mandelbrot (1982) gives
further generalizations and methods to produce such curves. Even though these
approaches are mathematically elegant and complete, in the next section we
propose a simpler model that is easier to attach a physical interpretation.

##  4 A Simple Toy Model for Angular Momentum Evolution

One dimensional Brownian motion for a particle’s position $P(t)$ can be
generated as follows. Let the motion over $\Delta t$ consist of $N$ steps with
equal duration $\delta t=\Delta t/N$. We choose an initial value $P_{0}$ and
at each step either increase or decrease the value of $P$ by $\delta P$. We
decide which action to take by generating a sequence of random numbers $X_{i}$
uniformly sampled from the interval $[-0.5,0.5]$ and choosing a threshold
value $q=0$. At each step, if $X_{i}<q$ we increase the value of $P$ and
otherwise decrease it. To make the variance of the motion independent of the
number of steps we choose $\delta P\propto 1/\sqrt{N}$. This scheme leads to
Brownian motion (Gaussian random walk) for small $\delta t$, i.e, a large
number of steps (Falconer, 2003, Chap. 16). It can be extended to a vector
variable by letting each component perform independent Brownian motions.

We can introduce correlations between the increments of the variable $P$ by
using a “repository”. For this, we keep track of the last $n$ values of
$X_{i}$ and let our threshold value be proportional to their average
$q=c\left<X_{i}\right>_{n}$, where $c$ is some constant. Here, $n$ determines
the length of the correlations and $c$ determines their strength. We generated
a few sets of data this way and measured the lengths of the resulting trails
with differing sampling durations.

Figure 4: Measured trail lengths as a function of duration of the sampling
intervals for various generated random walk data sets and the twice-broken
power-law fits made to them. The set in the middle (green) has $c=18$,
$n=4000$. The sets above it have varying $c=12,6,0$ and the sets below it have
varying $n=2000,1000,500$ (see the text for explanation of these parameters).
All sets have $N=200\,000$ points. To avoid confusion, sets are artificially
offset from each other; but to make comparison easier, top four and bottom
four sets are redrawn with dashed black lines at the top and bottom of the
figure with identical offsets. The gray lines in the background are to lead
the eye and have slopes $-1/2$.

The curves generated this way (Fig. 4) show very similar characteristics to
angular momentum evolution. They exhibit an intermediate regime with lowered
slope, whereas for long and short timescales the motion is more randomized.
The upper bound of the intermediate regime is determined by the parameter $n$:
the break occurs when the sampling interval matches $n\times\delta t$. The
other parameter $c$ determines the slope in the intermediate and short
timescale regime. Larger $c$ leads to a more coherent motion and a slope
closer to $0$.

Figure 5: Same as Fig. 4 but for bounded random walks. The slopes for long
sampling intervals are steeper than $-1/2$.

We also generated similar data with bounded random walks. For those, we
started with a (vector) variable of unit magnitude $|\mathbf{P}(t=0)|=1$, and
whenever a step led to $|\mathbf{P}|>1.5$, we took that step in the opposite
direction. This limit roughly corresponds to the limit experienced by an orbit
starting with initial eccentricity $e_{0}=0.75$,
$J_{\mathrm{circular}}/J_{0}\sim 1.5$. The results from this bounded random
walks are shown in Figure 5, exhibiting the long term correlations similar to
angular momentum evolution.

##  5 Discussion

In this work, we studied the evolution of angular momentum in near-Keplerian
systems by analyzing the outcome of $N$-body simulations. The simulation
method we use incorporates certain approximations. Our test stars are
massless, so they do not cause a back-reaction in the surrounding cluster.
Furthermore, the field stars all have the same mass. A mass spectrum would
change the granularity of the background potential, affecting the applied
torques and possibly the coherence timescale. Finally, since the field stars
do not see the granularity of their own potential, their angular momenta do
not evolve. This decreases the rate of randomization of the background
cluster, since vector resonant relaxation can change the orientation of the
orbits on timescales comparable to mass and GR precession for some
systems333We thank an anonymous referee for pointing out this shortcoming.
(for a comparison of these timescales for the Galactic centre, see Kocsis &
Tremaine (2010)). If the torques on a star mainly change because of the
rearrangement of the background cluster, rather than the reorientation of the
orbit, this decrease in randomization would alter the rate of angular momentum
evolution. All these approximations limit the domain of applicability of our
$N$-body simulation approach; however, none of them alter the essential
mechanism by which the angular momentum evolves.

We analyzed this evolution by calculating the fractal dimension of the angular
momentum trail. With this method it is possible to reliably determine the
onset of and the rate of evolution in different regimes. A key result of our
analysis is that the evolution of angular momentum is neither diffusive nor
ballistic in any regime, as was previously assumed in other studies (e.g.
Hopman & Alexander, 2006; Eilon et al., 2009). This seems to contradict the
results of the numerical experiments done by Rauch & Tremaine (1996). The
reason for this discrepancy is not clear, but we speculate that it arises from
the low number of stars used in that work.

We also developed a toy model that reproduces the features of angular momentum
evolution. This model has adjustable parameters with clear physical
interpretations. The relation between the appropriate values of these
parameters and the physical variables requires more detailed and extensive
analysis, which is currently underway. Apart from studying different initial
conditions, we also plan to analyze the components of the torque parallel and
perpendicular to the angular momentum separately. The nature of these torques
can be very different (Rauch & Tremaine, 1996; Gürkan & Hopman, 2007), so a
separate analysis should lead to a better understanding.

Koutsoyiannis (2002) developed fractional Gaussian noise generators (FGNGs)
similar to our repository model. In that work he compares autoregressive
moving average (ARMA) models to FGNGs, and finds that ARMA models are inferior
for describing long-term correlations. These models need to be modified to
have long term memory (Shumway & Stoffer, 2000), to describe angular momentum
evolution in near-Keplerian systems. Alternatively, they can be used to
describe short term correlations for torques, and long term correlations can
be introduced by taking physical bounds into account, as is done here. After
this paper was submitted, Madigan et al. (2010) also submitted a paper
containing their analysis of this problem with ARMA models. They use ARMA(1,1)
model, which fixes the value of the autocorrelation function for very large
timescales to zero (see their figure 4) and hence does not lead to any
correlations beyond a given time.

The computer programs used to generate and analyze the data are available from
the author upon request.

## Acknowledgments

This work is supported by a Netherlands Organization for Scientific Research
(NWO) Veni Fellowship. Most of the simulations were done on Lisa cluster,
maintained by SARA, the Dutch National High Performance Computing and
e-Science Support Center. I thank all SARA staff for doing an exceptional job
for maintaining this cluster and in particular to Walter Lioen for his help.
Usage of Lisa was possible through a grant (client number 10450) to Simon
Portegies Zwart. I am grateful to Clovis Hopman, Yuri Levin, Ann-Marie Madigan
and İnanç Adagideli for fruitful discussions on this topic, and an anonymous
referee for comments that improved this paper.

## References

  * Alexander (2008) Alexander T., 2008, in S. K. Chakrabarti & A. S. Majumdar ed., AIP Conference Series Vol. 1053, The Galactic Center as a laboratory for extreme mass ratio gravitational wave source dynamics. p. 79
  * Bahcall & Wolf (1976) Bahcall J. N., Wolf R. A., 1976, ApJ, 209, 214
  * Bahcall & Wolf (1977) Bahcall J. N., Wolf R. A., 1977, ApJ, 216, 883
  * Blanes & Moan (2000) Blanes S., Moan P. C., 2000, Journal of Computational and Applied Mathematics, 142, 313
  * Dehnen (2001) Dehnen W., 2001, MNRAS, 324, 273
  * Eilon et al. (2009) Eilon E., Kupi G., Alexander T., 2009, ApJ, 698, 641
  * Falconer (2003) Falconer K. J., 2003, Fractal Geometry: Mathematical Foundations and Applications. John Wiley & Sons, West Sussex
  * Gürkan & Hopman (2007) Gürkan M. A., Hopman C., 2007, MNRAS, 379, 1083
  * Hopman & Alexander (2006) Hopman C., Alexander T., 2006, ApJ, 645, 1152
  * Hut et al. (1995) Hut P., Makino J., McMillan S., 1995, ApJL, 443, L93
  * Kinoshita et al. (1991) Kinoshita H., Yoshida H., Nakai H., 1991, Celestial Mechanics and Dynamical Astronomy, 50, 59
  * Kocsis & Tremaine (2010) Kocsis B., Tremaine S., 2010, arXiv astro-ph, 1006.0001
  * Koutsoyiannis (2002) Koutsoyiannis D., 2002, Hydrological Sciences Journal, 47, 573
  * Madigan et al. (2010) Madigan A., Hopman C., Levin Y., 2010, arXiv astro-ph,1010.1535
  * Madigan et al. (2009) Madigan A., Levin Y., Hopman C., 2009, ApJL, 697, L44
  * Magorrian & Tremaine (1999) Magorrian J., Tremaine S., 1999, MNRAS, 309, 447
  * Mandelbrot (1967) Mandelbrot B. B., 1967, Science, 156, 636
  * Mandelbrot (1982) Mandelbrot B. B., 1982, The Fractal Geometry of Nature. Freeman, San Francisco
  * Mandelbrot & van Ness (1968) Mandelbrot B. B., van Ness J. W., 1968, SIAM Review, 10, 422
  * Merritt (2010) Merritt D., 2010, ApJ, 718, 739
  * Rauch & Tremaine (1996) Rauch K. P., Tremaine S., 1996, New Astronomy, 1, 149
  * Saha & Tremaine (1992) Saha P., Tremaine S., 1992, AJ, 104, 1633
  * Shumway & Stoffer (2000) Shumway R. H., Stoffer D. S., 2000, Time Series Analysis and Its Applications. Springer-Verlag, New York
  * Tremaine (1999) Tremaine S., 1999, in Impact of Modern Dynamics in Astronomy, IAU Colloquium 172, Resonant relaxation. p. 391
  * Tremaine (2005) Tremaine S., 2005, ApJ, 625, 143