# Generating Realistic Two-Line Elements for Notional Space Vehicles and
Constellations

Troy Rockwood, Greg Steeger, Matthew Stein

###### Abstract

As space becomes increasingly populated with new satellites and systems,
modeling and simulating existing and future systems becomes more important.
The two-line element set has been a standard format for sharing data about a
satellite’s orbit since the 1960s, and well-developed algorithms can predict
the future location of satellites based on this data. In order to simulate
potential future systems, especially when mixed with existing systems, data
must be generated to represent the desired orbits. We present a means to
create two-line element sets with parameters that closely resemble real
satellite behavior, and rely on a novel approach to calculate the mean motion
for even greater accuracy.

##  1 Introduction

The population of active satellites in orbit around Earth has increased
exponentially over the last few decades[1]. Commercial, civil, and government
entities continue to expand their role and influence in space. Modeling and
simulation of satellites and their behavior has likewise increased. The use
cases include constellation planning, optimization, performance assessment,
and conjunction analysis[2, 3, 4].

There are many software tools and algorithms to simulate satellites and
constellations of satellites. Among these is the well-known simplified
perturbations model, often referred to as SGP4[5]. These algorithms make use
of two-line element sets (TLEs)[6, 7] to describe the orbit of each satellite
and begin the simulation. Regardless of the algorithm employed, TLEs are
widely used and published for most active and inactive objects in orbit.

For the purposes of planning future space systems, it is useful to simulate
satellites and orbits that do not yet exist, possibly alongside real systems
already in orbit. Doing so requires the generation of TLEs for the notional
systems. A naïve approach based only on classical orbital elements might look
like the following:

  1. 1.

The desired inclination, eccentricity, altitude at perigee, and argument of
perigee are chosen according to the desired orbit

  2. 2.

The first and second derivative of mean motion and drag terms are set to
exactly zero

  3. 3.

The mean motion is calculated from Kepler’s third law using the altitude at
perigee and a mean radius of the Earth of 6,371 km

  4. 4.

The rest of the TLE elements are freely chosen, and the checksum is computed
for the last character in each line

The classical orbital elements in Step 1 are freely chosen, but Step 2 and 3
each present opportunities for significant error in the desired simulation
results. When propagating TLEs, some or all of the terms in Step 2 help model
realistic orbital decays depending on the algorithm chosen. Simulated notional
TLEs may thus behave quite differently than TLEs created for orbiting
satellites, making mixed simulations less realistic especially as simulation
time increases. Similarly, calculating mean motion from a fixed value for the
radius of the Earth presents potential for significant error in drift. Because
the Earth is not a perfect sphere, gravitational effects experienced by
satellites will be a function of the intended orbit, namely the inclination,
eccentricity, and argument of perigee.

In order to realistically simulate the behavior of satellites, the generated
TLEs should contain terms that accurately resemble those of real systems. We
describe a novel method for generating TLEs for almost any notional system
including realistic terms for drag and the first and second derivative of mean
motion. Additionally, we describe a novel way to determine the mean motion
(MM) for the desired orbit which is a function of inclination, eccentricity,
and the argument of perigee. Finally, we summarize our approach which is an
expansion of the above steps in Section 6, including a polynomial approach to
determining the mean radius of the Earth for computing the mean motion.

##  2 Essential TLE Terms

This study ignores the title line sometimes included in some TLEs, and focuses
only on lines one and two. A detailed visualization and explanation of each
TLE term is illustrated in Figure 1 of Ref [6]. The satellite catalog number,
classification, international designator value, element set number, and
revolution number can all be chosen at will, and do not affect the simulation.
The only exception is that satellite catalog numbers should not conflict with
other satellites in the simulation. Additionally, the ephemeris type should be
0.

The epoch of the TLE can be chosen at will, but in generating an epoch one
should consider the following implications. If a simulation only involves
notional systems, then it is only important that all TLEs in the simulation
have a similar epoch. If mixing notional with real satellites, it is best to
choose an epoch close to the real systems of interest.

The inclination, right ascension of the ascending node (RAAN), eccentricity,
argument of perigee, and mean anomaly are all chosen according to the desired
simulation. For example, a notional low Earth orbit constellation might
consist of 225 satellites in 15 planes with 15 satellites per plane. The
inclination might be chosen to be 90.0° to model a polar orbit. To provide
somewhat even distribution of coverage around the equator, the RAAN of each
plane could range from 0.0° to 168.0° in steps of 12.0°. Within the planes,
the mean anomaly might likewise be spaced from 0.0° to 336.0° in steps of
24.0°.

##  3 Catalog-Derived Terms

Three of the remaining terms–the first and second derivative of the mean
motion, and the drag term–can be derived from existing TLEs which are
accessible in open repositories on the internet. We use a collection of 18,477
published TLEs associated with satellites, rocket bodies, and other debris. In
order to assign terms to notional TLEs that represent realistic TLEs, the TLEs
are categorized into low-Earth orbit (LEO), medium Earth orbit (MEO),
geosynchronous (GEO), and highly elliptical orbit (HEO) populations.

The algorithm to assign TLEs to a population is given in Algorithm 1 and the
impact of the included criteria is illustrated in Figure 1. There are various
other criteria and algorithms, and the methods in this paper could be easily
adapted by modifying Algorithm 1, including additional population categories.

Algorithm 1 TLE Classification

if Eccentricity $\geq$ 0.5 then

HEO

else if Mean Motion $\geq$ 11.25 then

LEO

else if Mean Motion $\geq$ 1.2 then

MEO

else

GEO

Figure 1: The categorization of 18,477 published TLEs using the criteria from
Algorithm 1.

With populations assigned, a mean value can be computed for the first and
second derivative of the mean motion, and the drag term as shown in Table 1
and Figure 2.

Figure 2: The first and second derivative of the mean motion, and the drag term, highlighted by population type for all 18,477 TLEs in this study. Population | MM${}^{\prime}\,\times 10^{3}~{}(\textrm{Rev}/\textrm{day}^{2})$ | MM${}^{\prime\prime}\,\times 10^{6}~{}(\textrm{Rev}/\textrm{day}^{3})$ | Drag Term$\,\times 10^{3}~{}\left({R}_{\bigoplus}^{-1}\right)$
---|---|---|---
HEO | $0.048575\pm 0.384281$ | $~{}~{}0.0125888\pm 0.120716$ | $1.558450\pm 27.4125$
LEO | $0.154256\pm 3.778700$ | $~{}~{}0.0942242\pm 2.866240$ | $0.377655\pm 10.7492$
MEO | $0.154986\pm 1.640110$ | $-0.0166109\pm 0.381907$ | $1.295840\pm 48.6634$
GEO | $0.001190\pm 0.001467$ | $~{}~{}~{}~{}0\pm 0$ | $0.639138\pm 25.2697$
Table 1: The mean and standard deviation values for the first and second
derivative of the mean motion (MM′ and MM′′, respectively), and the drag term,
grouped by population.

##  4 Determining Mean Motion

A simple modification of Kepler’s third law allows the calculation of the mean
motion in revolutions per day from the semi-major axis of orbit, as shown in
Equation 1:

$\textrm{MM}=\frac{T_{\bigoplus}}{2\pi}\sqrt{\frac{\mu}{a^{3}}}$ (1)

where $T_{\bigoplus}$ is a solar day (86,400 seconds), $a$ is the semi-major
axis, and $\mu$ is the standard gravitational parameter of Earth
($3.986004418\times 10^{14}\textrm{m}^{3}/\textrm{s}^{2}$). The semi-major
axis can be computed from the radius of the Earth $R_{\bigoplus}$, the
altitude at perigee $A_{p}$, and the eccentricity $e$ as shown in Equation 2.

$a=\frac{A_{p}+R_{\bigoplus}}{1-e}$ (2)

The Earth’s radius $R_{\bigoplus}$ is the distance from the center of the
Earth to a point on its surface. Because the Earth is not a perfect sphere and
due to the Earth’s topography, this radius is not constant across the Earth’s
surface. The Earth can be modeled as an oblate spheroid with a minimum radius
of 6,357 km at the poles, and a maximum radius of 6,378 km at the equator. For
reference, analysis commonly relies on a nominal or “globally averaged” value
of 6,371 km[8, 9]. Globally averaged values are computed by either

  * •

Taking the mean of three radii measurements, two from the equator and one from
a pole;

  * •

Using the radius of a sphere with the same surface area as the Earth’s
(authalic radius); or

  * •

Using the radius of a sphere that has the same volume as the Earth spheroid
(volumetric radius).

While it is possible to simply use a single value for the radius of the Earth
such as a globally averaged value, this will present error in almost any
simulation as will be quantified below. Instead, the value used for
$R_{\bigoplus}$ in Equation 2 should be relevant to the mean gravitational
force that the satellite experiences as it travels around the non-spherical
Earth. In other words, the value for $R_{\bigoplus}$ should be the mean of the
distance between the center of the Earth ellipsoid and the surface of the
Earth directly under a satellite as it travels in its orbit. This study will
refer to that value as the mean radius under the satellite, or
$R_{\textrm{s}}$. Utilizing accurate numbers for this measure is imperative to
generating realistic MMs, and consequently realistic TLEs. The following
sections detail an approach to determine $R_{\textrm{s}}$ and provides values
that can be used in future simulations.

##  5 Mean Radius

This analysis assumes that the Earth is modeled as a WGS-84 ellipsoid[10]. The
Earth ellipsoid also has variation in its gravitational potential as defined
in the Earth Gravitational Model 2008[11] but this study focuses on SGP4
models that include the J2, J3, and J4 zonal harmonics coefficients. This
study relies on the use the Satrec package for python[12]. Satrec uses a
satellite’s TLE to compute the satellite’s position and velocity for a given
date and time.

To determine $R_{\textrm{s}}$, a satellite is propagated around the Earth and
the geodetic latitude directly under the satellite is determined for many
equally-spaced time steps throughout a simulation. Given each latitude, the
radius of the ellipsoid at that position is calculated and stored in an array.
Finally, the mean of these measurements over a simulation yields
$R_{\textrm{s}}$. In the following sections, we describe a method for
determining latitude directly from SGP4 propagation, and then our results for
computing $R_{\textrm{s}}$ from a large set of input inclinations,
eccentricities, and arguments of perigee.

###  5.1 Determining Radii from Propagation

Given Cartesian coordinates for a satellite, we can convert them into latitude
($\phi$), longitude ($\lambda$), and height ($h$). SGP4 propagation provides
the ECEF coordinates for a satellite at all desired time steps in the
simulation. The longitude is easily computed from the Cartesian coordinates as

$\lambda=ArcTan2\left(y,x\right)$ (3)

where ArcTan2 is the function that returns the angle from the origin, knowing
which quadrant the position $(x,y)$ lies in. The latitude and height can be
computed by application of Ferrari’s solution[13] which is not covered here.

For a given latitude $\phi$, the radius of the ellipsoid at that latitude
$R_{\phi}$ is given by

$R_{\phi}=\sqrt{\frac{\left(R_{\textrm{e}}^{2}cos(\phi)\right)^{2}+\left(R_{\textrm{p}}^{2}sin(\phi)\right)^{2}}{\left(R_{\textrm{e}}cos(\phi)\right)^{2}+\left(R_{\textrm{p}}sin(\phi)\right)^{2}}}$
(4)

For every time step in the simulation, we determine $\phi$ and therefore a
value for $R_{\phi}$. We then take the mean of all values for $R_{\phi}$ to
determine $R_{\textrm{s}}$.

###  5.2 Results for $R_{\textrm{s}}$

The results in this paper are presented for orbits with inclinations varying
from 0 to 90°, eccentricities ranging from 0 to 0.999, and arguments of
perigee ranging from 0 to 90°. The choice of the altitude at perigee does not
affect the results as long it creates a stable orbit. This analysis considered
an altitude at perigee of 605.736 km for all runs.

For each choice of inclination, eccentricity, and argument of perigee, we
determine $R_{\textrm{s}}$ from a single orbit that is propagated with 1,000
equally spaced time steps. The number of orbits and steps per orbit can be
increased but do not have a significant effect on the results. Furthermore,
significantly increasing the number of orbits can introduce unwanted error due
to the potentially degraded orbit modeled in SGP4 propagation.

