# Influence of Cross-section Shape on Granular Column Collapses111

Teng Man  Herbert E. Huppert  Zaohui Zhang  Sergio A. Galindo-Torres
<EMAIL_ADDRESS>
###### Abstract

We investigate granular column collapses with different column cross-sections
and associate the cross-section shape influence with a finite-size analysis.
Previous research, confined to initially circular configurations, reviewed the
importance of granular column collapse studies and concluded that the run-out
distance scales with the initial aspect ratio. In this work, granular columns
with three different types of initial column cross-sections (square,
equilateral triangular, and rectangular cross-sections) are simulated using
the discrete element method (DEM). We explore how non-circular cross-sections
lead to different run-out distances. Based on the previously obtained finite-
size analysis, we further link the initial radius in different directions to
the relative size of a column and perform the finite-size analysis to explain
the cross-section influence. In the end, a universal relationship, which
includes frictional properties, relative sizes, and the cross-section
influence, is proposed. Our results are shown to have direct relevance to
various natural and engineering systems.

###### keywords:

Granular materials , Discrete element method , Column collapse , Different
cross-sections

††journal: Powder Technology

[inst1]organization=Key Laboratory of Coastal Environment and Resources of
Zhejiang Province (KLaCER), School of Engineering, Westlake
University,addressline=18 Shilongshan St., city=Hangzhou, postcode=310024,
state=Zhejiang, country=China

[inst2]organization=Institute of Theoretical Geophysics, King’s College,
University of Cambridge,addressline=King’s Parade, city=Cambridge,
postcode=CB2 1ST, state=, country=United Kingdom

##  1 Introduction

Granular materials are ubiquitous in natural and engineering systems. The
dynamics of granular column collapses play an important role and can bring
insights into understanding the kinetics and rheology of complex granular
systems in civil engineering, chemical engineering, pharmaceutical
engineering, food processing, and geophysical flows [1]. In recent decades,
progress has been made in terms of the rheological behavior of granular
systems [2, 3, 4]. Pouliquen et al. [2, 5] concluded that the rheology of
granular systems was controlled by inertial numbers ($\mu-I$ rheology), and
continuum models based on $\mu-I$ rheology were proposed accordingly [6, 7,
8], which have made significant contributions to understanding and modeling
the behavior of granular assemblies.

Initially, the behavior of granular column collapses was investigated to
better understand the post-failure behavior of discrete systems, such as
geophysical flows with similar flow mechanisms [9, 10]. Roche et al. [9, 11]
linked the dam-break experiments to the physics of pyroclastic flows, and
argued that the behavior of granular flows is similar to pyroclastic flows
when they are in certain regimes. Lube et al. [10, 12] and Lajeunesse et al.
[13] independently discovered that both the run-out distance and the final
deposition height can be determined by the initial aspect ratio of the column.
In particular, the normalized run-out distance,
$\mathcal{R}=(R_{\infty}-R_{i})/R_{i}$, where $R_{\infty}$ is the final radius
of the granular pile and $R_{i}$ is the initial radius of the granular column,
proportionally scales with the initial aspect ratio, $\alpha=H_{i}/R_{i}$,
where $H_{i}$ is the initial height of the column, when $\alpha<\alpha_{c}$
(where $\alpha_{c}$ is a transition point determined by experimental results),
and proportionally scales with $\alpha^{0.5}$ when $\alpha>\alpha_{c}$.

Zenit [14] performed discrete element method (DEM) simulations on two
dimensional (2D) granular column collapses, and confirmed that the shape of
the final deposition was mainly determined by the initial aspect ratio. Staron
and Hinch [15, 16] further investigated 2D granular collapses with DEM, and
found that the inter-particle frictional coefficient played an important role
in the run-out distance. Lacaze and Kerswell [17] studied axisymmetric
granular collapses to test the viscoplasticity of granular materials. Lagrée
et al. [7] implemented the $\mu-I$ rheology to a Navier-Stokes solver to study
the behavior of granular column collapses using a continuum approach. The
results of the continuum approach with $\mu-I$ rheology showed good agreement
with the results of DEM simulations. Farin et al. [18] even associated the
characteristics of granular column collapses with high-frequency seismic
signals, and used it to further evaluate geological events. With experiments,
Cabrerra et al. [19] and Warnett et al. [20] discovered that the relative size
of the granular column could also influence the normalized run-out distance.
They also found that, when the system size was large enough, the size effect
can be neglected.

Based on previous studies, we implemented dimensional analysis, and obtained
an effective aspect ratio,
$\alpha_{\textrm{eff}}=\sqrt{1/(\mu_{w}+\beta\mu_{p})}(H_{i}/R_{i})$, where
$\mu_{w}$ is the basal frictional coefficient, $\mu_{p}$ the inter-particle
frictional coefficient, and $\beta=2.0$ a constant that can be physically
interpreted as the ratio of contributions between inter-particle frictions and
particle/boundary frictions [21]. The effective aspect ratio includes the
influence of both friction and the initial aspect ratio. Furthermore, we
investigated the finite-size scaling of granular column collapses [22], and
come up with a universal scaling equation to describe the run-out behavior. We
showed that changing the system size could not only influence the normalized
run-out distance but also shift the characteristic aspect ratio $\alpha_{c}$
that marks the transition of granular column collapses from a quasi-static
regime to an inertial regime. A general scaling equation for the run-out
distance of columns with different sizes was obtained [22] and given by,

$\begin{split}\mathcal{R}=\left({R_{i}}/{d}\right)^{-\beta_{1}/\nu}\mathcal{F}_{r}\left[(\alpha_{\textrm{eff}}-\alpha_{c\infty})\left({R_{i}}/{d}\right)^{1/\nu}\right]\
,\end{split}$ (1)

where $\mathcal{F}_{r}[\cdot]$ is a scaling function, scaling parameters
$\nu=1.39\pm 0.14$ and $\beta_{1}=0.28\pm 0.04$ are obtained to best collapse
all the data, $\alpha_{c\infty}$ is the transitional effective aspect ratio
when the system size goes to infinity, and $d$ is the average particle
diameter. The current $\mathcal{F}_{r}[\cdot]$ function still lacks a
functional form, which limits the application of the finite-size scaling and
the generalization of our findings to other granular flow cases. However,
there is no contradiction between the finite-size solution and the power-law
scaling. In fact, they are both related by the renormalization group formalism
of statistical mechanics. The power-law relationship between $\mathcal{R}$ and
$\alpha_{\textrm{eff}}$ (of systems with the same relative size) has a unique
critical point, $\mathcal{R}_{c}$ and $\alpha_{c}$, which marks the transition
from quasi-static collapses to inertial collapses. The finite-size solution
with $R_{i}/d$ is utilized to, in some sense, collapse $\mathcal{R}$ and
$\alpha_{c}$ for systems with different relative sizes. In other words, the
power-law relationship is still embedded inside the $\mathcal{F}_{r}[\cdot]$
function, and should be recovered as the domain size diverges. We are also
aware that it is of vital importance to obtain the exact function for
$\mathcal{F}_{r}[\cdot]$, but the current results did not lead us to that
stage, and we will continue to work on this problem.

So far, only granular column collapses of circular cylinders or 2D granular
column collapses have been considered. We were somewhat surprised by our early
toy experiments of the collapse of rectangular granular columns, where the
final deposition approaches a circular pattern when the initial height is
sufficiently large. Thus, we wondered how granular columns with different
initial shapes collapse and what their final deposition patterns are. In this
paper, we investigate the influence of cross-section shape on granular column
collapses. Previous research focuses on pseudo-2D or axisymmetric granular
column collapses, where the shape of the cross-section does not play a role in
determining the run-out distance. However, non-axisymmetric cross-sections
lead to different initial column radii in different directions. For instance,
for a square cross-section, the initial radius in the diagonal direction is
$\sqrt{2}$ times of that in the direction pointing from the cross-section
center to the center of the edge. A different initial radius results in a
different initial aspect ratio, which further leads to a different normalized
run-out distance. This paper systematically explores the relationship between
normalized run-out distance and effective aspect ratio in different directions
resulting from different types of cross-sections, including square, rectangle
, and equilateral triangle (We choose these types of cross-sections with right
angles or sharp angles so that the direction of vertex and edge can be clearly
and conveniently defined). The paper is organized in the following way. In
Section 2, the simulation setup and the associated numerical method are
introduced. In Section 3, we present the simulation results, and further
discuss the physical insight that follows the simulation results. We further
implement the finite-size analysis that we obtained from our previous work
[22]. An experimental validation in Section 4 will also be provided to show
experimental evidence of the influence of cross-section shapes, before
concluding remarks are provided in Section 5.

##  2 Simulation method and setup

###  2.1 Discrete element method

To explore the behavior of the collapse of granular columns, we implement a
Voronoi-based sphero-polyhedral discrete element method (DEM) [23, 24], so
that we could obtain detailed particle-scale information during column
collapses. We note that the particle shapes could significantly affect the
deposition morphology but, in this study, we focus on using Voronoi-based
particles to investigate a general behavior similar to that of sand particles.

The sphero-polyhedra method was initially introduced by Pourning [25] for the
simulation of complex-shaped DEM particles. Later, it was modified by Alonso
Marroquin [26], who introduced a multi-contact approach in 2D allowing the
modelling of non-convex shapes and was extended to 3D by Galindo-Torres et.
al. [24]. A sphero-polyhedron is a polyhedron that has been eroded and then
dilated by a sphere. The result is a polyhedron of similar dimensions but with
rounded corners.

An advantage of the sphero-polyhedra technique is that it allows for an easy
and efficient definition of contact detection and force calculation between
particles. This is due to the smoothing of edges of all geometric features by
circles (in 2D) or spheres (in 3D). A particle is defined as a polyhedron,
i.e. a set of vertices, edges and faces, where each one of these geometrical
feature is dilated by a sphere.

Since two contacting particles are dilated by their sphero-radii $R_{1}$ and
$R_{2}$, there exists a contact when the distance between two geometric
features is less than the addition of their corresponding dilating radii, and
the corresponding contact overlap $\delta$ can be calculated accordingly. The
advantage of the sphero-polyhedra technique becomes evident since this
definition is similar to the contact force calculation of two spheres[27]. In
our simulation, three types of contacts (vertex-vertex contact, edge-edge
contact, and vertex-face contact) are considered. For these types of contacts,
we implement a Hookean contact model with energy dissipation to calculate the
interactions between particles as illustrated in Ref. [21]. At each time step,
the overlap and the tangential relative displacement between adjacent
particles, $\delta$, are checked, and the normal and tangential contact forces
can be calculated accordingly.

In this study, since we use Voronoi-based particles, no rolling resistance
need be considered. The motion of particles is then calculated by step-wise
resolution of Newton’s second law with the normal and contact forces mentioned
before. The same neighbor detection and force calculation algorithms have
already been discussed and validated in previous studies. This DEM formulation
has been validated before with experimental data [28, 27] and is included in
the MechSys open source multi-physics simulation library [29].