The results provide a variation in $R_{\textrm{s}}$ from 6358.669 km to
6378.137 km, partly illustrated in Figure 3. When computing MM, using
$R_{\textrm{s}}$ presents up to a ${\sim}1.5\%$ correction as compared to
using the typical WGS-84 mean radius of 6,371 km.To measure the impact of such
a correction, TLEs are created using both the WGS-84 mean radius of the Earth
and $R_{\textrm{s}}$. This results in two TLEs that are identical except for
the MM. We repeat this process for the same ranges of inclinations,
eccentricities, and arguments of perigee. The time for each pair of TLEs to
drift 100 km apart is measured and illustrated in Figure 4.

Figure 3: Top: The mean radius of the Earth, under a satellite, as it
propagates throughout its orbit. The red contour represents the commonly used
mean radius of the Earth of 6,371.009 km. Bottom: The total percent error
between a standard WGS-84 mean Earth radius of 6371.009 km and the method
described in Section 5. The red contour represents where the difference is
exactly zero. These contour plots were created with an argument of perigee of
0°.  Figure 4: The time for two TLEs to drift 100 km, each with different
assumptions for the Earth radius: one with a standard WGS-84 mean Earth
radius, and one with an $R_{\textrm{s}}$ value computed from Section 5. This
contour plot was created with an argument of perigee of 0°.

In many cases, this drift can be greater than 100 km in less than only a few
minutes. For many analyses such as conjunction analysis, object tracking, and
coverage metrics, this accumulated error can present significant changes to
the simulation results, especially for longer duration simulations. The method
described above for determining $R_{\textrm{s}}$ presents an improved means of
positioning and modeling satellites through longer simulations.

##  6 Assembling Terms for Notional TLEs

With the catalog-derived terms ready, realistic TLEs can be created based on
the desired inclination, eccentricity, argument of perigee, RAAN, and altitude
at perigee. Table 2 lists all TLE terms, and whether they are derived from a
catalog, determined by calculation, or taken as direct inputs.

Term | Method
---|---
First Derivative of MM | Catalog
Second Derivative of MM | Catalog
Drag (B*) | Catalog
Inclination | Free
RAAN | Free
Eccentricity | Free
Argument of Perigee | Free
Mean Anomaly | Free
Mean Motion | Determined from $R_{\textrm{s}}$
Table 2: The TLE terms used in this study, and whether they are catalog-
derived, free variables, or determined from other elements including the
calculation of $R_{\textrm{s}}$ as described in Section 5. All other TLE terms
are trivially chosen except for the checksum which is computed based on all
other entries in each line.

The procedure for creating a TLE is as follows:

  1. 1.

The desired inclination, eccentricity, altitude at perigee, and argument of
perigee are chosen according to the desired orbit

  2. 2.

The $R_{\textrm{s}}$ value is determined based inclination, eccentricity, and
argument of perigee

  3. 3.

The semi-major axis is determined from Equation 2, substituting
$R_{\textrm{s}}$ in place of $R_{\bigoplus}$

  4. 4.

The MM is calculated from Equation 1

  5. 5.

The satellite population is classified according to Algorithm 1 using the
calculated MM and intended eccentricity

  6. 6.

The first and second derivative of mean motion and drag terms are assigned
from Table 1

  7. 7.

The rest of the TLE elements are freely chosen, and the checksum is computed
for the last character in each line

The simulated data for $R_{\textrm{s}}$ described in Section 5 can be fit by a
polynomial function whose parameters can be numerically determined. We fit a
fifth-order 3-dimensional polynomial to the data for $R_{s}$ with 56
parameters and determined a maximum difference of $5.47\times 10^{-3}$ %
compared to directly simulating the orbit. This is across all inclinations,
eccentricities, and arguments of perigee. Similarly, for an eighth-order
polynomial with 165 parameters, the maximum difference was $0.93\times
10^{-3}$ % compared to direct simulation, across all values. The values for
the parameters of the fifth-order polynomial are shown in Table 3. The
polynomial with associated parameters is shown in Equation 5 with $i$, $e$,
and $\omega$ the inclination, eccentricity, and argument of perigee,
respectively.

$\begin{split}\textrm{R}_{\textrm{s}}\left(i,e,\omega\right)=&c_{000}+c_{001}\omega+c_{002}\omega^{2}+c_{003}\omega^{3}+c_{004}\omega^{4}+\\\
c_{005}\omega^{5}+&c_{010}e+c_{011}{\omega}e+c_{012}{\omega}^{2}e+c_{013}{\omega}^{3}e+\\\
c_{014}{\omega}^{4}e+&c_{020}e^{2}+c_{021}{\omega}e^{2}+c_{022}{\omega}^{2}e^{2}+c_{023}{\omega}^{3}e^{2}+\\\
c_{030}e^{3}+&c_{031}{\omega}e^{3}+c_{032}{\omega}^{2}e^{3}+c_{040}e^{4}+c_{041}{\omega}e^{4}+\\\
c_{050}e^{5}+&c_{100}i+c_{101}i{\omega}+c_{102}i{\omega}^{2}+c_{103}i{\omega}^{3}+\\\
c_{104}i{\omega}^{4}+&c_{110}ie+c_{111}ie{\omega}+c_{112}ie{\omega}^{2}+c_{113}ie{\omega}^{3}+\\\
c_{120}ie^{2}+&c_{121}ie^{2}{\omega}+c_{122}ie^{2}{\omega}^{2}+c_{130}ie^{3}+\\\
c_{131}ie^{3}{\omega}+&c_{140}ie^{4}+c_{200}i^{2}+c_{201}i^{2}{\omega}+c_{202}i^{2}{\omega}^{2}+\\\
c_{203}i^{2}{\omega}^{3}+&c_{210}i^{2}e+c_{211}i^{2}e{\omega}+c_{212}i^{2}e{\omega}^{2}+\\\
c_{220}i^{2}e^{2}+&c_{221}i^{2}e^{2}{\omega}+c_{230}i^{2}e^{3}+c_{300}i^{3}+c_{301}i^{3}{\omega}+\\\
c_{302}i^{3}{\omega}^{2}+&c_{310}i^{3}e+c_{311}i^{3}e{\omega}+c_{320}i^{3}e^{2}+c_{400}i^{4}+\\\
c_{401}i^{4}{\omega}+&c_{410}i^{4}e+c_{500}i^{5}\end{split}$ (5)

Term | Value | Term | Value
---|---|---|---
$c_{000}$ | $6.377788\times 10^{6}$ | $c_{112}$ | $-9.972427\times 10^{-2}$
$c_{001}$ | $2.836106\times 10^{1}$ | $c_{113}$ | $7.383957\times 10^{-4}$
$c_{002}$ | $-4.915346\times 10^{-1}$ | $c_{120}$ | $9.150674\times 10^{1}$
$c_{003}$ | $-4.017214\times 10^{-3}$ | $c_{121}$ | $-1.33747\times 10^{0}$
$c_{004}$ | $1.273759\times 10^{-4}$ | $c_{122}$ | $1.155062\times 10^{-3}$
$c_{005}$ | $-5.643483\times 10^{-7}$ | $c_{130}$ | $4.645793\times 10^{0}$
$c_{010}$ | $2.120874\times 10^{3}$ | $c_{131}$ | $-1.909859\times 10^{0}$
$c_{011}$ | $-1.886565\times 10^{2}$ | $c_{140}$ | $5.290634\times 10^{1}$
$c_{012}$ | $5.241726\times 10^{0}$ | $c_{200}$ | $-3.926748\times 10^{0}$
$c_{013}$ | $-3.84173\times 10^{-2}$ | $c_{201}$ | $1.703196\times 10^{-2}$
$c_{014}$ | $-2.268633\times 10^{-6}$ | $c_{202}$ | $-8.207037\times 10^{-6}$
$c_{020}$ | $-3.447518\times 10^{3}$ | $c_{203}$ | $2.674298\times 10^{-8}$
$c_{021}$ | $8.553227\times 10^{1}$ | $c_{210}$ | $4.368207\times 10^{0}$
$c_{022}$ | $-4.719285\times 10^{0}$ | $c_{211}$ | $-9.744218\times 10^{-2}$
$c_{023}$ | $3.438466\times 10^{-2}$ | $c_{212}$ | $-8.298287\times 10^{-6}$
$c_{030}$ | $9.952275\times 10^{3}$ | $c_{220}$ | $2.116514\times 10^{-1}$
$c_{031}$ | $1.634122\times 10^{2}$ | $c_{221}$ | $-4.547244\times 10^{-4}$
$c_{032}$ | $5.259614\times 10^{-2}$ | $c_{230}$ | $-9.768337\times 10^{-2}$
$c_{040}$ | $-1.695787\times 10^{4}$ | $c_{300}$ | $1.323687\times 10^{-3}$
$c_{041}$ | $-8.008678\times 10^{1}$ | $c_{301}$ | $-1.327782\times 10^{-4}$
$c_{050}$ | $8.501297\times 10^{3}$ | $c_{302}$ | $5.032226\times 10^{-8}$
$c_{100}$ | $2.948006\times 10^{1}$ | $c_{310}$ | $-3.319618\times 10^{-2}$
$c_{101}$ | $-1.151884\times 10^{0}$ | $c_{311}$ | $7.264131\times 10^{-4}$
$c_{102}$ | $1.653916\times 10^{-2}$ | $c_{320}$ | $-5.846808\times 10^{-4}$
$c_{103}$ | $-1.276845\times 10^{-4}$ | $c_{400}$ | $4.655915\times 10^{-4}$
$c_{104}$ | $3.521573\times 10^{-}8$ | $c_{401}$ | $6.438237\times 10^{-8}$
$c_{110}$ | $-1.676297\times 10^{2}$ | $c_{410}$ | $4.999741\times 10^{-6}$
$c_{111}$ | $6.645437\times 10^{0}$ | $c_{500}$ |  $-2.08623\times 10^{-6}$
Table 3: The parameters for a fifth-order polynomial fit of the
$R_{\textrm{s}}$ data across all inclinations, eccentricities, and arguments
of perigee from Section 5.

##  7 Conclusion

We present a new means of creating TLEs for notional satellites with realistic
terms for mean motion, its first and second derivatives, and the drag term.
The only inputs required are the desired inclination, eccentricity, argument
of perigee, and altitude at perigee. The produced TLEs can be used in
simulation alone or alongside TLEs from real satellites, and will behave like
real systems. This increases the overall simulation accuracy for performance
evaluation, conjunction analysis, and more. The methods described can be
expanded to other classification schemes for different sets of populations.

Future work will involve investigating the satellite catalog and possibly
deriving correlations among populations that could reduce the standard
deviations in Table 1, though this may mean significantly increasing the
complexity of population categorization. This work could also be adapted to
include the EGM2008 terms and also regenerated for different satellite
classification schemes[11].

## Acknowledgments

This technical data was produced for the U. S. Government under Contract No.
FA8702-22-C-0001, and is subject to the Rights in Technical Data-Noncommercial
Items Clause DFARS 252.227-7013 (FEB 2014). © 2022 The MITRE Corporation. All
Rights Reserved. DISTRIBUTION A: Approved for public release; distribution
unlimited. Public Release Case Number 22-0304.