###  2.2 Simulation setup

Figure 1: (a) shows the collapse of a granular column with square cross-
section. Part of particles are cut from the figure to reveal the un-moving
core during the collapse. (b)-(e) show four types of cross-sections we use in
this work.

The DEM in this work has already been validated in previous studies to
simulate the behavior of granular materials in various conditions, such as
triaxial tests of granular soils [24], contact erosion phenomena [30], and
axisymmetric dry granular column collapses [21, 22]. We implement this DEM
model to granular column collapses with different initial cross-sections. As
shown in Fig. 1, we create granular columns with three different types of
initial cross-sections: 1) columns with square cross-sections, where the side
length $L_{s}=5$ cm; 2) columns with equilateral triangular cross-sections,
where the side length $L_{s}=6$ cm; 3) columns with rectangular cross-section
with side lengths 6 cm $\times$ 3 cm and 8 cm $\times$ 2 cm. The cross-section
size of granular columns are chosen so that the edge length is approximately
more than 10 times the particle size, and the area of the cross-section is
around 20 cm2. We exclude simulations with larger cross-section sizes because
of the correspondingly unacceptable computational time. For each type of
granular column, we vary the inter-particle frictional coefficient
($\mu_{p}=0.1,0.2,0.4$), while keeping the particle-boundary frictional
coefficient constant at $\mu_{w}=0.4$. For simulations with square and
triangular cross-sections, we also implement $\mu_{p}=0.6$. The height of
granular columns varies from 1 cm to 50 cm to obtain various initial aspect
ratios.

The original sphero-polyhedral granular packing is established using the 3D
Voronoi scheme with the Voro++ package [31]. The average particle size is
$d=0.2$ cm, so that the size of particles is similar to the size of typical
medium sized river sand and, at the same time, small enough compared to the
length of the cross-section of granular columns. Then 20% of particles are
randomly chosen and deleted to create an initial packing with a solid fraction
$\phi_{s}$ of 0.8 to keep the initial solid fraction the same as our previous
work [21]. The initial solid fraction of 80% is, on one hand, for reducing the
volume change during the collapse process, on the other hand, for introducing
randomness. We acknowledge the influence of the initial solid fraction on the
collapse and deposition of granular columns, and will conduct further studies
to address such influence in the future. Besides, we should be awared that,
due to the property of 3D Voronoi structure, the initial packing is more
similar to a fractured porous rock than to loosely packed sand.

Then, we remove the container in the simulation, release particles, and let
them flow under gravity. During the collapse of a granular column, the stored
potential energy will be transformed into kinetic energy and dissipated
through particle collisions. In the end, a stable granular pile can be
obtained with a final run-out distance, $R_{\infty}(\theta)$, where $\theta$
is the direction angle in the horizontal plane showing that, for a granular
column with non-circular cross-section, the final run-out distance varies with
respect to the direction in which we take the measurement. In detail, for
columns with square cross-sections, we present measurements in the eight
directions of 0∘, 45∘, 90∘, 135∘, 180∘, 225∘, 270∘, and 315∘. For columns with
equilateral triangular cross-sections, we present measurements in the six
directions of 30∘, 90∘, 150∘, 210∘, 270∘, and 330∘. For columns with
rectangular cross-section of length $6$ cm and width $3$ cm, we present
measurements in the eight directions of 0∘, 26.6∘, 90∘, 153.4∘, 180∘, 206.6∘,
270∘, and 333.4∘. Similarly, for columns with rectangular cross-section of
length $8$ cm and width $2$ cm, we present the measurement in the eight
directions of 0∘, 14∘, 90∘, 166∘, 180∘, 194∘, 270∘, and 346∘. These measuring
directions can be classified into two groups: i) edge directions (for
rectangular cross-section, the edge direction can be further classified into
longitudinal and width edge directions); and ii) vertex directions.

##  3 Results and discussions

###  3.1 Granular columns with square cross-sections

Figure 2: (a)-(d) show the polar distributions of particle numbers in
different directions for granular column collapses with square cross-sections
and faces parallel to the edge of the paper. The numbers on the vertical axis
represent the percentage of particles locating in certain directions, and the
numbers on the circumferential axis is the direction angle measured counter-
clockwise. (e) denotes the relationship between the normalized run-out
distance, $\mathcal{R}$, and the initial aspect ratio, $\alpha$, where markers
represent simulations $\mu_{p}=0.1$ in the edge direction , $\mu_{p}=0.2$ in
the edge direction , $\mu_{p}=0.4$ in the edge direction , $\mu_{p}=0.6$ in
the edge direction , $\mu_{p}=0.1$ in the vertex direction , $\mu_{p}=0.2$ in
the vertex direction , $\mu_{p}=0.4$ in the vertex direction , $\mu_{p}=0.6$
in the vertex direction . (f) plots $\mathcal{R}$ against the effective aspect
ratio, $\alpha_{\textrm{eff}}$, with the same markers as those in Fig. (e).

For granular columns with square cross-sections, we take measurements of
initial radius and final deposition radius in both edge and vertex directions.
The initial radius in the edge direction is $R_{i}^{e}=$ 2.5 cm, and the final
deposition radius in the edge direction is denoted as $R_{\infty}^{e}$.
Similarly, the initial radius in the vertex direction is
$R_{i}^{v}=\sqrt{2}R_{i}^{e}$, and the final deposition radius in the vertex
direction is denoted as $R_{\infty}^{v}$. Thus, we obtain the following key
parameters of initial aspect ratios and effective aspect ratios,

$\displaystyle\alpha^{e}$ $\displaystyle=H_{i}/R_{i}^{e},\ \
\alpha^{v}=H_{i}/R_{i}^{v},$ (2a) $\displaystyle\alpha_{\textrm{eff}}^{e}$
$\displaystyle=\alpha^{e}\sqrt{1/(\mu_{w}+\beta\mu_{p})},\
\alpha_{\textrm{eff}}^{v}=\alpha^{v}\sqrt{1/(\mu_{w}+\beta\mu_{p})},$ (2b)

where $\alpha^{e}$, $\alpha^{v}$ are the initial aspect ratios in the edge and
vertex directions, $\alpha_{\textrm{eff}}^{e}$ and $\alpha_{\textrm{eff}}^{v}$
are the effective aspect ratios in both edge and vertex directions, $\mu_{w}$
is the frictional coefficient between particles and the horizontal plane,
$\mu_{p}$ is the inter-particle frictional coefficient, and $\beta=2.0$ as
suggested by Ref. [21]. We could also obtain the normalized run-out distances
in both directions, $\mathcal{R}^{e}$ and $\mathcal{R}^{v}$, as

$\displaystyle\mathcal{R}^{e}$
$\displaystyle=(R_{\infty}^{e}-R_{i}^{e})/R_{i}^{e},$ (3a)
$\displaystyle\mathcal{R}^{v}$
$\displaystyle=(R_{\infty}^{v}-R_{i}^{v})/R_{i}^{v}.$ (3b)

The choice of only measuring run-out behaviors in both edge and vertex
directions for simplifying the analysis. In a more ideal condition, we can
measure the run-out distance in any arbitrary directions so that our analysis
can be more complete. In Fig. 2(a)-(d), we plot, for a range of initial
heights, the polar histogram of the number of particles in different
directions to show the initial and final pattern of the granular system. To
plot the polar histogram, we divide the round angle (360∘) into 100 pieces so
that d$\theta=3.6^{\circ}$. Each bar in the polar histogram is obtained by
counting the number of particles presenting within
$(\theta,\theta+\textrm{d}\theta]$, and dividing it by the total number of
particles in the system. Fig. 2(a) shows the initial state of the system,
where the majority of particles locate in the vertex direction, since the
initial radius of it is larger than that of the edge direction. For systems
with $H_{i}=$ 1 cm [Fig. 2(b)], the polar histogram still shows features of
squares, where more particles locate in the vertex direction. As we further
increase the initial height of the granular column, the final deposition
becomes approximately circular [Fig. 2(c) where $H_{i}=5$ cm] with an
increased horizontal length scale. This indicates that, as we increase the
initial aspect ratio of the granular column, the system evolves from a quasi-
static system[21], where most particles remain stationary during the collapse,
to an inertial system, where most particles participate in the collapse event,
and, during the evolution from a quasi-static system to an inertial collapsing
system, memory of its initial state is forgotten; hence, a circular deposition
pattern is developed. We note that, to be consistent with Ref. [21], the
quasi-static collapse is for describing the overall collapse and deposition
behaviors, and does not indicate that no particle inertia is playing a role
during the early stage of granular column collapses. Further increasing the
initial aspect ratio leads to more drastic changes to the deposition pattern,
where particles in the edge direction propagate much further than particles in
the vertex direction, as shown in Fig. 2(d) where $H_{i}=$ 25 cm. This usually
happens in the liquid-like collapse of granular columns in the phase diagram
in Ref. [21]. We obtain the general conclusion that, for a granular column
with square cross-section, the run-out distance in the edge direction is
larger than that in the vertex direction, especially when the initial aspect
ratio is large enough to trigger a liquid-like collapse. This is because the
initial aspect ratio in the edge direction is effectively larger than that in
the vertex direction for the same granular column.

We further explore the run-out distance in different directions in Fig. 2(e)
and (f), where we plot the relationship between $\mathcal{R}$ and $\alpha$,
and the relationship between $\mathcal{R}$ and $\alpha_{\textrm{eff}}$. In
Fig. 2(e), we show that the normalized run-out distance in both directions has
a similar behavior with increasing the initial aspect ratio. For granular
columns with the same initial aspect ratio but different inter-particle
frictional coefficients, a smaller inter-particle frictional coefficient
results in a larger run-out distance, which is the same for both edge and
vertex directions. Interestingly, the normalized run-out distance in an edge
direction is larger than that in a vertex direction, even when the initial
aspect ratio is the same.

We further investigate the relationship between $\mathcal{R}$ and
$\alpha_{\textrm{eff}}$, and expect that changing the $x-$axis from $\alpha$
to $\alpha_{\textrm{eff}}$ could lead to the collapse of normalized run-out
distance data of simulations with different frictional coefficients, according
to our previously published work [21]. Fig. 2(f) shows that, for simulation
results in either an edge or vertex direction, all the simulation data nicely
collapse onto one curve. However, the data of measurements in the vertex
direction are still smaller than those of measurements in the edge direction.
This is surprising, especially after we analyzed the finite-size scaling of
granular column collapses, where larger system sizes (larger $R_{i}/d$) lead
to larger normalized run-out distance when the effective aspect ratios,
$\alpha_{\textrm{eff}}$, are same [22]. The finite-size scaling of granular
column collapses implies that, with the same $\alpha_{\textrm{eff}}$, the
$\mathcal{R}$ in the vertex direction should be larger than that in the edge
direction. This paper aims to solve such a counter-intuitive situation (in
Sections 3.4 and 3.5).