## References

  * UCS [2021] Ucs satellite database, May 2021.  URL https://www.ucsusa.org/resources/satellite-database.  https://www.ucsusa.org/resources/satellite-database.
  * Parish [2004] Jason A. Parish.  Optimizing coverage and revisit time in sparse military satellite constellations; a comparison of traditional approaches and genetic algorithms, 2004.  URL https://calhoun.nps.edu/handle/10945/1209.
  * Hitomi and Selva [2018] Nozomi Hitomi and Daniel Selva.  Constellation optimization using an evolutionary algorithm with a variable-length chromosome.  In _2018 IEEE Aerospace Conference_ , pages 1–12, 2018.  doi: 10.1109/AERO.2018.8396743.
  * Budianto and Olds [2004] Irene A. Budianto and John R. Olds.  Design and deployment of a satellite constellation using collaborative optimization.  _Journal of Spacecraft and Rockets_ , 41(6):956–963, 2004.  doi: 10.2514/1.14254.  URL https://doi.org/10.2514/1.14254.
  * Vallado et al. [2012] David Vallado, Paul Crawford, Richard Hujsak, and T.S. Kelso.  Revisiting spacetrack report #3.  _AIAA/AAS Astrodynamics Specialist Conference and Exhibit_ , 2012\.  doi: 10.2514/6.2006-6753.  URL https://arc.aiaa.org/doi/abs/10.2514/6.2006-6753.
  * Kelso [2019] T S Kelso.  Norad two-line element set format, Dec 2019.  URL https://www.celestrak.com/NORAD/documentation/tle-fmt.php.  https://www.celestrak.com/NORAD/documentation/tle-fmt.php.
  * Vallado and Cefola [2012] David Vallado and Paul Cefola.  Two-line element sets - practice and use.  _Proceedings of the International Astronautical Congress, IAC_ , 7:5812–5825, 01 2012.
  * Mamajek et al. [2015] EE Mamajek, A Prsa, G Torres, P Harmanec, M Asplund, PD Bennett, N Capitaine, J Christensen-Dalsgaard, E Depagne, WM Folkner, et al.  Iau 2015 resolution b3 on recommended nominal conversion constants for selected solar and planetary properties, arxiv151007674 astro-ph, 2015.
  * Moritz [1980] Helmut Moritz.  Geodetic reference system 1980.  _Bulletin géodésique_ , 54(3):395–405, 1980.
  * McClain and Vallado [2001] W.D. McClain and D.A. Vallado.  _Fundamentals of Astrodynamics and Applications_.  Space Technology Library. Springer Netherlands, 2001.  ISBN 9780792369035.  URL https://books.google.com/books?id=PJLlWzMBKjkC.
  * Pavlis et al. [2012] Nikolaos K. Pavlis, Simon A. Holmes, Steve C. Kenyon, and John K. Factor.  The development and evaluation of the earth gravitational model 2008 (egm2008).  _Journal of Geophysical Research: Solid Earth_ , 117(B4), 2012.  doi: https://doi.org/10.1029/2011JB008916.  URL https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011JB008916.
  * Rhodes [2021] Brandon Rhodes.  python-sgp4.  https://github.com/brandon-rhodes/python-sgp4, 2021.
  * Zhu [1994] J. Zhu.  Conversion of earth-centered earth-fixed coordinates to geodetic coordinates.  _IEEE Transactions on Aerospace and Electronic Systems_ , 30(3):957–961, 1994.  doi: 10.1109/7.303772.


# Quench dynamics of mass-imbalanced three-body fermionic systems in a
spherical trap

A. D. Kerin  A. M. Martin  School of Physics, University of Melbourne,
Parkville, VIC 3010, Australia

###### Abstract

We consider a system of two identical fermions of general mass interacting
with a third distinguishable particle via a contact interaction within an
isotropic three-dimensional harmonic trap. We calculate time-dependent
observables of the system after it is quenched in s-wave scattering length. To
do this we use exact closed form mass-imbalanced hyperspherical solutions to
the static three-body problem. These exact solutions enable us to calculate
two time-dependent observables, the Ramsey signal and particle separation,
after the system undergoes a quench from non-interacting to the unitary regime
or vice-versa.

##  I Introduction

Investigating the dynamics of non-equilibrium quantum systems is relevant to
many areas in condensed matter physics. Areas of interest include the study of
superfluid turbulence, topological states, mesoscopic circuits and can even be
extended to neutron star dynamics. In these systems there are many interacting
entities, making exact calculations problematic. However, it is possible to
theoretically and experimentally examine quantum dynamics when there are only
a few quantum bodies in the system of interest. An idealised example of such a
system, where few-body quantum dynamics can be studied, is harmonically
trapped quantum gases, which can be constructed at the single to few atom
limit Serwane  _et al._ (2011); Murmann  _et al._ (2015); Zürn  _et al._
(2013, 2012); Stöferle  _et al._ (2006). In this work we focus on the
interaction quench dynamics of three fermionic atoms trapped in a spherically
symmetric harmonic trap.

We consider a system of two identical fermions which interact with a third
distinguishable particle, which can have a different mass. We consider the
scenario where the contact interaction, between the third particle and the two
identical fermions, is quenched either from the non-interacting regime to the
unitary regime or vice-versa. To do this we utilise the exact solutions for
the system Cui (2012); D’Incao  _et al._ (2018); Jonsell  _et al._ (2002);
Blume and Daily (2010); Kerin and Martin (2022); Kestner and Duan (2007).
These exact solutions have previously been used to elucidate the thermodynamic
properties of quantum gases Liu  _et al._ (2009, 2010); Stöferle  _et al._
(2006); Cui (2012); Rakshit  _et al._ (2012); Kaplan and Sun (2011); Mulkerin
_et al._ (2012a, b); Nascimbène  _et al._ (2010); Ku  _et al._ (2012);
Levinsen  _et al._ (2017); Daily and Blume (2010). Additionally there have
been similar calculations carried out for a two-dimensional system, Ref.
Bougas  _et al._ (2022) considers a system quenched from non-interacting to
various interaction strengths.

The paper is structured as follows. In Sec. II we briefly review the
stationary three-body problem of two identical fermions interacting with a
third distinguishable particle in a spherically symmetric harmonic trap. In
Sec. III we use the eigenstates of the system to investigate the quench
dynamics. In particular we focus on quenches from the non-interacting regime
to the strongly interacting (unitary) regime, a forwards quench, and vice-
versa, a backwards or reverse quench. For these two quenches we evaluate the
Ramsey signal, i.e. the overlap of the time evolving state with the initial
state, and post-quench evolution of the particle separation, as defined by the
hyperradius. For the Ramsey signal we find that it can be calculated semi-
analytically for any initial state for both quenches whilst for the particle
separation we find that it can also be calculated semi-analytically for any
initial state for the forwards quench but the reverse quench leads to non-
physical divergences, as is the case for two-body quench dynamics Kerin and
Martin (2020).

##  II Overview of the Three-Body Problem

Our starting point is the Hamiltonian of three non-interacting bodies in a
three-dimensional spherical harmonic trap:

$\displaystyle\hat{H}=\sum_{k=1}^{3}\left[\frac{-\hbar^{2}}{2m_{k}}\nabla_{k}^{2}+\frac{m_{k}\omega^{2}r_{k}^{2}}{2}\right],$
(1)

where $\vec{r}_{k}$ is the position of the $k^{\rm th}$ particle, $m_{k}$ is
its mass, and $\omega$ is the trapping frequency.

In this paper we consider the case of two identical fermions interacting with
a distinct third particle. We define particle one to be the impurity
$(m_{1}=m_{\rm i})$ and particles two and three to be identical
$(m=m_{2}=m_{3})$. For convenience we define the reduced mass $\mu=m_{\rm
i}m/(m_{\rm i}+m)$ and the lengthscale, $a_{\mu}=\sqrt{\hbar/\mu\omega}$.

For such a system the hyperspherical formulation Werner and Castin (2006a)
gives a closed form solution for the wavefunction. However, because the
interactions are enforced with the Bethe-Peierls boundary condition, the
wavefunction can only be fully specified in the non-interacting and strongly
interacting (unitary) regimes.

We define the hyperradius $R$ and hyperangle $\alpha$

$\displaystyle R^{2}=\sqrt{r^{2}+\rho^{2}},\quad\alpha=\arctan{(r/\rho)},$ (2)

where

$\displaystyle\vec{r}$ $\displaystyle=$
$\displaystyle\vec{r}_{2}-\vec{r}_{1},$ (3) $\displaystyle\vec{\rho}$
$\displaystyle=$ $\displaystyle\frac{1}{\gamma}\left(\vec{r}_{3}-\frac{m_{\rm
i}\vec{r}_{1}+m\vec{r}_{2}}{m_{\rm i}+m}\right),$ (4) $\displaystyle\gamma$
$\displaystyle=$ $\displaystyle\frac{\sqrt{m_{\rm i}(m_{\rm i}+2m)}}{m+m_{\rm
i}},$ (5)

and the centre-of-mass (COM) coordinate is

$\displaystyle\vec{C}$ $\displaystyle=$ $\displaystyle\frac{m_{\rm
i}\vec{r}_{1}+m\vec{r}_{2}+m\vec{r}_{3}}{m_{\rm i}+2m}.$ (6)

The centre-of-mass Hamiltonian is a simple harmonic oscillator (SHO)
Hamiltonian of a single particle of mass $M=2m+m_{\rm i}$ and position
$\vec{C}$. As such the centre-of-mass wavefunction is a SHO wavefunction. The
relative Hamiltonian is given

$\displaystyle\hat{H}_{\rm rel}$ $\displaystyle=$
$\displaystyle\frac{-\hbar^{2}}{2\mu}\Bigg{(}\frac{\partial^{2}}{\partial
R}+\frac{1}{R^{2}\sin(\alpha)\cos(\alpha)}\frac{\partial^{2}}{\partial\alpha^{2}}\left(\cos(\alpha)\sin(\alpha)\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\right)$
$\displaystyle+$ $\displaystyle\frac{5}{R}\frac{\partial}{\partial
R}-\frac{4}{R^{2}}-\frac{\hat{\Lambda}_{r}^{2}}{R^{2}\sin(\alpha)}-\frac{\hat{\Lambda}_{\rho}^{2}}{R^{2}\cos(\alpha)}\Bigg{)}+\frac{\mu\omega^{2}R^{2}}{2},$

where $\hat{\Lambda}_{r}^{2}$ and $\hat{\Lambda}_{\rho}^{2}$ are the angular
momentum operators in the $\hat{r}$ and $\hat{\rho}$ coordinate systems.

We define the trial wavefunction Werner (2008)

$\displaystyle\psi_{\rm 3b}^{\rm rel}$ $\displaystyle=$ $\displaystyle
N_{qls}\frac{F_{qs}(R)}{R^{2}}(1-\hat{P}_{23})\frac{\varphi_{ls}(\alpha)}{\sin(2\alpha)}Y_{lm}(\hat{\rho}),$
(8)

where $N_{qls}$ is the normalisation constant, $F_{qs}(R)$ is the hyperradial
wavefunction,
$\phi_{ls}(\alpha)=(1-\hat{P}_{23})\varphi_{ls}(\alpha)Y_{lm}(\hat{\rho})/\sin(2\alpha)$
is the hyperangular wavefunction, and $\hat{P}_{23}$ is the particle exchange
operator which swaps the positions of particles two and three.

Three conditions determine the functional forms of $F_{qs}(R)$ and
$\varphi_{ls}(\alpha)$,

$\displaystyle\varphi_{ls}\left(\frac{\pi}{2}\right)$ $\displaystyle=$
$\displaystyle 0,$ (9) $\displaystyle s^{2}\varphi_{ls}(\alpha)$
$\displaystyle=$
$\displaystyle-\varphi_{ls}^{\prime\prime}(\alpha)+\frac{l(l+1)}{\cos^{2}(\alpha)}\varphi_{ls}(\alpha),$
(10) $\displaystyle E_{\rm rel}$ $\displaystyle=$
$\displaystyle\frac{-\hbar^{2}}{4\mu}\left(F^{\prime\prime}(R)+\frac{F^{\prime}(R)}{R}\right)$
(11) $\displaystyle+$ $\displaystyle\left(\frac{\hbar^{2}s^{2}}{4\mu
R^{2}}+\mu\omega^{2}R^{2}\right)F(R).$

The first is enforced because a divergence at $\alpha=\pi/2$ is non-physical,
the second and third come from the Schrödinger equation. $l\in\mathbb{Z}_{\geq
0}$ is the angular momentum quantum number, $q\in\mathbb{Z}_{\geq 0}$ and
$s\in\mathbb{R}_{>0}$ are the energy eigenvalues and the energy is given
$E_{\rm rel}=(2q+l+s+1)\hbar\omega$. $F_{qs}(R)$ and $\varphi_{ls}(\alpha)$
are given Werner and Castin (2006a, b); Liu  _et al._ (2010)

$\displaystyle
F_{qs}(R)=\tilde{R}^{s}e^{-\tilde{R}^{2}/2}L_{q}^{s}\left(\tilde{R}^{2}\right),$
(12) $\displaystyle\varphi_{ls}(\alpha)=\cos^{l+1}(\alpha)$
$\displaystyle\qquad\qquad\times{}_{2}F_{1}\left(\frac{l+1-s}{2},\frac{l+1+s}{2};l+\frac{3}{2};\cos^{2}(\alpha)\right),$
(13)

where $\tilde{R}=R/a_{\mu}$.

The contact interactions are enforced by the Bethe-Peierls condition Bethe and
Peierls (1935)

$\displaystyle\lim_{r_{ij}\rightarrow
0}\left[\frac{d(r_{ij}\Psi)}{dr_{ij}}\frac{1}{r_{ij}\Psi}\right]=\frac{-1}{a_{\rm
s}},$ (14)

where $\Psi$ is the total three-body wavefunction,
$r_{ij}=|\vec{r}_{i}-\vec{r}_{j}|$, and $a_{\rm s}$ is the s-wave scattering
length.

In the non-interacting limit ($a_{\rm s}\rightarrow 0$) Eq. (14) implies, for
all values of $\kappa=m/m_{\rm i}$,

$\displaystyle s=\begin{cases}2n+4&l=0\\\ 2n+l+2&l>0\end{cases},$ (15)

where $n\in\mathbb{Z}_{\geq 0}$. In the unitary limit ($a_{\rm
s}\rightarrow\infty$) the Bethe-Peierls boundary condition gives the
transcendental equation

$\displaystyle 0=$
$\displaystyle\frac{d\varphi_{ls}}{d\alpha}\bigg{|}_{\alpha=0}-(-1)^{l}\frac{(1+\kappa)^{2}}{\kappa\sqrt{1+2\kappa}}\varphi_{ls}\left(\arctan(\frac{\sqrt{1+2\kappa}}{\kappa})\right).$
(16)

The values of $s$ at unitarity for $l=0$ and a variety of $\kappa$ are given
in Table 1.

$l=0$ | $\kappa=0.1$ | $\kappa=1$ | $\kappa=10$
---|---|---|---
$n$ | $s_{nl}$
0 | 2.004… | 2.166… | 3.316…
1 | 5.817… | 5.127… | 4.707…
2 | 6.195… | 7.114… | 6.747…
3 | 9.685… | 8.832… | 8.876…
Table 1: Some three-body $s$-eigenvalues at unitarity with $l=0$ and
$\kappa=0.1,1,10$, to three decimal places.

From this it is possible, in the non-interacting and unitary regimes, to
evaluate the eigenenergies and eigenstates of the three-particle system. From
this foundation in the following section we utilise these states to evaluate
the quench dynamics of this system.

##  III Quench Dynamics

Below we investigate the behaviour of the system after a quench in the s-wave
scattering length, $a_{\rm s}$. We are interested in the forwards quench (non-
interacting to unitary) and the backwards quench (unitary to non-interacting).
Specifically, we calculate the Ramsey signal, $S(t)$, and the particle
separation, $\langle\tilde{R}(t)\rangle$. In order to do this we need to
calculate various integrals involving the wavefunction. First we have the
Jacobian

$\displaystyle
dV=d\vec{r}_{1}d\vec{r}_{2}d\vec{r}_{3}=\frac{R^{5}}{4}\sin^{2}(2\alpha)\gamma^{3}dRd\alpha
d\vec{\Omega}_{r}d\vec{\Omega}_{\rho}d\vec{C}.\quad$ (17)

We make the definitions

$\displaystyle\bra{F_{qs}(R)}\ket{F_{qs}(R)}=\int_{0}^{\infty}RF_{qs}(R)^{*}F_{qs}(R)dR,$
(18) $\displaystyle\bra{\phi_{ls}(\alpha)}\ket{\phi_{ls}(\alpha)}=$
$\displaystyle\qquad\int\int\int_{0}^{\dfrac{\pi}{2}}\phi_{ls}(\alpha)^{*}\phi_{ls}(\alpha)2\sin^{2}(2\alpha)d\alpha
d\vec{\Omega}_{r}d\vec{\Omega}_{\rho}.\quad$ (19)

In this work we wish to calculate time-varying observables of a post-quench
system. To do this we need the time-dependent post-quench wavefunction. The
COM wavefunction is independent of $a_{\rm s}$ and so is unaffected by the
quench. As such we only need the time-dependent post-quench relative
wavefunction, it is given

$\displaystyle\ket{\psi_{\rm 3b}^{\rm rel}(t)}$ $\displaystyle=$
$\displaystyle e^{-i\hat{H}_{\rm rel}t/\hbar}\ket{F_{q_{\rm i}s_{\rm
i}}\phi_{l_{\rm i}s_{\rm i}}},$ $\displaystyle=$
$\displaystyle\sum_{q,s}\bra{F_{qs}\phi_{l_{\rm i}s}}\ket{F_{q_{\rm i}s_{\rm
i}}\phi_{l_{\rm i}s_{\rm i}}}e^{-iE_{ql_{\rm
i}s}t/\hbar}\ket{F_{qs}\phi_{l_{\rm i}s}},$

where quantum numbers with an i subscript are the initial quantum numbers and
the summation is over all eigenstates of the post-quench system.

###  III.1 Ramsey signal

The Ramsey signal is defined as the wavefunction overlap of the initial and
final states, it is given

$\displaystyle S(t)$ $\displaystyle=$ $\displaystyle\bra{\Psi_{\rm
i}(t)}\ket{\Psi^{\prime}(t)},$ $\displaystyle S(t)$ $\displaystyle=$
$\displaystyle\sum_{j}|\bra{\Psi_{\rm
i}(0)}\ket{\Psi^{\prime}_{j}}|^{2}e^{-i(E_{j}-E_{\rm i})t/\hbar},$ (21)

where $\Psi_{\rm i}$ is the initial state with energy $E_{\rm i}$,
$\Psi^{\prime}$ is the final state and the $\Psi^{\prime}_{j}$s are the
eigenstates of the post-quench system with energy $E_{j}$, and the sum over
$j$ is a sum over all post-quench states.

Since the COM wavefunction is unaffected by the quench the Ramsey signal is
given

$\displaystyle S(t)$ $\displaystyle=$ $\displaystyle\sum_{q,s}|\bra{F_{q_{\rm
i}s_{\rm i}}\phi_{l_{\rm i}s_{\rm i}}}\ket{F_{qs}\phi_{l_{\rm
i}s}}|^{2}e^{-i(E_{ql_{\rm i}s}-E_{q_{\rm i}l_{\rm i}s_{\rm
i}})t/\hbar},\quad$ (22)

where indices with subscript i are the eigenvalues of the initial state and
the unlabelled indices correspond to the post-quench eigenvalues. Note that
$\bra{\phi_{ls}}\ket{\phi_{l^{\prime}\neq ls^{\prime}}}=0$, i.e. hyperangular
states of different angular momenta are orthogonal.

To evaluate the Ramsey signal we need to evaluate the hyperradial integral
$\bra{F_{qs}}\ket{F_{q^{\prime}s^{\prime}}}$ and the hyperangular integral
$\bra{\phi_{ls}}\ket{\phi_{l^{\prime}s^{\prime}}}$.

The hyperradial integral is given Srivastava  _et al._ (2003)

$\displaystyle\bra{F_{qs}(R)}\ket{F_{q^{\prime}s^{\prime}}(R)}=$
$\displaystyle\qquad\frac{a_{\mu}^{2}}{2}\binom{q+s^{*}}{q}\binom{q^{\prime}+\frac{s^{\prime}-s^{*}}{2}-1}{q^{\prime}}\Gamma\left(\frac{s^{*}+s^{\prime}+2}{2}\right){}_{3}F_{2}\bigg{(}-q,\frac{s^{*}+s^{\prime}+2}{2},\frac{s^{*}-s^{\prime}+2}{2};s^{*}+1,\frac{s^{*}-s^{\prime}+2}{2}-q^{\prime};1\bigg{)}.\qquad$
(23)

For $s=s^{\prime}$ and $q\neq q^{\prime}$ the hyperradial integral vanishes.
In order to evaluate the hyperangular integral we need to perform a coordinate
rotation on the permuted part of the wavefunction to transform it into the
same Jacobi coordinate system as the unpermuted part Braaten and Hammer
(2006); Thøgersen (2009); Fedorov and Jensen (2001, 1993). However this
restricts us to the $l=0$ case due to the spherical harmonic term. For the
general mass case the overlap of $\phi_{0s}$ and $\phi_{0s^{\prime}}$ is
given,

$\displaystyle\bra{\phi_{0s}(\alpha)}\ket{\phi_{0s^{\prime}}(\alpha)}=8\pi\int_{0}^{\pi/2}\left((1+\hat{Q})\frac{\varphi_{0s}(\alpha)}{\sin(2\alpha)}\right)^{*}\left((1+\hat{Q})\frac{\varphi_{0s^{\prime}}(\alpha)}{\sin(2\alpha)}\right)\sin^{2}(2\alpha)d\alpha,$
(24)
$\displaystyle=16\pi\Bigg{[}\int_{0}^{\pi/2}\varphi_{0s}^{*}(\alpha)\varphi_{0s^{\prime}}(\alpha)d\alpha-\frac{(1+\kappa)^{2}}{2\kappa\sqrt{1+2\kappa}}\int_{0}^{\pi/2}\varphi_{0s}^{*}(\alpha)\Bigg{[}\int_{|\arctan(\sqrt{\frac{1}{\kappa}(2+\frac{1}{\kappa})})-\alpha|}^{\pi/2-|\pi/2-\arctan(\sqrt{\frac{1}{\kappa}(2+\frac{1}{\kappa})})-\alpha|}\varphi_{0s^{\prime}}(\alpha^{\prime})d\alpha^{\prime}\Bigg{]}d\alpha\Bigg{]}.\>\>$
(25)

For $l=0$ Eq. (13) reduces to Werner (2008); Fedorov and Jensen (2001)

$\displaystyle\varphi_{0s}(\alpha)\propto\sin(s\left(\frac{\pi}{2}-\alpha\right)).\quad$
(26)

Note that combining Eqs. (25) and (26) for $\kappa=1$ does not give the same
result as in Refs. Colussi (2019); Werner (2008) and only agrees when
$s^{\prime}=4$ and $s$ is one of the solutions to Eq. (16) for $\kappa=1$.
This is because the latter references appropriately substitute Eq. (16) for
$\kappa=1$ into the result for the hyperangular integral to simplify the
expression.

However, for a general value of $\kappa$ we have to be more careful. For a
specific $\kappa$ the $s$ eigenvalues produced by Eq. (16) are orthogonal to
one another when Eq. (25) is evaluated only for that specific value of
$\kappa$, i.e. the $s$-eigenspectrum for $\kappa=x$ produces an orthogonal set
of states only when Eq. (25) is evaluated using $\kappa=x$. However, the non-
interacting values of $s$ ($s=4,6,8,\dots$) are orthogonal to one another
whatever the value of $\kappa$.

Now that we can calculate the wavefunction overlaps we can evaluate the Ramsey
signal for the forwards and backwards quenches for any initial state with
$l=0$ and any mass imbalance. In Figs. 1 and 2 we plot the forwards and
backwards quenches respectively for a variety of $\kappa$ and initial states.

Before we discuss the properties of the Ramsey signals, in Figs. 1 and 2, it
is worth exploring what we may expect in a general sense. Assuming in Eq. (22)
that the sum is dominated by a few terms, the general form for $S(t)$ can be
represented by $S(t)\approx Ae^{-at}+Be^{-bt}+Ce^{-ct}$. The magnitude of the
signal oscillates with angular frequencies $(a-b),(b-c),(a-c)$ with the more
heavily weighted terms being more significant. However the phase of the signal
is dominated by the angular frequencies of the most heavily weighted terms not
the difference between angular frequencies of different terms. In our case the
angular frequencies ($a,b,c,\dots$) are the differences between a post-quench
eigenenergy and the initial energy, and the difference between the angular
frequencies ($(a-b),(b-c),(a-c),\dots$) is the difference between different
post-quench eigenenergies.

Figure 1: Ramsey signal for the forwards quench for a variety initial states
and mass imbalances. In all panels the dot-dashed red line corresponds to
$(q_{\rm i},s_{\rm i})=(0,4)$, the dashed green line corresponds to $(q_{\rm
i},s_{\rm i})=(1,4)$, and the solid blue line corresponds to $(q_{\rm
i},s_{\rm i})=(2,4)$. The upper, middle, and lower plots correspond to
$\kappa=0.1,1,10$ respectively. Each Ramsey signal is obtained by evaluating
Eq. (22) with 10 terms in each sum, i.e. 100 terms total. In this limit we
find that the sum is convergent.

In Fig. 1 we plot the Ramsey signals of the forwards quench for a variety of
initial states and mass imbalances. We find an irregularly repeating Ramsey
signal, this is in contrast to similar calculations performed for the two-body
problem Kerin and Martin (2020) where the signal is periodic with period
$\pi/\omega$. This irregular periodicity arises from the unitary
$s$-eigenspectrum. In this case angular frequencies of each term in Eq. (22)
are irrational and so are the differences between them, in general. These
irrational angular frequencies lead to the irregular phase and magnitude.