###  3.2 Granular columns with equilateral triangular cross-sections

Figure 3: (a)-(d) show the polar distributions of particle numbers in
different directions for granular column collapses with equilateral triangular
cross-sections with a vertex at the top of the page. (e) and (f) denote the
relationship between $\mathcal{R}$ and $\alpha$, and the relationship between
$\mathcal{R}$ and $\alpha_{\textrm{eff}}$, where markers represent simulations
with $\mu_{p}=0.1$ in the edge direction , $\mu_{p}=0.2$ in the edge direction
, $\mu_{p}=0.4$ in the edge direction , $\mu_{p}=0.6$ in the edge direction ,
$\mu_{p}=0.1$ in the vertex direction , $\mu_{p}=0.2$ in the vertex direction
, $\mu_{p}=0.4$ in the vertex direction , $\mu_{p}=0.6$ in the vertex
direction .

We investigate the collapse of granular columns with equilateral triangular
cross-section. Again, to simplify the analysis, we only take measurement in
edge and vertex direction. Such a cross-section has three edge directions and
three vertex directions, so we can define $\alpha^{e}$, $\alpha^{v}$,
$\alpha_{\textrm{eff}}^{e}$, $\alpha_{\textrm{eff}}^{v}$, $\mathcal{R}^{e}$,
and $\mathcal{R}^{v}$ as we did for columns with a square cross-section. In
Fig. 3(a)-(d), we plot polar histograms of the initial state and final
depositions of columns with different initial heights. The polar distributions
of the initial state shows that initially most particles locate in the vertex
direction. Similar to columns with square cross-sections, when the initial
height is relatively small [Fig. 3(b)], the final pattern is still somewhat
triangular. Increasing the initial height leads to changing the deposition
pattern from approximately a triangle to approximately a circle [Fig. 3(c)],
which implies the memory of its initial triangular cross-section has been
largely forgotten during the collapse. Fig. 3(d) shows that, when the initial
height of a granular column is large enough, the final deposition pattern
becomes no longer circular. The run-out distance in the edge direction is much
larger than that in the vertex direction, which results in a reverse of the
equilateral triangle so that the original edge directions become vertex
directions in the final deposition. It looks as if the triangle was flipped
over. This surprising, but important, behavior is exactly the same as that
when the cross-section is a square.

We plot the relationship between $\mathcal{R}^{e}$ (or $\mathcal{R}^{v}$) and
$\alpha^{e}$ (or $\alpha^{v}$), as well as between $\mathcal{R}^{e}$ (or
$\mathcal{R}^{v}$) and $\alpha_{\textrm{eff}}^{e}$ (or
$\alpha_{\textrm{eff}}^{v}$) in Fig. 3(e) and (f). Similar to the case for a
axisymmetric granular column collapse, either $\mathcal{R}^{e}$ or
$\mathcal{R}^{v}$ approximately scales proportional to $\alpha^{e}$ or
$\alpha^{v}$ when the aspect ratio is less than a threshold, or scales
proportional to $(\alpha^{e})^{0.5}$ or $(\alpha^{v})^{0.5}$ when the aspect
ratio is larger than that threshold. Meanwhile, the normalized run-out
distance is influenced by the frictional coefficient, as we predicted in our
previous work [21]. Decreasing the inter-particle frictional coefficient can
not only change the $\mathcal{R}-\alpha$ curve but also change the critical
initial aspect ratio. Additionally, for a same initial aspect ratio,
$\mathcal{R}^{e}$ is apparently larger than $\mathcal{R}^{v}$. We shift the
$x-$axis to $\alpha_{\textrm{eff}}$ in Fig. 3(f), and find that the behavior
is similar to that of square granular columns. The gap between measurements in
the edge direction and in the vertex direction is even larger than that when
the cross-section is square shaped. This is explained by the fact that for a
square cross-section, $R_{i}^{v}/R_{i}^{e}$ is only 1.414, while the initial
radius ratio is 2.0 for a triangular cross-section.

###  3.3 Granular columns with rectangular cross-sections

Figure 4: (a)-(d) show the polar distributions of particle numbers in
different directions for granular column collapses with an initial $6\times 3$
rectangular cross-section. (e) and (f) denote the relationship between
$\mathcal{R}$ and $\alpha$, and the relationship between $\mathcal{R}$ and
$\alpha_{\textrm{eff}}$, where markers represent simulations with
$\mu_{p}=0.1$ in the short-edge direction , $\mu_{p}=0.2$ in the short-edge
direction , $\mu_{p}=0.4$ in the short-edge direction , $\mu_{p}=0.1$ in the
long-edge direction , $\mu_{p}=0.2$ in the long-edge direction , $\mu_{p}=0.4$
in the long-edge direction , $\mu_{p}=0.1$ in the vertex direction $\times$,
$\mu_{p}=0.2$ in the vertex direction $\times$, and $\mu_{p}=0.4$ in the
vertex direction $\times$.

Another set of simulations we present here is the collapse of granular columns
with two different rectangular cross-sections ($6\times 3$ and $8\times 2$),
where we take measurements in three directions: (1) the vertex direction,
where its initial radius is $R_{i}^{v}$; (2) the long-edge direction, where
its initial radius is $R_{i}^{el}$; and (3) the short-edge direction, where
its initial radius is $R_{i}^{es}$. In this study, the long-edge direction is
the direction parallel to the longer edge of the cross-section, and the short-
edge direction is the one parallel to the short edge of the cross-section. For
granular columns with a $6\times 3$ rectangular cross-sections, $R_{i}^{el}=3$
cm, $R_{i}^{es}=1.5$ cm, and $R_{i}^{v}\approx 3.35$ cm. Similarly, for the
granular columns with an $8\times 2$ rectangular cross-section, $R_{i}^{el}=4$
cm, $R_{i}^{es}=1$ cm, and $R_{i}^{v}\approx 4.12$. We notice that the
difference between $R_{i}^{el}$ and $R_{i}^{v}$ is small, so we expect that
the run-out distance in these two directions will not differ much from one
another. Also, since the initial aspect ratio in the short-edge direction is
always the largest among three directions, we expect to obtain a larger run-
out distances in this direction.

Figure 5: (a)-(d) show the polar distributions of particle numbers in
different directions for granular column collapses with an initial $8\times 2$
rectangular cross-sections. (e) and (f) denote the relationship between
$\mathcal{R}$ and $\alpha$, and the relationship between $\mathcal{R}$ and
$\alpha_{\textrm{eff}}$, where markers represent simulations with
$\mu_{p}=0.1$ in the short-edge direction , $\mu_{p}=0.2$ in the short-edge
direction , $\mu_{p}=0.4$ in the short-edge direction , $\mu_{p}=0.1$ in the
long-edge direction , $\mu_{p}=0.2$ in the long-edge direction , $\mu_{p}=0.4$
in the long-edge direction , $\mu_{p}=0.1$ in the vertex direction +,
$\mu_{p}=0.2$ in the vertex direction +, and $\mu_{p}=0.4$ in the vertex
direction +.

We show the initial polar distribution of particles in Fig. 4(a) and Fig.
5(a). In the initial configuration, most particles locate in the vertex and
long-edge directions The pattern varies when we change the initial height of
the column. When the initial height is short, the final deposition pattern
approximately remains rectangular [Fig. 4(b) and Fig. 5(b)]. Surprisingly,
when $H_{i}\approx 12$ for $6\times 3$ rectangular columns and $8\times 2$
rectangular columns, the final deposition pattern already becomes circular,
which indicates that the increase of the run-out distance in the short-edge
direction is much stronger than the other two directions [Fig. 4(c) and Fig.
5(c)]. With further increase of the initial height of the column, we observe
that the particles in the short-edge direction out-run those in the other two
directions since the initial aspect ratio in the short-edge direction is
larger than that in both the long-edge and the vertex direction, as can be
seen in Fig. 4(d) and Fig. 5(d).

Figure 6: Relationship between the front position in different directions and
time of granular columns with different initial heights and different types of
cross-sections.

What we have observed can also be validated by the plot of the
$\mathcal{R}$-$\alpha$ relationship and the plot of the
$\mathcal{R}$-$\alpha_{\textrm{eff}}$ relationship in Figs. 4(e)-(f) and Figs.
5(e)-(f). On one hand, since the initial radius in the vertex direction is
similar to that in the long-edge direction, their results of normalized run-
out distances are also similar. One the other hand, since the initial radius
of the short-edge direction is much smaller than that of the other two
directions, we see that the $\mathcal{R}^{es}$ \- $\alpha^{es}$ relationships
are always above the $\mathcal{R}^{el}$ \- $\alpha^{el}$ relationships and
$\mathcal{R}^{v}$ \- $\alpha^{v}$ relationships, and it is the same for the
$\mathcal{R}^{es}$ \- $\alpha_{\textrm{eff}}^{es}$ relationships.

Failing to collapse all the data with even the effective aspect ratios leads
us to re-think the way we normalize the run-out distances. Particles are
connected with each other with inter-particle collisions, which form force
networks. Mehta et al. [32] and our previous work [22] suggested that
particles tend to move collectively, especially when the initial height is
large enough to generate strong inertia among particles. Non-local effect tend
to be strong, so that using the initial radius in a direction to normalize the
run-out distance in that direction might not be an optimal choice.

We plot the relationship between the front position in different directions
and the collapse time in Fig. 6. The left column of the figure shows the
collapse of granular columns with a small initial height, while the right
column plots those with a large initial height. For short columns, a larger
initial radius often leads to a larger final front position. However, we can
observe that the front position in the edge direction (or short-edge direction
for rectangular cross-sections) is catching up with the front positions in
other directions. In Fig. 6(e), we find that particles in the long-edge
direction, while lagging behind those in the vertex direction at the
beginning, quickly overrun particles in other directions even with such a
small initial height. This shows clearly that particles in the edge direction
have a better chance to acquire a larger run-out distance. This phenomenon is
not obvious when the cross-section is an $8\times 2$ rectangle. That is
because, when the cross-section is an $8\times 2$ rectangle, the difference
between $R_{i}^{v}$ and $R_{i}^{el}$ is almost negligible. Figs. 6(b), (d),
(f), and (h) also confirm what we have seen in Figs. 2, 3, 4, and 5 that, at
the end of the collapse, the front position in the edge direction (or in the
short-edge direction) is always larger than that in other directions.

Figure 7: Evolution of polar histograms of tall granular column collapses with
respect to time. (a) Square cross-section, $H_{i}=25$ cm; (b) triangular
cross-section, $H_{i}=25$ cm; (c) $6\times 3$ rectangular cross-section,
$H_{i}=32$ cm; and (d) $8\times 2$ rectangular cross-section, $H_{i}=32$ cm.
Sub-figures in the same row share the same legend.