For the $\kappa=0.1$ forwards quench the most heavily weighted terms in Eq.
(22) for the $(q_{\rm i},s_{\rm i})=(0,4)$ initial state (red line in the
upper panel of Fig. 1) are the
$|\bra{F_{0,4}\phi_{0,4}}\ket{F_{0,2.004}\phi_{0,2.004}}|^{2}\approx 0.539$,
$|\bra{F_{0,4}\phi_{0,4}}\ket{F_{0,5.817}\phi_{0,5.817}}|^{2}\approx 0.134$
and $|\bra{F_{0,4}\phi_{0,4}}\ket{F_{1,2.004}\phi_{0,2.004}}|^{2}\approx
0.179$ terms. The corresponding angular frequencies are $\approx-2\omega$,
$\approx 1.8\omega$, $\approx 0.004\omega$ respectively. Hence the phase
displays features with a period of $\approx\pi/\omega$, and the magnitude has
three main modes with periods $\approx 0.5\pi/\omega$, $\pi/\omega$ and
$\approx 1.1\pi/\omega$.

For the $\kappa=1$ forwards quench the most significant terms for the $(q_{\rm
i},s_{\rm i})=(0,4)$ initial state are
$|\bra{F_{0,4}\phi_{0,4}}\ket{F_{0,2.166}\phi_{0,2.166}}|^{2}\approx 0.499$,
$|\bra{F_{0,4}\phi_{0,4}}\ket{F_{0,5.127}\phi_{0,5.127}}|^{2}\approx 0.297$,
and $|\bra{F_{0,4}\phi_{0,4}}\ket{F_{1,2.166}\phi_{1,2.166}}|^{2}\approx
0.13$. The corresponding angular frequencies are $\approx-1.8\omega$, $\approx
1.1\omega$ and $\approx 0.17\omega$ respectively. The phase is then dominated
by a mode with period $\approx 1.1\pi/\omega$ and the magnitude is dominated
by modes with periods $\approx 2\pi/3\omega,\pi/\omega$ and $\approx
2\pi/\omega$.

For the $\kappa=10$ forwards quench the most significant terms for the
$(q_{\rm i},s_{\rm i})=(0,4)$ initial state are
$|\bra{F_{0,4}\phi_{0,4}}\ket{F_{0,3.3169}\phi_{0,3.3169}}|^{2}\approx 0.607$,
and $|\bra{F_{0,4}\phi_{0,4}}\ket{F_{0,4.707}\phi_{0,4.707}}|^{2}\approx
0.350$. The corresponding angular frequencies are $\approx-0.7\omega/\pi$ and
$\approx 0.7\omega/\pi$ respectively hence the phase has an approximate
periodicity of $2.8\pi/\omega$ and the magnitude has a period of $\approx
1.4\pi/\omega$.

Figure 2: Ramsey signal for the backwards quench for a variety of initial
states and mass imbalances. In all panels the dot-dashed red line corresponds
to $q_{\rm i}=0$, the dashed green line corresponds to $q_{\rm i}=1$, and the
solid blue line corresponds to $q_{\rm i}=2$. The upper, middle and lower
plotss correspond to $\kappa=0.1,1,10$ and $s_{\rm
i}=2.004\dots,2.166\dots,3.316\dots$ respectively. Each Ramsey signal is
obtained by evaluating Eq. (22) with 10 terms in each sum, i.e. 100 terms
total. In this limit we find that the sum is convergent.

In Fig. 2 we plot the Ramsey signals of the backwards quench for a variety of
initial states and mass imbalances. We find a strongly periodic Ramsey signal.
In the non-interacting regime we have $s=4,6,8,\dots$ for every value of
$\kappa$. This means that every oscillatory term in the magnitude has an
angular frequency that is a multiple of two, hence the period of $\pi/\omega$.
However the behaviour of the phase is still dominated by the largest terms in
the summation. For $\kappa=0.1$ and $(q_{\rm i}=0,s_{\rm i}=2.00449\dots)$ the
largest term is
$|\bra{F_{0,2.004}\phi_{0,2.004}}\ket{F_{0,4}\phi_{0,4}}|^{2}\approx 0.539$
and this defines the period of $\approx\pi/\omega$ we see in the phase.
Similarly for $\kappa=1$ and $(q_{\rm i}=0,s_{\rm i}=2.166\dots)$ the largest
term is $|\bra{F_{0,2.166}\phi_{0,2.166}}\ket{F_{0,4}\phi_{0,4}}|^{2}\approx
0.499$ defining a period of $\approx 1.1\pi/\omega$ for the phase, and for
$\kappa=10$ and $(q_{\rm i}=0,s_{\rm i}=3.3169\dots)$ the largest term is
$|\bra{F_{0,3.3169}\phi_{0,3.3169}}\ket{F_{0,4}\phi_{0,4}}|^{2}\approx 0.607$
defining a period of $\approx 2.8\pi/\omega$ for the phase.

###  III.2 Particle separation

The Ramsey signal is not the only quench observable we investigate. We can
also calculate the expectation value of $\tilde{R}$, i.e. the particle
separation.

The expectation value of $\tilde{R}(t)$ is given

$\displaystyle\langle\tilde{R}(t)\rangle=\bra{\Psi_{\rm
i}(t)}\tilde{R}\ket{\Psi(t)}$ $\displaystyle=\sum_{j,j^{\prime}}\bra{\Psi_{\rm
i}(0)}\ket{\Psi^{\prime}_{j}}\bra{\Psi^{\prime}_{j^{\prime}}}\ket{\Psi_{\rm
i}(0)}\bra{\Psi^{\prime}_{j}}\tilde{R}\ket{\Psi^{\prime}_{j^{\prime}}}e^{-i(E_{j}-E_{j^{\prime}})t/\hbar},$

where $\Psi_{\rm i}$ is the initial pre-quench state with energy $E_{\rm i}$
and $\Psi$ is the post-quench state. $\Psi^{\prime}_{j}$ and
$\Psi^{\prime}_{j^{\prime}}$ are eigenstates of the post-quench system with
eigenenergy $E_{j}$ and $E_{j^{\prime}}$ respectively, with the sum over $j$
and $j^{\prime}$ taken over all post-quench eigenstates.

The COM wavefunction is independent of the interparticle interaction and does
not change when the system is quenched and does not impact the post-quench
dynamics. Due to the hyperangular wavefunction’s orthogonality in $s$, two
sums over $s$ and $s^{\prime}$ collapse into a single sum over $s$. Hence
$\langle\tilde{R}(t)\rangle$ is given

$\displaystyle\langle\tilde{R}(t)\rangle$ $\displaystyle=$
$\displaystyle\sum_{q^{\prime},q}\sum_{s}\bra{F_{q_{\rm i}s_{\rm
i}}\phi_{l_{\rm i}s_{\rm i}}}\ket{F_{q^{\prime}s}\phi_{l_{\rm
i}s}}\bra{F_{qs}\phi_{l_{\rm i}s}}\ket{F_{q_{\rm i}s_{\rm i}}\phi_{l_{\rm
i}s_{\rm i}}}$ (28) $\displaystyle\times\bra{F_{q^{\prime}s}\phi_{l_{\rm
i}s}}\tilde{R}\ket{F_{qs}\phi_{l_{\rm i}s}}e^{-i(E_{ql_{\rm
i}s}-E_{q^{\prime}l_{\rm i}s})t/\hbar}.$

As in Eq. (22) indices with subscript i refer to the initial state and indices
with no subscript refer to the post-quench eigenstates. In Sec. III.1 all
integrals required to evaluate Eq. (28) are presented with the exception of
$\bra{F_{q^{\prime}s}}\tilde{R}\ket{F_{qs}}$.

All integrals required for calculating Eq. (28) are calculated in Sec. III.1
except for $\bra{F_{q^{\prime}s}}\tilde{R}\ket{F_{q^{\prime}s}}$. This is
given Srivastava  _et al._ (2003)

$\displaystyle\bra{F_{qs}}\tilde{R}\ket{F_{q^{\prime}s^{\prime}}}=$
$\displaystyle\frac{a_{\mu}^{2}}{2}\binom{q+s}{q}\binom{q^{\prime}+\frac{s^{\prime}}{2}-\frac{s}{2}-\frac{3}{2}}{q^{\prime}}\Gamma\left(\frac{s+s^{\prime}+3}{2}\right)$
$\displaystyle\times{}_{3}F_{2}\left(-q,\frac{s^{\prime}+s+3}{2},\frac{s-s^{\prime}+3}{2};s+1,\frac{s-s^{\prime}+3}{2}-q^{\prime};1\right).$

This result combined with the results of Sec. III.1 allows us to calculate
$\langle\tilde{R}(t)\rangle$ for the forwards and backwards quench for any
initial state with $l=0$ and any mass imbalance.

Figure 3: $\langle\tilde{R}(t)\rangle$ for the forwards quench from the non-
interacting ground state, $(q_{\rm i},s_{\rm i})=(0,4)$, for $\kappa=0.1$
(upper panel), $\kappa=1$ (middle panel) and $\kappa=10$ (lower panel). The
dot-dashed red line corresponds to $N_{\rm max}=3$, the dashed green line to
$N_{\rm max}=6$, the solid blue to $N_{\rm max}=12$ and the dotted black line
to $N_{\rm max}=24$. As can be seen from the plots the sum is convergent with
$N_{\rm max}$.

In Sec. III.1 we note that the magnitude and phase of the Ramsey signal for
the forwards quench is irregular and this is because the $s$ eigenvalues are
irrational at unitarity. Additionally the phase for the backwards quench is
also irregular and this is due to the irrationality of $s_{\rm i}$. However
the angular frequencies of the terms in Eq. (28) are always even integers
because the $s$ contributions to the energies cancel leaving angular
frequencies proportional to $(2q^{\prime}-2q)\omega$ and
$q,q^{\prime}\in\mathbb{Z}_{\geq 0}$. This results in the angular frequency of
every term in the summation being a multiple of $2\omega$ causing
$\langle\tilde{R}(t)\rangle$ to have a period of $\pi/\omega$ in both the
forwards and backwards quench.

In Fig. 3 we have plotted $\langle\tilde{R}(t)\rangle$ for the forwards quench
for $\kappa=0.1$, $\kappa=1$ and $\kappa=10$ in the upper, middle and lower
panels respectively. For each $\kappa$ we have taken the initial state to be
the ground state and calculated the sum in Eq. (28) over $q$, $q^{\prime}$ and
$s$ up to $N_{\rm max}=3,6,12$ and $24$, where $N_{\rm max}$ is the number of
terms in each individual sum, so there are $N_{\rm max}^{3}$ terms total. In
each case we find that the results for $\langle\tilde{R}(t)\rangle$ have
converged at $N_{\rm max}=24$. As discussed above $\langle\tilde{R}(t)\rangle$
oscillates with a period $\pi/\omega$. Additionally we observe that as
$\kappa$ increases the amplitude of the oscillations decreases. As
$\kappa\rightarrow\infty$ the unitary s-eigenspectrum approaches $s\rightarrow
4,6,8\dots$ Kerin and Martin (2022) and we expect the initial non-interacting
ground state to overlap perfectly with the unitary ground state. As such we
expect the oscillations to decrease as $\kappa\rightarrow\infty$.

Figure 4: $\langle\tilde{R}(t)$ for the backwards quench from the interacting
ground state for a variety of $\kappa$. In all panels $q_{\rm i}=0$, in the
upper panel ($\kappa=0.1$) $s_{\rm i}=2.004\dots$, in the middle panel
($\kappa=1$) $s_{\rm i}=2.166\dots$, and in the lower panel ($\kappa=10$)
$s_{\rm i}=3.3169$. The dot-dashed red line corresponds to $N_{\rm max}=3$,
the dashed green line to $N_{\rm max}=6$, the solid blue to $N_{\rm max}=12$
and the dotted black line to $N_{\rm max}=24$. This sum diverges
logarithmically.