Figure 7 shows the evolution of polar histograms with respect to time. We only
plot the polar histogram of tall columns with four different cross-sections
[(a) Square cross-section, $H_{i}=25$ cm; (b) triangular cross-section,
$H_{i}=25$ cm; (c) $6\times 3$ rectangular cross-section, $H_{i}=32$ cm; and
(d) $8\times 2$ rectangular cross-section, $H_{i}=32$ cm]. Figs. 7(a)-(d) show
similar behavior that, when the column is tall enough, the deposition pattern
will quickly forget its original shape. The polar histogram quickly develop to
be approximately a circle within 0.3 s. However, for tall columns, as we have
shown in the previous section, the circular distribution could not last long
until it evolves to its final deposition histogram shown in Figs. 2(d), 3(d),
4(d), and 5(d). This shows that tall granular columns, which can generate
tremendous inertia, often collapse in a way that most particles localize into
certain directions, and those localized directions are often not the initial
dominant direction in original configurations.

We note that polar histograms may exaggerate the percentage of particles
presented in certain directions. Thus, we also plot contour plots of the
deposition, in Fig. 8, for tall granular column collapses with four different
cross-sections. We can see from Fig. 8(a) and (b) that particles initially
presented in the edge direction tend to travel farther, and more particles
seem to rest in the edge direction in the end. The accumulation of particles
in the short-edge direction is more obvious for systems with rectangular
cross-sections shown in Fig. 8(c) and (d), and particles in short edge
direction have larger deposition radius. The final deposition pattern is not
so exaggerate as the polar histograms, but it still reflects the behavior
shown in Figs. 2(e & f), 3(e & f), 4(e & f), and 5(e & f), and does not weaken
our argument that particles in certain directions prefer to travel farther
than other directions, and the classical method for quantifying the run-out
behavior is not enough to describe systems with non-circular cross-sections.
The detailed deposition pattern of systems with different initial cross-
sections is also intriguing and is worth of further explorations in future
studies.

Figure 8: Contour plots of the final deposition of granular columns collapses
with four different initial cross-sections: (a) a column with square cross-
section and $H_{i}=30$ cm; (b) a column with triangular cross-section and
$H_{i}=30$ cm; (c) a column with $6\times 3$ rectangular cross-section and
$H_{i}=32$ cm; (d) a column with $8\times 2$ rectangular cross-section and
$H_{i}=32$ cm.

###  3.4 Analysis with equivalent column radius

Figure 9: Relationship between the equivalent normalized run-out distance,
$\tilde{\mathcal{R}}$, and the effective aspect ratio,
$\alpha_{\textrm{eff}}$, of granular columns with (a) square cross-sections,
(b) equilateral triangular cross-sections, (c) $6\times 3$ rectangular cross-
sections, and (d) $8\times 2$ rectangular cross-sections. The markers are the
same as those in Figs. 2-5.

In the previous section, we argue that it is not appropriate to use the
initial column radius in certain directions to normalize the run-out distance
in that direction. In this section, we conveniently suggest an equivalent
column radius, the same for any direction of a column, to be defined to obtain
a so-called equivalent normalized run-out distance, $\tilde{\mathcal{R}}$,
defined by

$\begin{split}\tilde{\mathcal{R}}=(R_{\infty}-R_{i})/R_{\textrm{equiv}},\ \
R_{\textrm{equiv}}=\sqrt{A_{c}/\pi},\end{split}$ (4)

where $R_{\textrm{equiv}}$ is the equivalent column radius and $A_{c}$ is the
area of the cross-section. The advantage of using such a definition of
$R_{\textrm{equiv}}$ is that, on one hand, the area-equivalent radius can, in
a sense, help us compare granular columns with various cross-sections to
axisymmetric columns. On the other hand, when the cross-section is circular,
$R_{\textrm{equiv}}$ becomes the same as $R_{i}$, which makes it convenient to
formulate a universal equation to describe the behavior of granular column
collapses with any types of cross-section.

Figure 9 shows the relationship between $\tilde{\mathcal{R}}$ and
$\alpha_{\textrm{eff}}$ of simulations with different cross-sections. In Fig.
9(a), we plot the results of columns with square cross-sections. After
changing the $y-$axis to $\tilde{\mathcal{R}}$, all the simulation data
collapse nicely onto one curve. The relationship also performs well in terms
of granular columns with equilateral triangular cross-sections. However, one
thing that concerns us is that the
$\tilde{\mathcal{R}}^{v}-\alpha_{\textrm{eff}}^{v}$ relationship is slightly
above the $\tilde{\mathcal{R}}^{e}-\alpha_{\textrm{eff}}^{e}$ relationship,
which indicates that, in terms of the equivalent run-out distance, particles
in the vertex direction may travel longer equivalent distances than those in
the edge direction. Our concern is confirmed by the results of granular column
collapses with rectangular cross-sections [Fig. 9(c) and (d)], where, even
though there is almost no difference between the
$\tilde{\mathcal{R}}^{v}-\alpha_{\textrm{eff}}^{v}$ relationship and the
$\tilde{\mathcal{R}}^{el}-\alpha_{\textrm{eff}}^{el}$ relationship, the
$\tilde{\mathcal{R}}^{es}-\alpha_{\textrm{eff}}^{es}$ relationship is below
the other two relationships, especially for columns with $8\times 2$
rectangular cross-sections. This can be seen as a failure of implementing the
equivalent column radius, yet also an opportunity to perform the finite-size
scaling to our work. Similar to what we obtained in our previous work [22],
larger initial radii in vertex and long-edge directions lead to larger
relative system sizes, $R_{i}/d$, and larger normalized run-out distances (in
this case, it is the equivalent normalized run-out distance,
$\tilde{\mathcal{R}}$). We should also be careful about $R_{equiv}$ that the
currently definition of the equivalent cross-section radius, $R_{equiv}$, is
rather geometrical than physical. In future works, we should try to link this
geometric parameter to the dynamics of particle in different direction.

###  3.5 Finite-size analysis

Figure 10: Relationship between $\tilde{\mathcal{R}}(R_{i}/d)^{\beta_{1}/\nu}$
and $\mathcal{A}=(\alpha_{\textrm{eff}}-\alpha_{c\infty})(R_{i}/d)^{1/\nu}$,
which considers the size effect of granular columns with different cross-
sections. The markers are the same as those in Figs. 2-5 and the dashed curve
is a fitting line

We summarized the finite-size scaling of axisymmetric granular column
collapses as Eq. 1 in Section 1. In our previous work, we analyzed the size
effect of granular column collapses[22]. The relative system size, $R_{i}/d$,
ranged from 2.0 to 30.0. We found that, as we increase the relative system
size, both the transitional $\alpha$ and the transitional $\mathcal{R}$
decreases in a power-law manner, which led us to perform a finite-size scaling
with results shown in Eq. 1. We also confirmed the size effect by measuring
the relative correlation length scale, which also showed a power-law decay
with respect to $\alpha_{\textrm{eff}}$ and decreased as we increased the
system size. In this study, the normalized run-out distance is calculated
using the equivalent initial radius in each direction. We plot the
relationship between $\tilde{\mathcal{R}}(R_{i}/d)^{\beta_{1}/\nu}$ and
$\mathcal{A}=(\alpha_{\textrm{eff}}-\alpha_{c\infty})(R_{i}/d)^{1/\nu}$ of all
the simulation results in Fig. 10 with the same scaling parameters as in Ref.
[22] that $\nu\approx 1.39$ and $\beta_{1}\approx 0.28$. The following
relationship is expected

$\begin{split}\tilde{\mathcal{R}}=\left({R_{i}^{\theta}}/{d}\right)^{-\beta_{1}/\nu}\mathcal{F}_{r}\left[(\alpha_{\textrm{eff}}^{\theta}-\alpha_{c\infty})\left({R_{i}^{\theta}}/{d}\right)^{1/\nu}\right]\
,\end{split}$ (5)

where the superscript $\theta$ denotes the direction in which we take the
measurement, and can be replaced with $e$, $v$, $el$, or $es$. The results in
Fig. 10 agree with our expectation that all the simulation data of
$\tilde{\mathcal{R}}(R_{i}/d)^{\beta_{1}/\nu}-\mathcal{A}$ relationship with
different cross-sections and different measuring directions form a function,
$\mathcal{F}_{r}(\cdot)$.

The collapse of the $\tilde{\mathcal{R}}(R_{i}/d)^{\beta_{1}/\nu}-\mathcal{A}$
relationship indicates that the finite-size scaling is still valid when we
consider a non-circular cross-section. In other works related to the size
effect of granular column collapses, sizes effects refer to different sizes of
different columns. This work shows that, even within one granular column, size
effect should still be considered. We realize that, without a proper
definition of $\mathcal{F}_{r}(\cdot)$, it is impossible to accurately
calculate the run-out distance with any cross-section. Nevertheless, this work
helps formulate a universal solution to describe the run-out behavior
considering frictional coefficients, initial conditions, and geometric
factors, which is convenient to utilize when evaluating some geo-hazards, such
as landslides and pyroclastic flows, where both the system size and the
geometric factor are important for understanding the mobility of them. We note
that, the analysis, based on equivalent cross-section radii, may fail to work
for systems with extremely slender cross-sections. A more proper definition of
$R_{equiv}$ need to be discussed with a clearer physical definition in future
works.

##  4 Experimental Evidence of the Cross-section influence

To show the evidence of the influence of changing cross-section shapes, we
performed two sets of experiments of granular column collapses from both
square cross-section and rectangular cross-section tubes with sand particles
of diameters ranging from 0.1 mm to 0.2 mm. The aim of this section is showing
the orientation anisotropy presented in granular column collapses with non-
circular cross-section, instead of presenting a comparison between simulations
and experiments. We considered two different types of granular columns: (1)
granular columns with 50 mm $\times$ 50 mm square cross-sections; and (2)
granular columns with 40 mm $\times$ 21 mm rectangular cross-sections. For
both experimental setups, after placing sand particles into the plastic tube,
we measured the initial height of the granular packing, $H_{i}$. Sand
particles are dropped in from the top of the tube so that the initial
condition resembles a randomly loose packing of granular system. Then, the
tube was manually lifted to release all the particles to form a granular pile.
For columns with square cross-sections, we measured run-out distances in both
vertex and edge directions. For column with rectangular cross-sections, we
only measured run-out distances in both long edge and short edge directions as
previously defined in Fig. 1(d) and (e). The corresponding normalized run-out
distances, $\mathcal{R}$ and $\tilde{\mathcal{R}}$, can be obtained
accordingly. Since the frictional coefficient was kept as constant during
experiments, we use the initial aspect ratio, $\alpha$, as the $x-$axis in the
figures (Fig. 11), without considering the corresponding effective aspect
ratio, $\alpha_{\textrm{eff}}$.

Figure 11: (a) and (b) shows the initial condition and the final granular pile
of a granular column collapse from a tube with a square cross-section; (c) and
(d) show the experimental results of granular columns with 50 mm $\times$ 50
mm square cross-sections; and (e) and (f) show experimental results of
granular columns with 40 mm $\times$ 21 mm rectangular cross-sections.

Fig. 11(a) and (b) shows the initial condition and the final granular pile of
a granular column collapse from a tube with a square cross-section. We can see
that, even thought the final deposition pattern still shows some clue of the
initial square cross-section, Particles in the edge direction are catching up
to make the final pattern approximately circular. Figure 11(c) shows the
relationship between the normalized run-out distance and the initial aspect
ratio of granular columns with square cross-sections. For the same initial
aspect ratio, the normalized run-out distance in the edge directions is
apparently larger than that in the vertex directions which shows similar
behavior to that in Fig. 2(e). As we convert the $y-$axis into the equivalent
normalized run-out distance, $\tilde{\mathcal{R}}$ in Fig. 11(d), the data
collapse onto one curve, which is the same as that in Fig. 9(a). No obvious
size effect is shown in Fig. 11(d) because that $R_{i}^{v}$ is only
approximately 1.4 times of $R_{i}^{e}$.

Figure 11(e) shows experimental results of granular columns with 40 mm
$\times$ 21 mm rectangular cross-sections. Similar to Fig. 4(e), the
normalized run-out distance in the long-edge direction is larger than that in
the short-edge direction. As we shift the measuring direction from the edge
direction to the vertex direction, the transitional initial aspect ratio also
changes accordingly. In Fig. 11(f), we plot the relationship between
$\tilde{\mathcal{R}}$ and $\alpha$ in both short-edge and long-edge
directions. The equivalent normalized run-out distance in the long-edge
direction becomes larger than that in the short-edge direction, and
significant size effect can be observed since the initial radius in the long-
edge direction is almost twice the initial radius in the short-edge direction.
The experimental work shows significant similarities in both cross-section
shape influence and size effects with the simulation results we presented
earlier. This set of experiments confirms that, due to the influence of cross-
sesction anisotropy, particles tend to flow towards certain directions, and
size effect caused by different initial radii for granular columns with non-
circular cross-sections should be considered.

##  5 Conclusions

In this paper, we explore the influence of cross-section shapes on the
collapse of granular columns. Four different cross-sections are considered:
square, equilateral triangular, $6\times 3$ rectangular, and $8\times 2$
rectangular. We show that the normalized run-out distance in the edge
direction tends to be larger than that in the vertex directions when they have
similar initial aspect ratios or effective aspect ratios. It indicates that
the geometric factor plays an important role in determining the run-out
distance of granular column collapses. This might be important when dealing
with other problems in granular physics.

We propose that using an equivalent initial radius, instead of using the
initial radius in the measuring direction, to normalize the run-out distance
in that direction, performs better in terms of collapsing the simulation data.
In this study, we defined the equivalent initial radius as
$R_{\textrm{equiv}}=\sqrt{A_{c}/\pi}$, which is the radius of the circle with
the same area as the cross-section. Using the equivalent radius, we could
obtain the equivalent normalized run-out distance,
$\tilde{\mathcal{R}}=(R_{\infty}-R_{i})/R_{\textrm{equiv}}$ in any measuring
directions. However, the $\tilde{\mathcal{R}}-\alpha_{\textrm{eff}}$
relationship fails when considering the columns with rectangular cross-
sections, especially when the length of the cross-section is much larger than
the width. This leads us to perform the finite-size analysis to obtain the
$\tilde{\mathcal{R}}-\alpha_{\textrm{eff}}$ relationship, and results in a
universal relationship which could simultaneously consider the influence of
inter-particle friction, particle/boundary friction, initial aspect ratio, and
geometric factors.

We conclude that size effect should be considered not only when dealing with
different granular columns but also when calculating the run-out distance in
different direction within the same granular column. This becomes important
when we have a geo-hazard source with an irregular initial geometry. We also
note that, due to the nature of Voronoi diagrams, the initial packing of the
system is denser and more ordered than natural granular systems. We also
believe that the run-out behavior of granular columns should be linked to the
rheological properties of granular systems, which implies that the rheology of
granular systems with different system sizes and different frictional
coefficients should be included in the $\mathcal{F}_{r}[\cdot]$ function so
that the scaling law could become more physics-based. In fact, the body of
study we have introduced here offer clues on the true form of the rheology law
governing this behavior. Furthermore, a similar critical point should be
observed at some micro-mechanical quantity and should be included in the true
rheological law. Further investigations to link the behavior of idealized
granular system and realistic geophysical flows and include granular rheology
into the analysis of granular column collapses are still needed, and will be
presented in future publications.

## Acknowledgements

The authors acknowledge the financial support from the General Program (NO.
12172305) of the National Natural Science Foundation of China and Westlake
University, and thank the Westlake University Supercomputer Center for
computational resources and related assistance. The authors would like to
thank Prof. Ling Li for continuing helpful discussions related to this study.
All the simulations were conducted using MechSys (mechsys.nongnu.org)[33].

## Conflict of interest

The authors declare that they have no conflict of interest.

## References

  * [1] B. Andreotti, Y. Forterre, O. Pouliquen, Granular Media: Between Fluid and Solid, Cambridge University Press, 2013.  doi:10.1017/cbo9781139541008.
  * [2] G. D. R. MiDi, On dense granular flows, The European Physical Journal E 14 (4) (2004) 341–365.  doi:10.1140/epje/i2003-10153-0.
  * [3] M. Trulsson, B. Andreotti, P. Claudin, Transition from the viscous to inertial regime in dense suspensions, Physical Review Letter 109 (11) (2012) 118305.  doi:10.1103/PhysRevLett.109.118305.
  * [4] T. Man, K. Hill, Granular-slurry rheology and asphalt compaction, in: EPJ Web of Conferences, Vol. 249, EDP Sciences, 2021, p. 09010.  doi:10.1051/epjconf/202124909010.
  * [5] O. Pouliquen, C. Cassar, P. Jop, Y. Forterre, M. Nicolas, Flow of dense granular material: towards simple constitutive laws, Journal of Statistical Mechanics: Theory and Experiment 2006 (07) (2006) P07020.  doi:10.1088/1742-5468/2006/07/P07020.
  * [6] P. Jop, Y. Forterre, O. Pouliquen, A constitutive law for dense granular flows, Nature 441 (7094) (2006) 727.  doi:10.1038/nature04801.
  * [7] P.-Y. Lagrée, L. Staron, S. Popinet, The granular column collapse as a continuum: validity of a two-dimensional Navier–Stokes model with a $\mu$ (I)-rheology, Journal of Fluid Mechanics 686 (2011) 378–408.  doi:10.1017/jfm.2011.335.
  * [8] T. Barker, D. G. Schaeffer, P. Bohórquez, J. Gray, Well-posed and ill-posed behaviour of the $\mu({I})$-rheology for granular flow, Journal of Fluid Mechanics 779 (2015) 794–818.  doi:10.1017/jfm.2015.412.
  * [9] O. Roche, M. Gilbertson, J. C. Phillips, R. S. J. Sparks, Experiments on deaerating granular flows and implications for pyroclastic flow mobility, Geophysical Research Letters 29 (2002) 40–1–40–4.  doi:10.1029/2002GL014819.
  * [10] G. Lube, H. E. Huppert, R. S. J. Sparks, M. A. Hallworth, Axisymmetric collapses of granular columns, Journal of Fluid Mechanics 508 (2004) 175–199.  doi:10.1017/S0022112004009036.
  * [11] O. Roche, S. Montserrat, Y. Niño, A. Tamburrino, Experimental observations of water‐like behavior of initially fluidized, dam break granular flows and their relevance for the propagation of ash‐rich pyroclastic flows, Journal of Geophysical Research 113 (2008).  doi:10.1029/2008JB005664.
  * [12] G. Lube, H. E. Huppert, R. S. J. Sparks, A. Freundt, Collapses of two-dimensional granular columns, Physical Review E 72 (4) (2005) 041301.  doi:10.1103/PhysRevE.72.041301.
  * [13] E. Lajeunesse, J. Monnier, G. Homsy, Granular slumping on a horizontal surface, Physics of Fluids 17 (10) (2005) 103302.  doi:10.1063/1.2087687.
  * [14] R. Zenit, Computer simulations of the collapse of a granular column, Physics of Fluids 17 (3) (2005) 031703.  doi:10.1063/1.1862240.
  * [15] L. Staron, E. Hinch, Study of the collapse of granular columns using two-dimensional discrete-grain simulation, Journal of Fluid Mechanics 545 (2005) 1–27.  doi:10.1017/S0022112005006415.
  * [16] L. Staron, E. Hinch, The spreading of a granular mass: role of grain properties and initial conditions, Granular Matter 9 (3-4) (2007) 205.  doi:10.1007/s10035-006-0033-z.
  * [17] L. Lacaze, R. R. Kerswell, Axisymmetric granular collapse: a transient 3d flow test of viscoplasticity, Physics Review Letter 102 (10) (2009) 108305.  doi:10.1103/PhysRevLett.102.108305.
  * [18] M. Farin, A. Mangeney, J. D. Rosny, R. Toussaint, P. Trinh, Relations between the characteristics of granular column collapses and resultant high‐frequency seismic signals, Journal of Geophysical Research 124 (2019) 2987–3021.  doi:10.1029/2019JF005258.
  * [19] M. Cabrera, N. Estrada, Granular column collapse: Analysis of grain-size effects, Physical review. E 99 1-1 (2019) 012905.  doi:10.1103/PhysRevE.99.012905.
  * [20] J. Warnett, P. Denissenko, P. Thomas, E. Kiraci, M. Williams, Scalings of axisymmetric granular column collapse, Granular Matter 16 (1) (2014) 115–124.  doi:10.1007/s10035-013-0469-x.
  * [21] T. Man, H. E. Huppert, L. Li, S. A. Galindo-Torres, Deposition morphology of granular column collapses, Granular Matter 23 (3) (2021) 1–12.  doi:10.1007/s10035-021-01112-7.
  * [22] T. Man, H. E. Huppert, L. Li, S. A. Galindo-Torres, Finite-size analysis of the collapse of dry granular columns, Geophysical Research Letters 48 (24) (2021) e2021GL096054.
  * [23] P. A. Cundall, O. D. L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65.  doi:10.1680/geot.1979.29.1.47.
  * [24] S. A. Galindo-Torres, D. M. Pedroso, Molecular dynamics simulations of complex-shaped particles using voronoi-based spheropolyhedra, Physical Review E 81 (6) (2010) 061303–1–061303–9.  doi:10.1103/PhysRevE.81.061303.
  * [25] L. Pournin, T. M. Liebling, A generalization of distinct element method to tridimensional particles with complex shapes, Powders and Grains 2005 II (2005) 1375–1378.