In Fig. 4 we have considered $\langle\tilde{R}(t)\rangle$ for the backwards
quench for the $\kappa=0.1$, $\kappa=1$ and $\kappa=10$ cases in the upper,
middle and lower panels respectively. For each $\kappa$ we have taken the
initial state to be the ground state and calculated the sum in Eq. (28) over
$q$, $q^{\prime}$ and $s$ up to $N_{\rm max}=3,6,12$ and $24$. Similar to
$r=|\vec{r}_{1}-\vec{r}_{2}|$ in the two body case Kerin and Martin (2020) we
observe in the backwards quench that the particle separation diverges away
from $\omega t=n\pi$. For the three-body case we find that
$\langle\tilde{R}(t=(n+1/2)\pi/\omega)\rangle$ diverges logarithmically with
$N_{\rm max}$. More specifically we find that if the number of terms in the
sums over $q$ and $q^{\prime}$ is fixed then as more terms are included in the
sum over $s$ $\langle\tilde{R}(t)\rangle$ converges. The divergence in
$\langle\tilde{R}(t)\rangle$ comes from the sums over $q$ and $q^{\prime}$ and
therefore from the hyperradial wavefunction.

In the two-body case it has been suggested that the divergence in $\langle
r\rangle=\langle|\vec{r}_{1}-\vec{r}_{2}|\rangle$ is due to the $1/r$
divergence in the two-body relative wavefunctionKerin and Martin (2020). In
the three-body case the hyperradial part of Eq. (8) does not have a $1/R$
divergence and has a cusp at $\tilde{R}\rightarrow 0$ for the unitary
$\kappa=0.1$ and $\kappa=1$ ground states but not for the $\kappa=10$ ground
state. Nonetheless the divergence in $\langle\tilde{R}(t)\rangle$ is present
in all three cases. More specifically we find that for $s_{\rm i}<3$ the
initial wavefunction exhibits a cusp, whilst for $s_{\rm i}\geq 3$ there is no
cusp in the initial wavefunction. Regardless of which regime the initial state
is in the logarithmic divergence of $\langle\tilde{R}(t)\rangle$ persists.

However, it is clear that the finite range of the interaction in a real system
provides a natural cut-off in the sum in Eq. (28). A finite range of
interaction defines a minimum de Broglie wavelength which in turn defines a
maximum energy which defines a maximum number of terms in the sums of Eq. (28)
and thus a maximum $\langle\tilde{R}(t)\rangle$. For a system of three sodium
atoms (i.e. $\kappa=1$) in a 1kHz trap and an interaction cut-off of
$10^{-9}$m we obtain a cut off energy of $E_{\rm rel}\approx 8.7\times
10^{6}\hbar\omega$ and thus expect a maximum $\langle\tilde{R}(t)\rangle$ of
$\approx 11$. This is a significantly larger oscillation than in the forwards
quench.

##  IV Conclusion

In this work we examine the time evolution of quenched systems. We consider a
harmonically trapped system of two identical fermions plus a distinct particle
interacting via a contact interaction where the system is quenched from non-
interacting to strongly interacting or vice versa. We calculate the static
wavefunction in both the non-interacting and strongly interacting regimes for
general mass and use these solutions to calculate two time-dependent post-
quench observables; the Ramsey signal, $S(t)$, and the particle separation,
$\langle\tilde{R}(t)\rangle$. These observables are calculated for both the
forwards and backwards quenches.

For the Ramsey signal we find an irregularly repeating signal for the forwards
quench. This irregularity is due to the irrationality of the unitary energy
spectrum and this irregular behaviour is more pronounced for small $\kappa$.
For the reverse quench the magnitude of the Ramsey signal oscillates with
period $\pi/\omega$ as the non-interacting energies are integer multiples and
the phase has an irregular period due the the irrational initial energy.

For the particle separation the forwards quench yields the expected
oscillating result, however the period is $\pi/\omega$ as the irrational
contributions to the unitary eigenenergies cancel. For the backwards quench we
find that the particle separation diverges logarithmically, similar to
divergence in $r=|\vec{r}_{2}-\vec{r}_{1}|$ for the same quench performed on
the two-body system Kerin and Martin (2020). By enforcing a cut-off based on
the van der Waals interaction range we estimate a very large amplitude of
oscillations $\approx 5a_{\mu}$ for $\kappa=1$. From a physical perspective it
is unclear why the divergence occurs for the backwards but not forwards
quench.

Finally we note that experimental testing of these theoretical predictions is
within reach. Few-atom systems can be reliably constructed with modern
techniques Serwane  _et al._ (2011); Murmann  _et al._ (2015); Zürn  _et al._
(2013, 2012); Stöferle  _et al._ (2006), the quench in s-wave scattering
length can be achieved using Feshbach resonance Fano (1935); Feshbach (1958);
Tiesinga  _et al._ (1993); Chin  _et al._ (2010), and the Ramsey signal can be
measured using Ramsey interferometry techniques Cetina  _et al._ (2016).
Notable Ref. Guan  _et al._ (2019) measured the particle separation of two
harmonically trapped 6Li atoms following a quench in trap geometry rather than
s-wave scattering length.

##  V Acknowledgements

A.D.K. is supported by an Australian Government Research Training Program
Scholarship and by the University of Melbourne.

With thanks to Victor Colussi for illuminating discussions regarding the
evaluation of the hyperangular integral.

## References

  * Serwane  _et al._ (2011) F. Serwane, G. Zürn, T. Lompe, T. Ottenstein, A. Wenz, and S. Jochim, Science 332, 336 (2011).
  * Murmann  _et al._ (2015) S. Murmann, A. Bergschneider, V. M. Klinkhamer, G. Zürn, T. Lompe, and S. Jochim, Physical Review Letters 114, 080402 (2015).
  * Zürn  _et al._ (2013) G. Zürn, A. N. Wenz, S. Murmann, A. Bergschneider, T. Lompe, and S. Jochim, Physical Review Letters 111, 175302 (2013).
  * Zürn  _et al._ (2012) G. Zürn, F. Serwane, T. Lompe, A. N. Wenz, M. G. Ries, J. E. Bohn, and S. Jochim, Physical Review Letters 108, 075303 (2012).
  * Stöferle  _et al._ (2006) T. Stöferle, H. Moritz, K. Günter, M. Köhl, and T. Esslinger, Physical Review Letters 96, 030401 (2006).
  * Cui (2012) X. Cui, Few Body Systems 52, 65 (2012).
  * D’Incao  _et al._ (2018) J. P. D’Incao, J. Wang, and V. E. Colussi, Physical Review Letters 121, 023401 (2018).
  * Jonsell  _et al._ (2002) S. Jonsell, H. Heiselberg, and C. Pethick, Physical Review Letters 89, 250401 (2002).
  * Blume and Daily (2010) D. Blume and K. Daily, Physical review letters 105, 170403 (2010).
  * Kerin and Martin (2022) A. Kerin and A. Martin, arXiv preprint arXiv:2204.09205 (2022).
  * Kestner and Duan (2007) J. Kestner and L.-M. Duan, Physical Review A 76, 033611 (2007).
  * Liu  _et al._ (2009) X.-J. Liu, H. Hu, and P. D. Drummond, Physical Review Letters 102, 160401 (2009).
  * Liu  _et al._ (2010) X.-J. Liu, H. Hu, and P. D. Drummond, Physical Review A 82, 023619 (2010).
  * Rakshit  _et al._ (2012) D. Rakshit, K. M. Daily, and D. Blume, Physical Review A 85, 033634 (2012).
  * Kaplan and Sun (2011) D. B. Kaplan and S. Sun, Physical Review Letters 107, 030601 (2011).
  * Mulkerin  _et al._ (2012a) B. Mulkerin, C. Bradly, H. Quiney, and A. Martin, Physical Review A 85, 053636 (2012a).
  * Mulkerin  _et al._ (2012b) B. C. Mulkerin, C. J. Bradly, H. M. Quiney, and A. M. Martin, Physical Review A 86, 053631 (2012b).
  * Nascimbène  _et al._ (2010) S. Nascimbène, N. Navon, F. Jiang, K. Chevy, and C. Salomon, Nature 463, 1057 (2010).
  * Ku  _et al._ (2012) M. Ku, A. Sommer, L. Cheuk, and M. Zwierlein, Science 335, 563 (2012).
  * Levinsen  _et al._ (2017) J. Levinsen, P. Massignan, S. Endo, and M. M. Parish, Journal of Physics B: Atomic, Molecular and Optical Physics 50, 072001 (2017).
  * Daily and Blume (2010) K. Daily and D. Blume, Physical Review A 81, 053615 (2010).
  * Bougas  _et al._ (2022) G. Bougas, S. Mistakidis, P. Giannakeas, and P. Schmelcher, arXiv preprint arXiv:2205.02015 (2022).
  * Kerin and Martin (2020) A. Kerin and A. Martin, Physical Review A 102, 023311 (2020).
  * Werner and Castin (2006a) F. Werner and Y. Castin, Physical Review Letters 97, 150401 (2006a).
  * Werner (2008) F. Werner, Theses, Université Pierre et Marie Curieâ Paris VI, Paris France (2008).
  * Werner and Castin (2006b) F. Werner and Y. Castin, Physical Review A 74, 053604 (2006b).
  * Bethe and Peierls (1935) H. Bethe and R. Peierls, Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 148, 146 (1935).
  * Srivastava  _et al._ (2003) H. Srivastava, H. A. Mavromatis, and R. S. Alassar, Applied Mathematics Letters 16, 1131 (2003).
  * Braaten and Hammer (2006) E. Braaten and H.-W. Hammer, Physics Reports 428, 259 (2006).
  * Thøgersen (2009) M. Thøgersen, arXiv preprint arXiv:0908.0852 (2009).
  * Fedorov and Jensen (2001) D. V. Fedorov and A. Jensen, Journal of Physics A: Mathematical and General 34, 6003 (2001).
  * Fedorov and Jensen (1993) D. Fedorov and A. Jensen, Physical review letters 71, 4103 (1993).
  * Colussi (2019) V. E. Colussi, Atoms 7, 19 (2019).
  * Fano (1935) U. Fano, Il Nuovo Cimento 12, 154 (1935).
  * Feshbach (1958) H. Feshbach, Annals of Physics 5, 357 (1958).
  * Tiesinga  _et al._ (1993) E. Tiesinga, B. J. Verhaar, and H. T. C. Stoof, Physical Review A 47, 4114 (1993).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Reviews of Modern Physics 82, 1225 (2010).
  * Cetina  _et al._ (2016) M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. Walraven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt, M. Knap, _et al._ , Science 354, 96 (2016).
  * Guan  _et al._ (2019) Q. Guan, V. Klinkhamer, R. Klemt, J. H. Becher, A. Bergschneider, P. M. Preiss, S. Jochim, and D. Blume, Physical Review Letters 122, 083401 (2019).


# Do new investment strategies
take existing strategies’ returns
– An investigation into agent-based models – ††thanks: Note that the opinions
contained herein are solely those of the authors and do not necessarily
reflect those of SPARX Asset Management Co., Ltd.

Takanobu Mizuta  SPARX Asset Management Co. Ltd.,
Tokyo, Japan
https://orcid.org/0000-0003-4329-0645

###### Abstract

Commodity trading advisors (CTAs), who mainly trade commodity futures, showed
good returns in the 2000s. However, since the 2010’s, they have not performed
very well. One possible reason of this phenomenon is the emergence of short-
term reversal traders (STRTs) who prey on CTAs for profit. In this study, I
built an artificial market model by adding a CTA agent (CTAA) and STRT agent
(STRTA) to a prior model and investigated whether emerging STRTAs led to a
decrease in CTAA revenue to determine whether STRTs prey on CTAs for profit.
To the contrary, my results showed that a CTAA and STRTA are more likely to
trade and earn more when both exist. Therefore, it is possible that they have
a mutually beneficial relationship.

###### Index Terms:

CTA, Preying Investment Strategy, Agent-Based Model, Artificial Market Model

##  I Introduction

Commodity trading advisors (CTAs), who mainly trade commodity futures, showed
good returns in the 2000s. However, since the 2010’s, they have not performed
very well[1].

Clenow [1] reported the following points as reasons for this: (1) the yield of
government bonds decreased, (2) most CTAs utilized similar trend-following
strategies and tended to trade at the same time, leading to the increased
number of CTAs further impacting market prices, (3) the variation of investors
and investment strategies was diversified by alternative investments
spreading, and (4) short-term reversal traders (STRTs) who preyed on CTAs for
profit emerged.

In regard to Point (2), Nagao and Kobayashi [2] reported that CTAs utilized
similar and simple trend-following strategies to enable people to build an
imitated fund. On the basis of Point (4), Raschke and Connors [3] proposed a
short-term reversal strategy to prey on CTAs for profit.