URL http://infoscience.epfl.ch/record/77616

  * [26] F. Alonso-Marroquin, Spheropolygons: A new method to simulate conservative and dissipative interactions between 2d complex-shaped rigid bodies, EPL (Europhysics Letters) 83 (1) (2008) 14001.  doi:10.1209/0295-5075/83/14001.
  * [27] N. Belheine, J.-P. Plassiard, F.-V. Donzé, F. Darve, A. Seridi, Numerical simulation of drained triaxial test using 3d discrete element modeling, Computers and Geotechnics 36 (1-2) (2009) 320–331.  doi:10.1016/j.compgeo.2008.02.003.
  * [28] J. Cabrejos-Hurtado, S. Galindo Torres, D. Pedroso, Assessment of the mechanical behaviour of granular media by dem-based true triaxial tests, in: Advances of Computational Mechanics in Australia, Vol. 846 of Applied Mechanics and Materials, Trans Tech Publications Ltd, 2016, pp. 428–433.  doi:10.4028/www.scientific.net/AMM.846.428.
  * [29] S. Galindo-Torres, A coupled discrete element lattice boltzmann method for the simulation of fluid–solid interaction with particles of general shapes, Computer Methods in Applied Mechanics and Engineering 265 (2013) 107–119.  doi:10.1016/j.cma.2013.06.004.
  * [30] S. Galindo-Torres, A. Scheuermann, H. Mühlhaus, D. Williams, A micro-mechanical approach for the study of contact erosion, Acta Geotechnica 10 (2015) 357–368.  doi:10.1007/s11440-013-0282-z.
  * [31] C. Rycroft, Voro++: a three-dimensional voronoi cell library in c++., Chaos 19 4 (2009) 041111.  doi:10.1063/1.3215722.
  * [32] A. Mehta, Granular Physics, Cambridge University Press, 2007.  doi:10.1017/CBO9780511535314.
  * [33] MechSys: Multi-physics simulation library, http://mechsys.nongnu.org/, accessed: 2021-06-20 (2021).


A Python Code to Determine Orbital Parameters of Spectroscopic Binaries
Nicholas Milson, Caroline Barton, and Philip D. Bennett
Department of Physics and Atmospheric Science, Dalhousie University, Halifax,
NS B3H 4R2, Canada

###### Abstract

We present the open source Python code BinaryStarSolver that solves for the
orbital elements of a spectroscopic binary system. Given a time-series of
radial velocity measurements, six orbital parameters are determined: the long-
term mean, or systemic, radial velocity, the velocity amplitude, the argument
of periastron, the eccentricity, the epoch of periastron, and the orbital
period referred to by $\\{{\gamma,K,\omega,e,T_{0},P}\\}$ respectively. Also
returned to the user is the projected length of the semi-major axis,
$a_{1}\sin(i)$, and the mass function, $f(M)$. The determination of
spectroscopic orbits and masses is an example of another important area of
astrophysics, once the domain of professional astronomers, to which amateurs
can now make significant contributions. This code, available from GitHub, is
provided in support of that work, and should be of general use to the amateur
and professional astronomical community.

## 1\. Introduction

With the advent of commercially available research-grade spectrographs, the
amateur astronomy community, especially in Europe, is now acquiring
professional quality spectroscopic and radial velocity data for variable and
binary stars, see, e.g, Pollmann & Bennett (2020), and this trend is likely to
grow over the coming years. The determination of spectroscopic orbits and
masses is an example of another important area of astrophysics, once the
domain of professional astronomers, to which amateurs can now make significant
contributions. Here we present a summary of the problem of determining orbital
parameters from a time series of radial velocity observations of a
spectroscopic binary star, and provide the astronomical community with a new
open source Python code, BinaryStarSolver, to carry out this orbit solution.

The orbit of a star in a binary system is uniquely specified by seven orbital
parameters. These are the eccentricity $e$, the semi-major axis $a$, the time
(or epoch) of periastron $T_{0}$, the orbital period $P$, the argument of
periastron $\omega$, the position angle of the ascending node $\Omega$, and
the inclination angle $i$. For a spectroscopic binary system, radial
velocities of the stars can be inferred from analysis of their Doppler-shifted
spectral lines. The radial velocity, $V(t)$, of the primary star in a binary
system as a function of time $t$ is given by (Green 1985),

$V(t)=\gamma+\frac{na_{1}\sin
i}{\sqrt{1-e^{2}}}\,\,\big{\\{}\cos[v(t)+\omega]+e\cos\omega\big{\\}}$ (1)

where $n=2\pi/P$ is the mean motion, $a_{1}$ is the semi-major axis of the
primary star’s orbit, $\gamma$ is the long term mean or systemic velocity of
the binary, and $v(t)$ is the true anomaly as a function of time. It can be
seen in this expression that $a_{1}$ and $i$ are coupled. A separate
determination of each parameter is not possible without additional
information, independent of radial velocities. Therefore, it is standard
practice to define the new parameter, the radial velocity semi-amplitude of
the primary star $K_{1}$, as

$K_{1}=\frac{na_{1}\sin i}{\sqrt{1-e^{2}}}.$ (2)

The velocity semi-amplitude of the secondary star in the binary, $K_{2}$, is
defined in an analogous manner using $a_{2}$. The expression for radial
velocity as a function of time for the $j$-th star ($j=1,2$) in the binary
system then becomes

$V_{j}(t)=\gamma+K_{j}\big{\\{}\cos[v(t)+\omega]+e\cos\omega\big{\\}},\hskip
7.22743ptj=1,2$ (3)

If the spectra of both stars can be observed and separated, then we have the
case of an SB2 (double-lined) spectoscopic binary. In this case, knowledge of
$K_{1}$ and $K_{2}$ of both stars suffices to determine the ratio of stellar
masses: $M_{1}/M_{2}=K_{2}/K_{1}=a_{2}/a_{1}$. Without additional information
about the inclination angle $i$, individual stellar masses can not be found;
the best that can be accomplished is the determination of the mass function
$f(M)=(M_{1}+M_{2})\sin^{3}i$.

An additional constraint is provided by Kepler’s third law (Green 1985).

$M_{1}+M_{2}=\frac{4\pi^{2}a^{3}}{GP^{2}}=3.985\times
10^{-20}\,\frac{(a_{1}+a_{2})^{3}}{P^{2}}$ (4)

where $a=a_{1}+a_{2}$ is the semi-major axis of the relative orbit of the two
stars, and the constants have been evaluated to yield masses in solar units
($M_{\odot}$) given $a_{j}$ in km and $P$ in days. Then, the mass function can
be written

$f(M)=(M_{1}+M_{2})\sin^{3}i=3.985\times 10^{-20}\,\frac{(a_{1}\sin
i+a_{2}\sin i)^{3}}{P^{2}}$ (5)

where

$a_{j}\sin i=\frac{K_{j}\sqrt{1-e^{2}}}{n},\hskip 7.22743ptj=1,2$ (6)

For the special case of eclipsing binaries, the inclination $i$ must be close
to $90^{\circ}$ for eclipses to occur, and in that case we can typically
assume $\sin i=1$. Then, the total mass $M_{1}+M_{2}$, and the mass ratio
$M_{1}/M_{2}$ are both known, and the individual stellar masses $M_{1}$ and
$M_{2}$ are determined, as are $a_{1}$ and $a_{2}$. Therefore, a complete
determination of the orbit and masses of the stars in an eclipsing, double-
lined spectroscopic binary (SB2) is possible.

For convenience, we will from now on normally use $K$ without a subscript to
refer to either star, and if only one star is implied, the context will make
this clear. Notice that time only appears in this expression implicitly. Since
the measurements made are a function of time, not true anomaly, a relationship
between the two variables is needed. This follows from two equations: first,
an expression between the true anomaly and the eccentric anomaly, $E$, and
second, between the eccentric anomaly and time (the latter is known as
Kepler’s Equation, Green 1985).

$\tan\frac{v}{2}=\frac{1+e}{1-e}\tan\frac{E}{2}$ (7) $t=\frac{E-e\sin
E}{n}\,+\,T_{0}$ (8)

Note that $\Omega$ does not appear in any of the previous expressions. This is
to be expected since $\Omega$ specifies the orientation of the orbit in the
plane of the sky, and so has no effect of the observed radial velocity. This
also means that it cannot be determined solely from radial velocity
measurements. Therefore, we do not concern ourselves further with $\Omega$
here. This leaves six parameters to be determined:
$\\{\gamma,K,\omega,e,T_{0},P\\}$.

In the case of an eclipsing binary, the period $P$ may be found more
accurately through the observation of eclipses. If this is the case, in the
provided code, the user can opt to provide the known period and solve for only
$\\{\gamma,K,\omega,e,T_{0}\\}$. Generally though, $P$ is not predetermined
and the following discussion assumes this case.

## 2\. Discussion of procedure

### 2.1. Minimization

To determine the orbital parameters, the given radial velocity data is fit
with a curve of the form of (3), using a nonlinear least-squares approach to
minimize the residual:

$\chi^{2}=\sum_{i=1}^{N}w_{i}\bigg{(}V_{i}-\Big{\\{}\gamma+K\big{[}\cos(v(t_{i})+\omega)+e\cos\omega\big{]}\Big{\\}}\bigg{)}^{2},$
(9)

where $V_{i}$ is the observed radial velocity at time $t_{i}$, and $w_{i}$ is
the weight of the measurement.

To carry out this minimization, the Levenberg-Marquardt algorithm (also known
as the damped least-squares method) is used. This algorithm combines Newton’s
method and the gradient descent method through the introduction of a dampening
parameter $\lambda$ to give both a more efficient and reliable convergence
than either method on their own (Gavin 2011, Marquardt 1963).

Specifically, in this method, the entries along the diagonal of the Hessian
matrix, $\boldsymbol{H}$, are multiplied by $(1+\lambda)$. Refer to Appendix A
for the full Hessian matrix. When the parameter $\lambda$ is large
($\lambda\gg 1$), then the Hessian matrix
${\boldsymbol{H}}\approx\lambda{\boldsymbol{I}}$, where $\boldsymbol{I}$ is
the $6\times 6$ diagonal identity matrix, and then the change in parameter
vector $\Delta\boldsymbol{p}$ follows the gradient descent direction. As
$\lambda\rightarrow 0$, the change $\Delta\boldsymbol{p}$ approaches that of
Newton’s method. Since a change in the direction of the negative of the
gradient ensures a reduction in the local value of $\chi^{2}$, an initial
estimate should be made using a suitably large value of $\lambda$.

Apart from the introduction of the $\lambda$ coefficient, however, the
minimization procedure closely follows Newton’s method. So, where the current
estimate of parameters is the vector, $\boldsymbol{p}_{n}$,
$\boldsymbol{p}_{n+1}$ may be found by solving for
$\Delta\boldsymbol{p}=\boldsymbol{p}_{n+1}-\boldsymbol{p}_{n}$ from the system
of equations