However, the four points are not only inadequately founded but also contradict
each other. Point (2) implies that CTAs did not perform well due to the number
of CTAs increasing, but Point (3) implies that the number of CTAs has
relatively decreased. Point (4) implies that STRTs utilize reversal strategies
to trade on the opposite side and reduce market impacts by CTAs, who utilize
trend-following strategies. Therefore, Point (4) contradicts Point (2). When
one investment strategy fails to work, new strategies, especially those that
use new technologies or faster tradings, have tended to be blamed for a long
time. Since the second half of the 20th century, this tendency has become
stronger, and the media have often criticized new investment strategies using
new technologies, electronic trades, program trades, algorithm trades, high-
frequency trades, and AI traders. However, these critics seem to be missing
the point. Therefore, arguments that new strategies are used to impede
existing strategies, as mentioned in Point (4), should be discussed with
particular suspicion.

Empirical studies cannot be conducted to investigate situations that have
never occurred in actual financial markets, for example, comparing whether a
certain investment strategy exists. As so many factors affect price formation
in actual markets, empirical studies cannot be conducted to isolate the direct
effect of whether a certain investment strategy exists. In contrast, an
artificial market model, which is an agent-based model for a financial
market,111Excellent reviews include [4, 5, 6, 7, 8, 9]. can isolate the pure
contributions of an investment strategy. These are the advantages of
artificial market simulation studies. Articles in both Nature [10] and Science
[11] have argued that artificial market studies are expected to contribute to
a greater understanding of actual markets. Bookstaber [12] also argued that an
agent-based model can investigate financial crises better than classical
economics.

Many previous artificial market studies have contributed to explaining the
nature of financial market phenomena such as bubbles and crashes. Recent
artificial market studies have also contributed to discussions of appropriate
financial regulations and rules [7, 8, 9]. The JPX Working Paper series
includes various studies that contributed to such
discussions222https://www.jpx.co.jp/english/corporate/research-study/working-
paper/index.html.

In this study, I built an artificial market model by adding a CTA agent (CTAA)
and STRT agent (STRTA) to the prior model of Mizuta et al.[13] and
investigated whether emerging STRTAs led to a decrease in CTAA revenue to
determine whether STRTs prey on CTAs for profit, as stated in Point (4).

##  II Model

The model of Chiarella and Iori [14] is very simple but replicates long-term
statistical characteristics observed in actual financial markets: a fat tail
and volatility clustering. In contrast, that of Mizuta et al. [13] replicates
high-frequency micro structures, such as execution rates, cancel rates, and
one-tick volatility, that cannot be replicated with the model of Chiarella and
Iori [14]. Only fundamental and technical analysis strategies that exist
generally for any market at any time333Many empirical studies using
questionnaires found these strategies to be the majority, which are
comprehensively reviewed by Menkhoff and Taylor [15]. The empirical study
using market data by Yamamoto[16] showed that investors are switching
fundamental and technical analysis strategies. are implemented into the agent
model.

The simplicity of the model is very important for this study because
unnecessary replication of macro phenomena leads to models that are overfitted
and too complex. Such models prevent the understanding and discovery of
mechanisms affecting price formation because of the increase in related
factors. I explain the basic concept for constructing my artificial market
model in the review article [8, 9] and the Appendix “Basic Concept for
Constructing Model.”

In this study, the artificial market model was built by adding a CTAA and
STRTA to the prior model of Mizuta et al.[13].

This model contains one stock. The stock exchange uses a continuous double
auction to determine the market price. In the auction mechanism, multiple
buyers and sellers compete to buy and sell financial assets in the market, and
transactions can occur at any time whenever an offer to buy and one to sell
match. The minimum unit of price change is $\delta P$. The buy-order price is
rounded down to the nearest fraction, and the sell-order price is rounded up
to the nearest fraction.

The model includes $n$ normal agents (NAs), one CTAA, and one STRTA. The
agents always place an order for only one share. I implemented the variable
“tick time” $t$, which increases by one when an NA orders.

####  II-1 Normal agent (NA)

To replicate the nature of price formation in actual financial markets, I
introduced the NA to model a general investor. The number of NAs is $n$. NAs
can short sell freely. The holding positions are not limited, so NAs can take
an infinite number of shares for both long and short positions. First, at time
$t=1$, NA No. $1$ places an order to buy or sell its risk asset; then, at
$t=2,3,,,n$, NAs No. $2,3,,,n$ respectively place buy or sell orders. At
$t=n+1$, the model returns to the first NA and repeats this cycle. An NA
determines the order price and buys or sells using a combination of a
fundamental and technical analysis strategies to form an expectation of a risk
asset’s return. The expected return of agent $j$ for each risk asset at $t$ is

$r^{t}_{e,j}=(w_{1,j}\ln{\frac{P_{f}}{P^{t-1}}}+w_{2,j}\ln{\frac{P^{t-1}}{P^{t-\tau_{j}-1}}}+w_{3,j}\epsilon^{t}_{j})/\Sigma_{i}^{3}w_{i,j}$
(1)

where $w_{i,j}$ is the weight of term $i$ for agent $j$ and is independently
determined by random variables uniformly distributed on the interval
$(0,w_{i,max})$ at the start of the simulation for each agent. $\ln$ is the
natural logarithm. $P_{f}$ is a fundamental value and is a constant. $P^{t}$
is a mid-price (the average of the highest buy-order price and the lowest
sell-order price) at $t$, and $\epsilon^{t}_{j}$ is determined by random
variables from a normal distribution with average $0$ and variance
$\sigma_{\epsilon}$. $\tau_{j}$ is independently determined by random
variables uniformly distributed on the interval $(1,\tau_{max})$ at the start
of the simulation for each agent444When $t<\tau_{j}$, the second term of Eq.
(1) is zero..

The first term in Eq. (1) represents a fundamental strategy: the NA expects a
positive return when the market price is lower than the fundamental value, and
vice versa. The second term represents a technical analysis strategy using a
historical return: the NA expects a positive return when the historical market
return is positive, and vice versa. The third term represents noise.

After the expected return has been determined, the expected price is

$P^{t}_{e,j}=P^{t}\exp{(r^{t}_{e,j})}.$ (2)

Order prices are scattered around the expected price $P^{t}_{e,j}$ to
replicate many waiting limit orders. An order price $P^{t}_{o,j}$ is
determined by random variables uniformly distributed on the interval
$(P^{t}_{e,j}-P_{d},P^{t}_{e,j}+P_{d})$ where $P_{d}$ is a constant and means
dispersion of the order prices scattered. Whether the agent buys or sells is
determined by the magnitude relationship between $P^{t}_{e,j}$ and
$P^{t}_{o,j}$:

when $P^{t}_{e,j}>P^{t}_{o,j}$, the NA places an order to buy one share, and

when $P^{t}_{e,j}<P^{t}_{o,j}$, the NA places an order to sell one
share555When $t<t_{c}$, to generate enough waiting orders, the agent places an
order to buy one share when $P_{f}>P^{t}_{o,j}$, or to sell one share when
$P_{f}<P^{t}_{o,j}$. . The remaining order is canceled $t_{c}$ after the order
time.

####  II-2 CTA agent(CTAA) and STRT agent(STRTA)

One CTAA and STRTA exist. They are limited to hold only one share or minus one
share. The CTAA orders at $n\times\delta t$ ($n$ is a natural number and
$\delta t$ is constant), and the STRTA orders at $n\times\delta t+\delta t/2$.
They order only when market prices and their holdings meet the following
conditions666When $t<t_{c}$, they do not place an order to generate a
sufficient number of waiting orders.(see also Fig. 1). The tick time $t$ does
not increase whether they order.

Trades of the CTAA are as follows. If the CTAA has no share, it places a
market buy (sell) order777When an agent orders to buy (sell), if there is a
lower sell-order price (a higher buy-order price) than the agent’s order,
dealing immediately occurs. Such an order is called a “market order.” of one
share when $P^{t}$ reaches the maximum (minimum) within the last $\Delta
T_{1}$ (and then it will have one (minus one) share). If the CTAA has one
(minus one) share, it places a market sell (buy-back) order of one share when
$P^{t}$ reaches the minimum (maximum) within the last $\Delta T_{2}$ (and then
it will have no share).

Trades of the STRTA are as follows. If the STRTA has no share, it places a
market buy (sell) order of one share when $P^{t}$ reaches the minimum
(maximum) within the last $\Delta T_{2}$ (and then it will have one (minus
one) share). If the STRTA has one (minus one) share, it places market sell
(buy-back) order of one share when $P^{t}$ reaches the maximum (minimum)
within the last $\Delta T_{3}$ (and then it will have no share).

Figure 1: Order cycle of CTAA and STRTA Figure 2: Time evolution of prices when both CTAA and STRTA exist and not TABLE I: Returns of CTAA |  | STRTA
---|---
| exist | not exist
CTAA | exist | 105% | 90%
not exist | 93% | 79%

TABLE II: Trading volume of CTAA | STRTA
---|---
| exist | not exist
CTAA | exist | 892 | 865
not exist | 891 | 866

TABLE III: Returns of STRTA |  | STRTA
---|---
| exist | not exist
CTAA | exist | 293% | 205%
not exist | 265% | 185%

TABLE IV: Trading volume of STRTA | STRTA
---|---
| exist | not exist
CTAA | exist | 1739 | 1633
not exist | 1685 | 1587

TABLE V: Standard deviations of returns for $100$ tick times. |  | STRTA
---|---
| exist | not exist
CTAA | exist | 0.047% | 0.038%
not exist | 0.041% | 0.032%

TABLE VI: Standard deviations of returns for $20000$ tick times. | STRTA
---|---
| exist | not exist
CTAA | exist | 0.33% | 0.28%
not exist | 0.30% | 0.25%

Figure 3: Mechanism of co-prosperity

##  III Simulation results

In this study, I set the same parameters as in [13]888I explain how they
validated their model in the Appendix “Validation of the Model.”.
Specifically, I set
$n=1000,w_{1,max}=1,w_{2,max}=10,w_{3,max}=1,\tau_{max}=10000,\sigma_{\epsilon}=0.03,P_{d}=1000,t_{c}=10000,\delta
P=0.01$, and $P_{f}=10000$. The simulations ran to $t=1000000$. I simulated
four cases where the CTAA and STRTA existed or not, respectively, where the
other parameters were fixed and used the same random number table. I then
simulated these runs 100 times, changing the random number table each time,
and reported the results as the averaged statistical values.

Fig. 2 shows the time evolution of the prices when both the CTAA and STRTA
exist and when both do not exist. The two cases are not very different, and
therefore trades of them do not impact the prices very much.

Tables VI and VI show the returns and trading volume of the CTAA for the four
cases, respectively, where returns are defined as final earnings per the
fundamental price ($P_{f}$). The results for when the CTAA does not exist are
calculated on the basis that the CTAA would virtually trade with the best buy
or sell price of waiting orders, and those orders would not change. As the
CTAA does not exist, it cannot change the prices nor affect other agents, but
the figures can be estimated.

In Table VI, when the CTAA exists, it earns more when the STRTA exists than
when it does not even thought the time evolution of prices are not so
different as Fig. 2. Table VI shows that the CTAA can trade more if the STRTA
exists than without it. This means that the CTAA has more chances of trades
and can earn more when the STRTA exists.

Tables VI and VI show the returns and trading volume of the STRTA for the four
cases, respectively. The STRTA trades and earns more when the CTAA exists
similarly to the case of the CTAA in Tables VI and VI. In summary, the CTAA
and STRTA have more chances of trades and earning more when both exist.
Therefore, they have a mutually beneficial relationship.

Tables VI and VI show the standard deviations of returns for 100 (Table VI)
and 20000 (Table VI) tick times for the four cases, respectively. The standard
deviations are larger when both the CTAA and STRTA exist. The larger variation
of prices leads to an increase in frequency that prices reach the maximum or
minimum and that the CTAA and STRTA have more chances to trade.

Fig. 3 shows the mechanism derived from the aforementioned discussion. The
CTAA and STRTA place a market order, which leads to a waiting order being
removed from the order book. Decreasing waiting orders leads to other buy
(sell) market orders executing with a waiting sell (buy) order of a higher
(lower) price. This leads to a larger variation of prices and an increase in
frequency that prices reach the maximum or minimum. Furthermore, the CTAA and
STRTA have more chances to trade and earn more when each other exist, and
increasing their trades leads to more market orders and the strengthening of
this cycle.

##  IV Summary

In this study, I built an artificial market model by adding a CTAA and STRTA
to the prior model of Mizuta et al.[13] and investigated whether emerging
STRTAs led to a decrease in CTAA revenue to determine whether STRTs prey on
CTAs for profit.

To the contrary, my results showed that CTAAs and STRTAs have more chances of
trade and earn more when each other exist. Therefore, they have a mutually
beneficial relationship while some people says STRTs prey on CTAs for profit.

The mechanism of mutually beneficial relationship is following. The CTAA and
STRTA place a market order, which leads to a waiting order being removed from
the order book. Decreasing waiting orders leads to other buy (sell) market
orders executing with a waiting sell (buy) order of a higher (lower) price.
This leads to a larger variation of prices and an increase in frequency that
prices reach the maximum or minimum. Furthermore, the CTAA and STRTA have more
chances to trade and earn more when each other exist, and increasing their
trades leads to make more market orders and the strengthening of this cycle.

When one investment strategy fails to work, new strategies, especially those
that use new technologies or faster tradings, have tended to be blamed for a
long time. Therefore, such arguments that new strategies are used to impede
existing strategies should be discussed with particular suspicion. Since the
second half of the 20th century, this tendency has become stronger, and the
media have often criticized new investment strategies using new technologies,
electronic trades, program trades, algorithm trades, high-frequency trades,
and AI traders. However, these critics seem to be missing the point. These
criticisms should be investigated using an artificial market model, which are
considered future works.

TABLE VII: Stylized facts | CTAA | exist | not exist
---|---|---|---
| STRTA | exist | not exist | exist | not exist
kurtosis or returns | $3.46$ | $3.80$ | $3.67$ | $3.99$
| lag |  |  |  |
| 1 | $0.131$ | $0.127$ | $0.138$ | $0.128$
autocorrelation | 2 | $0.081$ | $0.084$ | $0.078$ | $0.080$
coefficient for | 3 | $0.061$ | $0.069$ | $0.059$ | $0.060$
square returns | 4 | $0.054$ | $0.058$ | $0.047$ | $0.051$
| 5 | $0.046$ | $0.053$ | $0.042$ | $0.043$

## Appendix

###  IV-A Basic concept for constructing a model

An artificial market model, which is a kind of agent-based model, can be used
to investigate situations that have never occurred, handle regulation changes
that have never been made, and isolate the pure contribution of these changes
to price formation and liquidity [8, 9]. These are the advantages of an
artificial market simulation.

However, the outputs of this simulation would not be accurate or be credible
forecasts of the actual future. The simulation needs to reveal possible
mechanisms that affect price formation through many simulation runs, e.g.,
searching for parameters or purely comparing the before and after states of
changes. The possible mechanisms revealed by these runs provide new
intelligence and insight into the effects of the changes on price formation in
actual financial markets. Other methods of study, e.g., empirical studies,
would not reveal such possible mechanisms.

Artificial markets should replicate the macro phenomena that exist generally
for any asset at any time. Price variation, which is a kind of macro
phenomenon, is not explicitly modeled in artificial markets. Only micro
processes, agents (general investors), and price determination mechanisms
(financial exchanges) are explicitly modeled. Macro phenomena emerge as the
outcome interactions of micro processes. Therefore, the simulation outputs
should replicate existing macro phenomena to generally prove that simulation
models are probable in actual markets.

However, it is not the primary purpose for an artificial market to replicate
specific macro phenomena only for a specific asset or a specific period.
Unnecessary replication of macro phenomena leads to models that are overfitted
and too complex. Such models would prevent us from understanding and
discovering mechanisms that affect price formation because the number of
related factors would increase.

In addition, artificial market models that are too complex are often
criticized because they are very difficult to evaluate[5]. A model that is too
complex not only would prevent us from understanding mechanisms but also could
output arbitrary results by overfitting too many parameters. It is more
difficult for simpler models to obtain arbitrary results, and these models are
easier to evaluate.

Therefore, I constructed an artificial market model that is as simple as
possible and does not intentionally implement agents to cover all the
investors who would exist in actual financial markets.

As Michael Weisberg mentioned[17], “Modeling, (is) the indirect study of real-
world systems via the construction and analysis of models.” “Modeling is not
always aimed at purely veridical representation. Rather, they worked hard to
identify the features of these systems that were most salient to their
investigations.”

Therefore, effective models are different depending on the phenomena they
focus on. Thus, my model is effective only for the purpose of this study and
not for others. The aim of my study is to understand how important properties
(behaviors, algorithms) affect macro phenomena and play a role in the
financial system rather than representing actual financial markets precisely.

As Weisberg also mentioned[17], “When one invokes a computational model to
explain some phenomenon, one is typically using transition rules or algorithm
as the explanans. Schelling explained segregation by pointing out that small
decisions reflecting small amounts of bias will aggregate to massively
segregated demographics[18]. Neither the time sequence of the model’s states
nor the final, equilibrium state of the model carries the explanatory force;
the algorithm itself is needed.”

An aim of my study is to understand how important properties (behaviors,
algorithms) affect the investigation of macro phenomena and play a role in the
financial system rather than representing actual financial markets precisely.

The above discussion holds not only for artificial markets but also for agent-
based models used in other fields besides financial markets. For example,
Thomas Schelling, who was received the Nobel prize in economics, used an
agent-based model to discuss the mechanism of racial segregation. The model
was built very simply compared with an actual town in order to focus on the
mechanism[18]. Indeed, while it was not able to predict the segregation
situation in the actual town, it was able to explain the mechanism of
segregation as a phenomenon.

Harry Stevens, a newspaper writer, simulated an agent-based model to explain
the spread of COVID-19 and to find a way how to prevent infection [19]. The
model is too simple to replicate the real world, but its simplicity enabled it
to reveal the mechanism behind the spread.

Michael Weisberg studied what mathematical and simulation models are in the
first place and cited the example of a map[17]. Needless to say, a map models
geographical features on the way to a destination. With a simple map we can
easily understand the way to the destination. On the other hand, while a
satellite photo replicates actual geographical features very well, we can not
easily find the way to the destination.

The title page of the book [17] cited a passage from a short story by Jorge
Borges[20], “In time, those Unconscionable Maps no longer satisfied, and the
Cartographers Guilds struck a Map of the Empire whose size was that of the
Empire, and which coincided point for point with it . . . In the Deserts of
the West, still today, there are Tattered Ruins of that Map, inhabited by
Animals and Beggars.” The story in which a map was enlarged to the same size
as the real Empire to become the most detailed of any map is an analogy to
that too detailed a model is not useful. This story give us one of the most
important lessons for when we build and use any model.

###  IV-B Validation of the model

In many previous artificial market studies, the models were validated to
determine whether they could explain stylized facts, such as a fat tail or
volatility clustering [4, 5, 8, 9]. A fat tail means that the kurtosis of
price returns is positive. Volatility clustering means that square returns
have a positive autocorrelation, and this autocorrelation slowly decays as its
lag becomes longer.

Many empirical studies, e.g., that of Sewell [21], have shown that both
stylized facts (fat tail and volatility clustering) exist statistically in
almost all financial markets. Conversely, they also have shown that only the
fat tail and volatility clustering are stably observed for any asset and in
any period because financial markets are generally unstable. This leads to the
conclusion that an artificial market should replicate macro phenomena that
exist generally for any asset at any time, fat tails, and volatility
clustering. This is an example of how empirical studies can help an artificial
market model.

The kurtosis of price returns and the autocorrelation of square returns are
stably and significantly positive, but the magnitudes of these values are
unstable and very different depending on the asset and/or period. Very broad
magnitudes of about $1\sim 100$ and about $0\sim 0.2$, respectively, have been
observed [21].

For the aforementioned reasons, an artificial market model should replicate
these values as significantly positive and within a reasonable range. It is
not essential for the model to replicate specific values of stylized facts
because the values of these facts are unstable in actual financial markets.

Table VII lists the statistics showing the stylized facts, kurtosis of price
returns for $100$ tick times, and autocorrelation coefficient for square
returns for $100$ tick times for the four cases of whether CTAAs and STRTAs
exist. This shows that this model replicated the statistical characteristics,
fat-tails, and volatility clustering observed in real financial markets.

## References

  * [1] F. A. Clenow, _Following the Trend: Diversified Managed Futures Trading_. Wiley, 2012.
  * [2] S. Nagao and M. Kobayashi, “Developmental case study of hedge fund clone program as knowledge creation means: Action research using soft systems methodology and social learning,” _Review of Integrative Business and Economics Research_ , vol. 5, no. 4, pp. 307–328, 2016.
  * [3] B. L. Raschke and A. L. Connors, _Street Smarts: High Probability Short-Term Trading Strategies_. M. Gordon Publishing Group, 1996.
  * [4] B. LeBaron, “Agent-based computational finance,” _Handbook of computational economics_ , vol. 2, pp. 1187–1233, 2006.
  * [5] S.-H. Chen, C.-L. Chang, and Y.-R. Du, “Agent-based economic models and econometrics,” _Knowledge Engineering Review_ , vol. 27, no. 2, pp. 187–219, 2012.
  * [6] A. Todd, P. Beling, W. Scherer, and S. Y. Yang, “Agent-based financial markets: A review of the methodology and domain,” in _2016 IEEE Symposium Series on Computational Intelligence (SSCI)_ , 2016, pp. 1–5.
  * [7] T. Mizuta, “A brief review of recent artificial market simulation (agent-based model) studies for financial market regulations and/or rules,” _SSRN Working Paper Series_ , 2016.
  * [8] T. Mizuta, “An agent-based model for designing a financial market that works well,” in _2020 IEEE Symposium Series on Computational Intelligence (SSCI)_ , 2020, pp. 400–406.
  * [9] T. Mizuta, “Artificial intelligence (ai) for financial markets: a good ai for designing better financial markets and a bad ai for manipulating markets,” _Digital Designs for Money, Markets, and Social Dilemmas, Springer, In press_ , 2022.
  * [10] J. D. Farmer and D. Foley, “The economy needs agent-based modelling,” _Nature_ , vol. 460, no. 7256, pp. 685–686, 2009.
  * [11] S. Battiston, J. D. Farmer, A. Flache, D. Garlaschelli, A. G. Haldane, H. Heesterbeek, C. Hommes, C. Jaeger, R. May, and M. Scheffer, “Complexity theory and financial regulation,” _Science_ , vol. 351, no. 6275, pp. 818–819, 2016.
  * [12] R. Bookstaber, _The End of Theory: Financial Crises, the Failure of Economics, and the Sweep of Human Interaction_. Princeton University Press, 2017.
  * [13] T. Mizuta, S. Kosugi, T. Kusumoto, W. Matsumoto, K. Izumi, I. Yagi, and S. Yoshimura, “Effects of price regulations and dark pools on financial market stability: An investigation by multiagent simulations,” _Intelligent Systems in Accounting, Finance and Management_ , vol. 23, no. 1-2, pp. 97–120, 2016.
  * [14] C. Chiarella and G. Iori, “A simulation analysis of the microstructure of double auction markets,” _Quantitative Finance_ , vol. 2, no. 5, pp. 346–353, 2002.
  * [15] L. Menkhoff and M. P. Taylor, “The obstinate passion of foreign exchange professionals: Technical analysis,” _Journal of Economic Literature_ , vol. 45, no. 4, pp. 936–972, 2007.
  * [16] R. Yamamoto, “Predictor choice, investor types, and the price impact of trades on the tokyo stock exchange,” _Computational Economics_ , 2021.
  * [17] M. Weisberg, _Simulation and Similarity: Using Models to Understand the World_. Oxford Studies in the Philosophy of Science, 2012.
  * [18] T. C. Schelling, _Micromotives and Macrobehavior_. W. W. Norton & Company, 2006.
  * [19] H. Stevens, “Why outbreaks like coronavirus spread exponentially, and how to ”flatten the curve”,” 2020. [Online]. Available: https://www.washingtonpost.com/graphics/2020/world/corona-simulator/
  * [20] J. L. Borges, “Del rigor en la ciencia,” in _Historia universal de la infamia_. Emecé, 1954.
  * [21] M. Sewell, “Characterization of financial time series,” _Research Note, University College London, Department of Computer Science_ , no. RN/11/01, 2011\.