$\big{[}\boldsymbol{H}\chi^{2}(\boldsymbol{p}_{n})\big{]}\Delta\boldsymbol{p}=-\nabla\chi^{2}(\boldsymbol{p}_{n}),$
(10)

where $\nabla\chi^{2}(\boldsymbol{p}_{n})$ is the gradient of $\chi^{2}$ at
$\boldsymbol{p}_{n}$ (Press et al. 1992). Complete details of the Hessian
matrix and the partial derivatives used in the minimization are given in
Appendix A.

Upon running the minimization procedure, an initial of $\lambda=3$ is used.
This initial value of $\lambda$ was found to be a reasonable compromise
between ensuring convergence while retaining accuracy of the final parameter
values. If after one iteration the residual decreases, $\lambda$ is reduced by
a factor of 9. If instead the residual increases, $\lambda$ is increased by a
factor of 11. In this case, the current parameter vector
$\boldsymbol{p}_{n+1}$ is rejected and the next iteration starts using the
previous parameter vector $\boldsymbol{p}_{n}$ again (Gavin 2011, Press et al.
1992, Transtrum & Sethna 2012). The minimization procedure is then left to run
until the change in the residual becomes negligible (either an absolute change
$<0.01$, or a fractional change $<10^{-3}$ ), at which point sufficient
convergence is achieved (Press et al. 1992).

The parameters returned are then used to calculate $a_{1}\sin(i)$ and the mass
function (see Appendix A) and the stellar parameters
$\\{\gamma,K,\omega,e,T_{0},P,a_{1}\sin i,f(M)\\}$ are returned to the user
with their respective uncertainties, as discussed in Section 2.4.

### 2.2. Companion star

Included in this code, as well, is a separate function that solves for the
orbital parameters of a star given the parameters of its companion (and radial
velocity data). This allows for a more accurate determination of the orbital
parameters of the star if less data is available for it compared to its
companion.

In a binary star system, the orbital elements of the companion star are nearly
all identical to that of the primary star. The eccentricity, inclination
angle, orbital period, and epoch of periastron passage are all the same for
both stars and the argument of periastron of the companion star only differs
from the argument of periastron for the primary star by 180°(Green 1985). This
leaves the semi-major axis, $a_{2}$, (and by extension the velocity amplitude
$K_{2}$) as the only undetermined parameter. Since $K_{2}$ is the only
variable parameter, it is determined using Newton’s method.

As arguments, the companion function takes the determined orbital elements
$\\{{\gamma,K_{1},\omega_{1},e,T_{0},P}\\}$ of the primary star, and returns
the orbital elements $\\{\gamma,K_{2},\omega_{2},e,T_{0},P,a_{2}\sin
i,f(M)\\}$ of the companion star, requiring only a few radial velocity points
and their respective times to accurately determine $K_{2}$.

### 2.3. Initial estimates

Like other minimization methods, the Levenberg-Marquardt algorithm will often
fail to find the global minimum if the initial estimates are not accurate
enough. Thus, unless initial estimates are provided by the user, the
minimization routine begins by carefully approximating the parameters
$\\{\gamma,K,\omega,e,P,T_{0}\\}$ from the given data.

First, $\gamma$ can be estimated using the average of the radial velocity
data, $V_{avg}=\sum V_{i}/N$, where $N$ is the number of data points. To
approximate the other parameters, equations (3), (7), and (8) must be
manipulated. Beginning with equation (3), it can be seen that the maximum
radial velocity occurs when $\cos(v+w)=1$ and the minimum radial velocity
occurs when $\cos(v+w)=-1$ (Green 1985). At $V_{max}$ and $V_{min}$, (3) then
reduces to two equations

$V_{max}=\gamma+K\big{[}1+e\cos w\big{]}$ (11)
$V_{min}=\gamma+K\big{[}-1+e\cos w\big{]}$ (12)

By adding/subtracting equations (11) and (12), it can be shown that

$\frac{V_{max}-V_{min}}{2}=K$ (13) $\frac{V_{max}+V_{min}}{2}=\gamma+Ke\cos w$
(14)

From equation (13), $K$ can be estimated using $V_{max}$ and $V_{min}$ from
the data.

$T_{0}$ and $\omega$ can be related to each other through two equations. First
by rearranging equation (7),

$E_{max}=2\arctan\bigg{[}{\Big{(}\frac{1+e}{1-e}\Big{)}}^{-\frac{1}{2}}\tan\Big{(}\frac{v_{max}}{2}\Big{)}\bigg{]}$
(15)

Where $E_{max}$ and $v_{max}$ are the eccentric and true anomaly when
$V=V_{max}$ (which may be found from equations (7) and (8)). However,
$\cos(v+\omega)=1$ at $V_{max}$, $v_{max}=-\omega$ and so,

$E_{max}=2\arctan\bigg{[}{\Big{(}\frac{1+e}{1-e}\Big{)}}^{-\frac{1}{2}}\tan\Big{(}\frac{-\omega}{2}\Big{)}\bigg{]}$
(16)

Next, rearranging equation (8), we find

$T_{0}=t_{max}-\frac{1}{n}\big{[}E_{max}-e\sin E_{max}\big{]}$ (17)

Here $t_{max}$ is the time corresponding to $V_{max}$.

Now, armed with these equations, initial estimates may be found using the
method outlined in the following list.

  1. 1.

Estimate $\gamma$ by averaging the radial velocity data

  2. 2.

Estimate $K$ using (13)

  3. 3.

Unless $P$ is supplied by the user, find $P$ using two methods (further
described in 2.3.1)

    1. (a)

From the two points lying close to the same $V$ axis, separated by an integer
number of periods

    2. (b)

By averaging the separation between maxima, minima, and points crossing the
same $V$ axis

  4. 4.

Iterate $e$ from 0.01 to 0.99 with a step size of 0.01 and for each value

    1. (a)

Find $\omega$ using (14)

    2. (b)

Find $T_{0}$ using (15) and (16)

    3. (c)

Using the estimates for $\\{\gamma,K,\omega,e,T_{0}\\}$ and

      1. i.

$P$ from 3(a), calculate the sum of the squared deviations

      2. ii.

$P$ from 3(b), calculate the sum of the squared deviations

    4. (d)

If either of 4(c)i. or 4(c)ii. result in a smaller residual than the previous
best estimate of the parameters, take the $\\{\gamma,K,\omega,e,T_{0},P\\}$
used as the new best estimate

#### 2.3.1. Estimating the period

If $P$ is not provided by the user, it must also be estimated. Both estimates
of the period mentioned in 3(a) and 3(b) begin by finding each point where the
data crosses a given line. The line used here is the $V_{avg}$ axis, also the
estimate for $\gamma$. A point $V_{i}$ is deemed a ‘crossing point’ if
$V_{i-1}$ and $V_{i+1}$ lie on opposite sides of the axis formed by $V_{avg}$.
Thus, unless the data is messy where it crosses the axis or $V_{i}=V_{avg}$,
crossing points will occur in pairs. A common error caused by messy data is
that the data crosses $V_{avg}$ twice but not in a sequential pair, as
demonstrated in Figure 1. It is assumed that there will be sufficient data
that if the crossing points are on opposite sides of an extremum in the true
radial velocity curve, at least one point in between will exceed a
displacement of $\frac{\text{1}}{\text{6}}^{\text{th}}$ the amplitude of the
data, $\frac{V_{max}-V_{min}}{2}$, from $V_{avg}$. So, for all crossing points
that do not occur sequentially in the data, the code verifies that there is at
least one data point in between that meets the above criteria. If this is not
the case, whichever of the two crossing points lay further from the $V_{avg}$
axis is discarded and replaced by the data point directly before (if the first
point is discarded) or after (if the second point is discarded) the point that
was kept.

After finding the crossing points, the first method estimates the period by
finding two crossing points, separated by an integer number of periods, that
have the closest radial velocity values. The period is then found from these
time values and the integer number of periods between the two points.

In the second method, once the crossing points are found, if the points span
over one period the extremum between every pair of points is found. Now $P$
may be found by averaging all estimates of the period from the crossing data
as well as the extrema data.

Now both estimates for $P$ may be used in the iterations of $e$.

Figure 1: Example of two non-consecutive points that cross the $V_{avg}$ axis.
Here the second red point would be kept and the first would be replaced by the
purple point.

### 2.4. Error analysis

Once the parameters $\boldsymbol{p}$ have been determined, the uncertainties
for each of the six (or five if $P$ is known) fitted parameters are computed
and returned as well.

First, the covariance matrix, $\boldsymbol{C}$, is given by inverting the
Hessian matrix, i.e.

$\boldsymbol{C}=\boldsymbol{H}^{-1}$ (18)

The diagonal elements of the covariance matrix, $\boldsymbol{C}_{kk}$, give
the variance of parameters $\boldsymbol{p}_{k}$. Thus, the square root of the
diagonal elements gives the standard error in the associated parameters (Press
et al. 1992). I.e.,

$\sqrt{\boldsymbol{C}_{kk}}=\sigma_{k}$ (19)

The off diagonal elements, $\boldsymbol{C}_{mk}$, where $m\neq k$, gives the
covariance of parameters $\boldsymbol{p}_{k}$ and $\boldsymbol{p}_{m}$. The
linear coefficient of correlation for two given parameters is defined as

$r_{mk}=\frac{\sigma_{mk}}{\sigma_{m}\sigma_{k}}$ (20)

Therefore, the correlation coefficient $r_{mk}$ can be determined from the
covariance matrix by

$r_{mk}=\frac{\boldsymbol{C}_{mk}}{\sqrt{\boldsymbol{C}_{mm}}\sqrt{\boldsymbol{C}_{kk}}}$
(21)

(Press et al. 1992).

## 3\. Comparison of results

Here the BinaryStarSolver code is used to determine the orbital elements of 31
Cygni, 32 Cygni, and Tau Persei using R. F. Griffin’s 2008 and R. E. M.
Griffin’s 1992 data, with the same weights as originally published. The
results are compared to the original solutions to show the effectiveness of
this code.

Table 1: Orbital Elements of 31 Cygni, from Griffin (2008)

Orbital Element | Griffin Solution | Solution, this paper
---|---|---
$P$ (days) | $3784.3$ | $3784.3$
$\gamma$ (km/s) | $-6.421\pm 0.034$ | $-6.425\pm 0.032$
$K$ (km/s) | $13.94\pm 0.04$ | $13.94\pm 0.04$
$\omega$ (°) | $204.5\pm 1.0$ | $204.5\pm 1.0$
$e$ | $0.2084\pm 0.0031$ | $0.2080\pm 0.0029$
$T_{0}$ (RJD) | $52346\pm 9$ | $52346\pm 9$
$a\sin{i}$ (Gm) | $709.5\pm 2.3$ | $709.3\pm 2.1$
$f(m)$ $(\textup{M}_{\odot})$ | $0.995\pm 0.010$ | $0.993\pm 0.009$

Note – Here, the period is kept fixed, as it is in Griffin (2008).

Figure 2: 31 Cygni radial velocity observations of Griffin (2008) plotted as a function of time, with the velocity curve using the parameters found by this code drawn through the points. Table 2: Orbital Elements of 32 Cygni, from Griffin (2008) Orbital Element | Griffin Solution | Solution, this paper
---|---|---
$P$ (days) | $1147.51\pm 0.31$ | $1147.58\pm 0.31$
$\gamma$ (km/s) | $-6.389\pm 0.032$ | $-6.393\pm 0.031$
$K$ (km/s) | $16.77\pm 0.05$ | $16.76\pm 0.05$
$\omega$ (°) | $221.4\pm 0.5$ | $221.5\pm 0.5$
$e$ | $0.3041\pm 0.0027$ | $0.3040\pm 0.0027$
$T_{0}$ (RJD) | $52647.4\pm 1.4$ | $52647.6\pm 1.4$
$a\sin{i}$ (Gm) | $252.0\pm 0.8$ | $252.0\pm 0.8$
$f(m)$ $(\textup{M}_{\odot})$ | $0.486\pm 0.005$ | $0.484\pm 0.005$
Figure 3: 32 Cygni radial velocity observations of Griffin (2008) plotted as a function of time, with the velocity curve using the parameters found by this code drawn through the points. Table 3: Orbital Elements of $\tau$ Persei, from Griffin et al. (1992) Orbital Element | Griffin et al. Solution | Solution, this paper
---|---|---
$P$ (days) | $1516.1\pm 1.9$ | $1517.0\pm 2.6$
$\gamma$ (km/s) | $2.37\pm 0.11$ | $2.34\pm 0.15$
$K$ (km/s) | $18.7\pm 0.3$ | $18.7\pm 0.4$
$\omega$ (°) | $234.7\pm 1.3$ | $233.7\pm 1.5$
$e$ | $0.721\pm 0.007$ | $0.719\pm 0.009$
$T_{0}$ (RJD) | $47524.1\pm 1.4$ | $46006.7\pm 3.4$
$a\sin{i}$ (Gm) | $271\pm 5$ | $270\pm 7$
$f(m)$ $(\textup{M}_{\odot})$ | $0.344\pm 0.020$ | $0.342\pm 0.025$
Figure 4: $\tau$ Persei radial velocity observations of Griffin et al. (1992)
plotted as a function of time, with the velocity curve using the parameters
found by this code drawn through the points. This star, with its highly
eccentric ($e=0.72$) orbit, provides a challenging test of convergence for the
BinaryStarSolver code.

## 4\. Conclusions

In this paper, we present a new, open source Python code BinaryStarSolver to
determine spectroscopic binary orbits, that is readily accesible to both
amateur and professional astronomers. The determination of spectroscopic
binary orbits and stellar masses is an example of another important area of
astrophysics to which amateurs can now make significant contributions, but
analysis software accessible to the amateur community remains a constraint.
The present code is intended to provide a robust, easy-to-use package to
determine binary star orbits given a time series of radial velocity
observations.

We describe the solution procedure used by the code to determine binary star
orbital elements, and present results for three binary systems with recent
analyses in the literature. One of these binaries ($\tau$ Per) has a highly
eccentric orbit ($e=0.72$), which renders the problem less stable numerically.
Our orbital solutions for all three binaries are in extremely close agreement
with the previously published orbits, demonstrating the accuracy and
robustness of the method used by the BinaryStarSolver code.

The BinaryStarSolver code described in this paper is available, with complete
documentation, from
GitHub111https://github.com/NickMilsonPhysics/BinaryStarSolver. We also plan
to have the code included on the American Association of Variable Stars
Observers (AAVSO) software directory page.

##  Appendix A Appendix: Summary of Equations

The Hessian matrix222If $P$ is known, the Hessian reduces to a 5x5 matrix
where the 6th row and 6th column are deleted.:

$\boldsymbol{H}=\begin{bmatrix}\frac{\partial^{2}\chi^{2}}{\partial\gamma^{2}}&\frac{\partial^{2}\chi^{2}}{\partial\gamma\partial
K}&\frac{\partial^{2}\chi^{2}}{\partial\gamma\partial\omega}&\frac{\partial^{2}\chi^{2}}{\partial\gamma\partial
e}&\frac{\partial^{2}\chi^{2}}{\partial\gamma\partial
T_{0}}&\frac{\partial^{2}\chi^{2}}{\partial\gamma\partial P}\\\
\frac{\partial^{2}\chi^{2}}{\partial
K\partial\gamma}&\frac{\partial^{2}\chi^{2}}{\partial
K^{2}}&\frac{\partial^{2}\chi^{2}}{\partial
K\partial\omega}&\frac{\partial^{2}\chi^{2}}{\partial K\partial
e}&\frac{\partial^{2}\chi^{2}}{\partial K\partial
T_{0}}&\frac{\partial^{2}\chi^{2}}{\partial K\partial P}\\\
\frac{\partial^{2}\chi^{2}}{\partial\omega\partial\gamma}&\frac{\partial^{2}\chi^{2}}{\partial\omega\partial
K}&\frac{\partial^{2}\chi^{2}}{\partial\omega^{2}}&\frac{\partial^{2}\chi^{2}}{\partial\omega\partial
e}&\frac{\partial^{2}\chi^{2}}{\partial\omega\partial
T_{0}}&\frac{\partial^{2}\chi^{2}}{\partial\omega\partial P}\\\
\frac{\partial^{2}\chi^{2}}{\partial
e\partial\gamma}&\frac{\partial^{2}\chi^{2}}{\partial e\partial
K}&\frac{\partial^{2}\chi^{2}}{\partial
e\partial\omega}&\frac{\partial^{2}\chi^{2}}{\partial
e^{2}}&\frac{\partial^{2}\chi^{2}}{\partial e\partial
T_{0}}&\frac{\partial^{2}\chi^{2}}{\partial e\partial P}\\\
\frac{\partial^{2}\chi^{2}}{\partial
T_{0}\partial\gamma}&\frac{\partial^{2}\chi^{2}}{\partial T_{0}\partial
K}&\frac{\partial^{2}\chi^{2}}{\partial
T_{0}\partial\omega}&\frac{\partial^{2}\chi^{2}}{\partial T_{0}\partial
e}&\frac{\partial^{2}\chi^{2}}{\partial
T_{0}^{2}}&\frac{\partial^{2}\chi^{2}}{\partial T_{0}\partial P}\\\
\frac{\partial^{2}\chi^{2}}{\partial
P\partial\gamma}&\frac{\partial^{2}\chi^{2}}{\partial P\partial
K}&\frac{\partial^{2}\chi^{2}}{\partial
P\partial\omega}&\frac{\partial^{2}\chi^{2}}{\partial P\partial
e}&\frac{\partial^{2}\chi^{2}}{\partial P\partial
T_{0}}&\frac{\partial^{2}\chi^{2}}{\partial P^{2}}\end{bmatrix}$

All second order partial derivatives are found numerically using central
differences, from the following analytic first order derivatives

$\frac{\partial\chi^{2}}{\partial\gamma}=2\sum_{i=1}^{N}{[V(t_{i})-V_{i}]}\,w_{i}$
$\frac{\partial\chi^{2}}{\partial
K}=2\sum_{i=1}^{N}{[V(t_{i})-V_{i}][\cos(v+\omega)+e\cos\omega]}\,w_{i}$
$\frac{\partial\chi^{2}}{\partial\omega}=-2K\sum_{i=1}^{N}{[V(t_{i})-V_{i}][\sin(v+\omega)+e\sin\omega]}\,w_{i}$
$\frac{\partial\chi^{2}}{\partial
e}=2K\sum_{i=1}^{N}{[V(t_{i})-V_{i}][\cos(v+\omega)-\sin(v+\omega)\frac{\partial
v}{\partial e}]}\,w_{i}$ $\frac{\partial\chi^{2}}{\partial
T_{0}}=-2K\sum_{i=1}^{N}{[V(t_{i})-V_{i}]\sin(v+\omega)\frac{\partial
v}{\partial T_{0}}}\,w_{i}$ $\frac{\partial\chi^{2}}{\partial
P}=-K\sum_{i=1}^{N}{[V(t_{i})-V_{i}]\sin(v+\omega)\frac{\partial v}{\partial
P}}\,w_{i}$

Here $\frac{\partial v}{\partial e}$, $\frac{\partial v}{\partial T_{0}}$, and
$\frac{\partial v}{\partial P}$ are given by

$\begin{split}\frac{\partial v}{\partial
e}=2\cos^{2}\Big{(}\frac{v}{2}\Big{)}\Bigg{[}\sqrt{\frac{1-e}{1+e}}\frac{1}{(1-e)^{2}}\tan\Big{(}\frac{E}{2}\Big{)}\\\
+\frac{1}{2}\sec^{2}\Big{(}\frac{E}{2}\Big{)}\sqrt{\frac{1+e}{1-e}}\frac{\sin{E}}{(1-e\cos{E})}\Bigg{]}\end{split}$
$\frac{\partial v}{\partial
T_{0}}=\frac{1+\cos{v}}{1+\cos{E}}\sqrt{\frac{1+e}{1-e}}\frac{-n}{(1-e\cos{E})}$
$\frac{\partial v}{\partial
P}=-\frac{2\pi}{P^{2}}\frac{1+\cos{v}}{1+\cos{E}}\sqrt{\frac{1+e}{1-e}}\frac{t_{i}-T_{0}}{1-e\cos{E}}$

## References

  * 1 Gavin, H. P. 2011, The Levenberg-Marquardt algorithm for nonlinear least squares curve-fitting problems, Dept. of Civil and Environmental Engineering, Duke Univ. Tech. Rep., 28
  * 2 Green, R. M. 1985, Spherical Astronomy. Cambridge Univ. P., Cambridge
  * 3 Griffin, R. E. M., Schröder, K.-P., Misch, A., and Griffin, R. F. 1992, Astron. Astrophys., 254, 289
  * 4 Griffin, R. F. 2008, The Observatory, 128, 362
  * 5 Marquardt, D. W. 1963, Journal of the Society for Industrial and Applied Mathematics, 11(2), 431
  * 6 Pollmann, E. & Bennett, P. D. 2020, J. Amer. Assoc. Var. Star Obs., 48, 2
  * 7 Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. 1992, Numerical Recipes in Fortran: The Art of Scientific Computing, 2nd ed., Cambridge Univ. Press, Cambridge
  * 8 Transtrum, M. K., and Sethna, J. P. 2012, Improvements to the Levenberg-Marquardt algorithm for nonlinear least-squares minimization, arXiv:1201.5885