# Discrete time crystals in the presence of non-Markovian dynamics

Bandita Das  Department of Physical Sciences, Indian Institute of Science
Education and Research Berhampur, Berhampur 760010, India  Noufal Jaseem
Department of Physical Sciences, Indian Institute of Science Education and
Research Berhampur, Berhampur 760010, India  Victor Mukherjee  Department of
Physical Sciences, Indian Institute of Science Education and Research
Berhampur, Berhampur 760010, India

###### Abstract

We study discrete time crystals (DTCs) in periodically driven quantum systems,
in the presence of non-Markovian dissipation. In contrast to DTCs observed in
earlier works in the presence of Markovian dynamics, using the open Dicke
model in presence of Jaynes-Cummings-like dissipation, we show that non-
Markovian regime can be highly beneficial for stabilizing DTCs over a wide
range of parameter values. This may be attributed to periodically varying
dissipation rates even at long times in the case of non-Markovian dynamics.
Further the Markovian and non-Markovian regimes show sharp distinctions for
intermediate strengths of the dissipator coefficient, with a time-independent
steady-state in the Markovian regime being replaced by varied dynamical
phases, including DTC order, in the non-Markovian regime. We also verify the
robustness of the DTC phase in the non-Markovian regime by introducing errors
both in the Hamiltonian as well as in the dissipation. Our study shows the
possibility of using DTC as a probe for non-Markovian dynamics in periodically
modulated open quantum systems, at long times.

##  I Introduction

Many-body quantum systems driven out of equilibrium present exciting fields of
research, owing to the different non-trivial behaviors exhibited by such
systems, such as Kibble-Zurek mechanism [1, 2], dynamical localization [3, 4],
and time-translational symmetry breaking [5, 6], to name a few. Recent
developments in quantum simulators have enabled researchers to experimentally
study the behavior of many-body systems driven out of equilibrium as well [7,
8]. Consequently, theoretical and experimental research on dynamics of many-
body systems have received a lot of attention in the recent years. In
particular, time-translational symmetry breaking in the form of continuous [9,
10, 11, 12, 13] or discrete [14, 15, 16] time crystals have received wide
interest from the community in the last decade, and also realised
experimentally [17, 18, 19, 20]

Discrete time crystals (DTCs) are associated with spontaneous symmetry
breaking in time, and are formed in periodically driven many-body quantum
systems. In closed quantum systems DTCs have been studied widely, both
theoretically [21, 14, 22, 23, 15] as well as experimentally [7], in many-body
localized systems. More recently, DTCs have been shown to exist in
periodically driven clean quantum systems as well [24, 25], and also in the
presence of dissipation [26, 27, 28]. Such dissipative DTCs in periodically
driven many-body open quantum systems can be viewed as quantum engines,
wherein a part of the energy supplied through periodic drive flows to a cold
bath, while the rest is obtained as output work [29, 6]. Analogous to
different phases of matter, DTCs are robust to small perturbations. However,
large perturbations, such as errors in the Hamiltonian [30], or strong rates
of dissipation [27], can destroy a DTC. Consequently, finding scenarios which
can result in stable DTCs is a fundamental question in the field of time-
translational symmetry breaking in many-body systems. Furthermore, open
questions remain regarding the existence and behavior of DTCs for different
types of dissipative dynamics. Here we address the above two crucial issues by
focusing on a periodically modulated many-body quantum system in the presence
of non-Markovian dynamics. We show that suitably controlled non-Markovian
dynamics can result in the generation of DTCs as well as other dynamical
phases.

Dissipative dynamics can be classified as Markovian or non-Markovian depending
on the absence or presence of memory effects, respectively [31]. Markovian
dynamics can be described by time-independent Lindblad superoperators, such
that the system approaches a long-time steady-state monotonically. On the
other hand, non-Markovian dynamics are associated with memory effects, which
may result in the system moving away from the steady state for some intervals
of time [32, 33, 34, 35].

In this work, we focus on an open Dicke model, and present results in the
Markovian limit, following Ref. [30]. In case of Markovian dynamics, DTC phase
obtained for weak disspation is replaced by a time-independent steady-state
(TISS) for intermediate or higher rates of dissipation [27]. In order to study
the dynamics in the non-Markovian regime, we focus on a damped Jaynes-
Cummings-like dissipator with tunable parameters, which allows one to traverse
between the Markovian and the non-Markovian regimes [31]. We show that in
contrast to the Markovian regime, the non-Markovian regime can be associated
with a stable DTC for a wide range of parameter values, thereby significantly
expanding the regime of parameters allowing the existence of time
translational symmetry breaking. We present phase diagrams for the Markovian
and non-Markovian regimes; striking differences emerge in the response of a
system in the presence of Markovian and non-Markovian dissipation, for
intermediate strengths of the dissipator coefficient. This presents the
intriguing possibility of using dissipative DTC as a probe for detecting non-
Markovianity at long times. Furthermore, we show that the DTC phase persists
in the non-Markovian regime even in the presence of a random noise in the
dissipator.

We begin by discussing the model and dynamics Sec. II; we start by discussing
the dynamics in a generic open quantum system in Sec. II.1, and then focus on
the specific example of an open Dicke model in Sec. II.2. We consider
Markovian dynamics in Sec. II.2.1, while we address the non-Markovian regime
in Sec. II.2.2, and dissipation in the presence of a random noise in Sec.
II.2.3. Finally we conclude in Sec. III.

##  II Model and dynamics

###  II.1 Dissipative discrete time crystals

We consider a system $\mathcal{S}$ described by generic periodically modulated
Hamiltonian $H(t)=H(t+T)$. The system is coupled to a dissipative environment,
such that the state $\rho$ of the system evolves following the master equation

$\displaystyle\dot{\rho}=-\frac{i}{\hbar}[H(t),\rho(t)]+\kappa(t)\mathcal{L}[\rho(t)].$
(1)

Here $\mathcal{L}$ is a Lindblad superoperator, while $\kappa(t)$ determines
the rate of thermalization with a bath, which we assume to be time-dependent
in general. Complete possitivity demands $\int^{t}_{0}\kappa(t)dt\geq 0$ for
all $t$. Furthermore, in case of Markovian dynamics, $\kappa(t)\geq 0$ for all
times, which ensures that the system approaches the long-time steady state at
all times. As shown in [30] for a time-independent $\kappa(t)=\kappa_{0}$
$\forall~{}t$, a periodic modulation in $H(t)$ can result in a DTC in case of
Markovian dynamics, for small values of $\kappa_{0}$. On the other hand, a
large $\kappa_{0}$ is associated with rapid thermalization with the bath in a
time scale $\tau_{\rm th}\sim~{}\kappa_{0}^{-1}$, and a destruction of the DTC
phase [27].

In contrast to Markovian dynamics, $\kappa(t)$ can assume negative values for
some time intervals in the case of non-Markovian dynamics [32]. This results
in the non-Markovian regime being associated with information backflow, such
that a negative $\kappa(t)$ may drive $\mathcal{S}$ away from the long-time
steady state for some time-intervals [35]. The distinct properties of
Markovian and non-Markovian dynamics raises questions regarding the existence
and characteristics of DTCs in the non-Markovian regime. Below we focus on the
specific setup of a Dicke model in the presence of a Jaynes-Cummings-like
dissipation, to show that indeed, the behavior of DTC in the non-Markovian
regime can be significantly different from that seen in the Markovian regime.

###  II.2 Modulated open Dicke model

We consider an open Dicke model comprising $N$ two-level atoms in a cavity,
which is described by the Hamiltonian [36]

$\displaystyle{\hat{H}(\lambda)=\omega\hat{a}^{\dagger}\hat{a}+\omega_{0}\hat{J}_{z}+\frac{2\lambda_{t}}{\sqrt{N}}(\hat{a}+\hat{a}^{\dagger})\hat{J}_{x}}.$
(2)

Here $a$ and $a^{\dagger}$ are respectively the Bosonic annihilation and
creation operator for the photons, $J_{\mu}=\sum_{i=1}^{N}\sigma_{i}^{\mu}$,
where $\sigma_{i}^{\mu}$ denotes the Pauli matrix corresponding to the $i$th
spin along the $\mu=x,y,z$ axis, $\omega$ is the frequency of photon field,
$\omega_{0}$ represents the transition frequency of the two level atoms, and
$\lambda_{t}$ is the atom-photon coupling strength, which we consider to be
time-dependent in general.

The setup considered here posseses a $\mathcal{Z}_{2}$ symmetry; the
Hamiltonian (2) commutes with the parity operator
$P=e^{i\pi(\hat{a}^{\dagger}\hat{a}+\hat{J_{z}}+N/2)}$, such that
$\hat{H}(\lambda)$ remains invariant under the transformation
$\hat{a}\rightarrow-\hat{a},$ and $\hat{J_{x}}\rightarrow-\hat{J_{x}}$.
However, in the thermodynamic limit $N\to\infty$, the $\mathcal{Z}_{2}$
symmetry is spontaneously broken leading to superradiant phase transition at
the critical value of
$\lambda=\lambda_{c}=\frac{1}{2}\sqrt{(\omega_{0}/\omega)(\omega^{2}+\kappa_{0}^{2}/4)}$
[37, 30].

In order to study the possibility of time-translational symmetry breaking in
the presence of non-Markovian dynamics, we now consider a phenomenological
model, wherein the above atom-photon setup (see Eq. (2)) evolves in the
presence of a dissipative environment, following the master equation (1) [31].
In analogy with the damped Jaynes-Cummings model, here we focus on a time
dependent photon decay rate $\kappa(t)$ given by

${\kappa(t)}=\begin{cases}\frac{2m\kappa_{0}\sinh(td/2)}{d\cosh(td/2)+m\sinh(td/2)}&\text{$|\kappa(t)|<\kappa_{max}$}\\\
\kappa_{max}&\text{ $|\kappa(t)|\geq\kappa_{max}$},\\\ \end{cases}$ (3)

where $d=\sqrt{m^{2}-2m\kappa_{0}}$ and
$\mathcal{L}[\hat{a}]\hat{\rho}=\hat{a}\hat{\rho}\hat{a}^{\dagger}-\frac{1}{2}\\{{\hat{a}^{\dagger}\hat{a},\hat{\rho}}\\}$
(see Eq. (1)). Here $\kappa_{max}>\kappa_{0}$ is a parameter which can be
tuned to control the maximum possible rate of dissipation. In case of the
damped Jaynes-Cummings model, $\kappa_{max}\to\infty$, $m$ denotes the
spectral width of the bath, while $\kappa_{0}$ is related to the system-bath
coupling strength [31].

The above form of $\kappa(t)$ (Eq. (3)) allows us to tune between the
Markovian regime ($\kappa_{0}<m/2$) in which case $\kappa(t)>0~{}\forall~{}t$,
and the non-Markovian regime ($\kappa_{0}>m/2$), in which case $\kappa(t)$
assumes an oscillatory form realized by replacing $\sinh$ ($\cosh$) by $\sin$
($\cos$) in Eq. (3), and can can take negative values for some time intervals
[31, 35].

Following Ref. [30], one can use Eq. (1) to arrive at time-evolution equations
for the scaled variables
$x=\langle\hat{a}+\hat{a}^{\dagger}\rangle/\sqrt{2N\omega}$,
$p=i\langle\hat{a}-\hat{a}^{\dagger}\rangle/\sqrt{2N\omega}$, and
$\mathbb{j}=(j_{x},j_{y},j_{z})$ with $j_{\mu}=\langle\hat{J_{\mu}}\rangle/N$,
in thermodynamic limit of $N\gg 1$ (see Appendix).

####  II.2.1 Markovian regime

We start by focusing on the extreme Markovian limit of
$m\gg\kappa_{0},~{}m\to\infty$, in which case $\kappa(t)$ reduces to the time-
independent form $\kappa(t)\approx\kappa_{0}$ for all times. For a time-
independent $\lambda_{t}=\lambda_{0}$, the above setup has two-symmetry broken
steady states $\rho_{ss}$ and $\rho_{ss}^{\prime}$ (see Appendix) [38, 30]

In order to study the emergence of time-translational symmetry breaking in the
above system, we introduce a period doubling dynamics aimed at periodically
evolving between $\rho_{ss}$ and $\rho_{ss}^{\prime}$. To that end, we apply a
periodically modulated $\lambda_{t}$ of the form

${\lambda_{t+T}=\lambda_{t}}=\begin{cases}\lambda_{0}&\text{$0\leq t<T/2$}\\\
0&\text{ $T/2\leq t<T$}\\\ \end{cases}$ (4)

where, $T=\frac{2\pi}{\omega_{T}}$ , $\lambda_{0}>\lambda_{c}$ and
$\omega_{T}=d/2$ [30]. In the resonant case $\omega=\omega_{0}$ and in absence
of dissipation ($\kappa_{0}=0$), the above form of $\lambda_{t}$ results in
the parity operation $\hat{P}=\exp[-i(T/2)H(0)+i\pi N/2]$ during the second
half period $T/2\leq t<T$, up to a global phase which we can ignore.

We start with the system in one of the symmetry broken steady states
$\rho_{s}$ at the start of a time period. The system stays at $\rho_{s}$
during the first half of the time period, during which time
$\lambda_{t}=\lambda_{0}$. Thereafter, the modulation in $\lambda$ results in
a parity operation $\hat{P}$ during the second half period, which, in the
absence of a dissipative bath (i.e., $\kappa_{0}=0$), would take the system to
the other symmetry broken steady state
$\rho_{s}^{\prime}=\hat{P}\rho_{s}\hat{P}^{\dagger}$ at the end of the second
half period. Consequently, the above modulation Eq. (4) results in
$\rho(t)=\rho_{s}$ ($\rho(t)=\rho_{s^{\prime}}$) at the end of even (odd)
number of time periods for $\kappa_{0}=0$, thus giving rise to a DTC behavior.

One can study the robustness of the DTC phase by introducing an error
parameter $\epsilon$, defined by
$\omega=\left(1-\epsilon\right)\omega_{T};~{}\omega_{0}=\left(1+\epsilon\right)\omega_{T}$.
As shown in Ref. [30], the DTC order is robust to small values of $\epsilon$,
while a larger $\epsilon$ may result in destruction of the same.

Figure 1: The mean magnetizations $j_{x}$, $j_{z}$, $j_{z}$ and the
corresponding Bloch sphere representations are shown for the Markovian with
$m\to\infty$ (a, b, c, d) and non-Markovian (e, f) regimes. We get a DTC phase
for (a) small $\kappa_{0}=0.05$, represented by two dots in the (b) Bloch
sphere. In contrast, a TISS is obtained for large $\kappa_{0}$ (c, d) in the
Markovian regime. The DTC phase is preserved for intermediate values of
$\kappa_{0}$ in the non-Markovian regime (e, f). Here $\omega_{T}=1$,
$T=2\pi$, $\epsilon=0.02$, $\lambda_{0}=1$, a) $\kappa_{0}=0.05$, c)
$\kappa_{0}=2.7$ c) $\omega_{T}=d/2$, $\kappa_{0}=2.7$ and $m=\kappa_{0}/4$,
$T=2\pi/\omega_{T}$.

Let us now focus on the behavior of the DTC order with increasing values of
the dissipator coefficient $\kappa_{0}$. The parameter $\kappa_{0}$ sets the
time-scale of thermalization $\tau_{th}\sim\kappa_{0}^{-1}$. Consequently, a
low $\kappa_{0}$ (i.e., $\tau_{th}\gg T$) may be expected to facilitate the
emergence of DTC (see Fig. 1a), as can also be verified in the Bloch sphere
representation of the stroboscopic dynamics, where a DTC corresponds to two
distinct states, shown by the blue dots in Fig. 1b. On the other hand, in the
limit of large $\kappa_{0}$ (i.e., $\tau_{th}\ll T$), the dissipative
mechanism may be expected to dominate the dynamics, thus leading to
destruction of the DTC phase, and emergence of a TISS (see Figs. 1c - 1d)
[27]. Therefore it is crucial to study scenarios where a DTC order might be
robust over a wide range of parameter values. Below, we go beyond the
Markovian approximation to show that remarkably, non-Markovian dynamics may
allow us to achieve the above aim, thus resulting in DTC which is robust to
both $\epsilon$ and $\kappa_{0}$.

####  II.2.2 Non-Markovian regime

In this section, we focus on the non-Markovian regime, obtained for
$\kappa_{0}>m/2$. As discussed above, this regime corresponds to a
periodically varying $\kappa(t)$ with $\kappa(t)<0$ for some time intervals
(see Eq. (3)), which eventually results in the so called information backflow
in the system [39].

In contrast to the behavior reported for Markovian dynamics [30] where DTC
phase is present only for small $\kappa_{0}$ (see Figs. 1a - 1d), numerical
analysis shows that non-Markovian dynamics makes DTC more robust against
$\kappa_{0}$, as signified by the presence of time-crystalline order for
intermediate values of $\kappa_{0}$ and $\kappa_{max}$; this is verified both
for the Markovian limit with time-independent $\kappa(t)=\kappa_{0}$, obtained
for $m/\kappa_{0}\to\infty$ (see Fig. 1), as well as for finite values of
$m/\kappa_{0}$ (see Figs. 2 and 3). This robustness of the DTC phase w.r.t.
$\kappa_{0}$ may be attributed to the periodically varying $\kappa(t)$ in the
non-Markovian regime, even at long times [31, 40]. However, we note that the
DTC phase is replaced by a TISS or a thermal phase for large $\kappa_{max}$
(see Appendix).

Figure 2: Stroboscopic dynamics (a,b,c,d), of the atomic average angular
momentum $j_{\mu}$ for atom-photon coupling $\lambda=1$ and time dependent
photon-loss rate $\kappa(t)$. As $\epsilon$ is varied, various dynamical
phases emerges over long periods of time. Only the last 30 periods is shown in
the figure. Last 200 periods (blue) is projected onto the bloch sphere. Here
$\kappa_{0}=2.7$ and $\kappa_{max}=5$. Figure 3: Phase diagram representing
TISS, DTC, other dynamical phases (DP) and thermal phase (TP) as functions of
$m$ and and $\kappa_{0}=2.7a)\kappa_{max}=5,b)\kappa_{max}=3$ , the dashed
line at $m=2\kappa_{0}$, represents the transition line from NM to M regime.

In addition to robustness w.r.t. large $\kappa_{0}$, the DTC phase shows
resilience even in the presence of a non-zero detuning parameter $\epsilon$.
Analogous to that seen in the Markovian regime [30], as we vary $\epsilon$,
several dynamical phases emerge, viz. thermal (Fig. 2a), DTC (Fig. 2b), sextet
(Fig. 2c) and limit cycle (Fig. 2d) for these parameters.

In order to have a deeper understanding of the effect of non-Markovian
dynamics on the presence of DTC, we study the behavior of the system for
different values of $\epsilon$ and $m/\kappa_{0}$. The results are summarized
in the phase diagram Fig. 3. We navigate between the Markovian
($m/\kappa_{0}>2$) and the non-Markovian ($m/\kappa_{0}<2$) regimes by varying
$m$, for a constant $\kappa_{0}$. Interestingly, as seen in Fig. 3, one can
clearly distinguish the Markovian and non-Markovian regimes from the stark
difference in behavior of the system across the transition. For a large enough
$\kappa_{0}$ ($\kappa_{0}\gg T^{-1}$), the system approaches a TISS in the
Markovian regime. However, varied dynamical phases, including DTC, appear on
undergoing transition to the non-Markovian regime at $m=2\kappa_{0}$
(highlighted by a dashed black line in Fig. 3). Interestingly, the DTC phase
shows the most robustness w.r.t. $\epsilon$ close to the Markovian to non-
Markovian transition line $m=2\kappa_{0}$. This might be owing to small values
of the parameter $d$ for $m\to 2\kappa_{0}^{-}$, which results in a slowly
varying $\kappa(t)$, and an effectively longer time scales $\sim~{}d^{-1}$. As
we move away from the transition line (i.e., $m\ll 2\kappa_{0}$), we find rich
dynamical phases as we vary $\epsilon$. However, we note that this stark
difference between Markovian and non-Markovian regimes vanishes for small
$\kappa_{0}$, in which case varied dynamical phases may exist for different
values of $\epsilon$ in the Markovian regime as well.

####  II.2.3 Non-Markovian regime with noisy aperiodic $\kappa(t)$

We next focus on the question whether the DTC phase in the non-Markovian
regime is stable in the presence of random noise in $\kappa(t)$. To this end,
we consider an aperiodic $\kappa(t)$, realized by introducing random
fluctuations in the dissipator coefficient, given by
$\kappa(t)^{\prime}=\kappa(t)+a_{0}f(t)$.

$\displaystyle\kappa(t)^{\prime}=\kappa(t)+a_{0}f(t).$ (5)

Here $f(t)$ ($-1\leq f(t)\leq 1$) is a random function applied at every time
step, while $a_{0}$ gives the strength of the fluctuations.

Figure 4: Figure showing (a) an aperiodic $\kappa(t)$, formed through the
introduction of random fluctuations in the dissipation (see Eq. (5)), and (b)
the corresponding DTC phase. Here
$\kappa_{0}=2.7,m=\kappa_{0}/4,\epsilon=0.03$,$a_{0}=0.5$ and $\kappa_{max}=5$

As shown in Fig. 5, the DTC order persists in this case, thereby showing the
robustness of DTC phase in the non-Markovian regime, in the presence of
intermediate values of $\kappa_{0}$, non-zero detuning error $\epsilon$ as
well as random noise $f(t)$ in the dissipator coefficient.

##  III Conclusion

We study the emergence of DTC in the presence of non-Markovian dynamics by
considering an open Dicke model with Jaynes-Cummings-like dissipation
$\kappa(t)$ (Cf. Eq. (3)). Our analysis shows that in contrast to Markovian
dynamics, the DTC phase is more robust to a wide range of parameter values in
the non-Markovian regime ($m<2\kappa_{0}$), which may be attributed to
periodically varying $\kappa(t)$ even at long times. We present a dynamical
phase diagram w.r.t. $\epsilon$ and $m$ (Cf. Fig. 3). The transition from
Markovian to non-Markovian regime at $m=2\kappa_{0}$ is marked by TISS for all
values of $\epsilon$ in the Markovian regime for large enough $\kappa_{0}$,
which changes to varied dynamical phases, including DTC, in the non-Markovian
regime. Furthermore, the DTC phase emerging for $m$ close to $2\kappa_{0}$
shows substantial robustness w.r.t. $\epsilon$, while this robustness is
reduced, even though does not vanish, as we move deeper into the non-Markovian
regime of $m<2\kappa_{0}$. Furthermore, we verify the presence of robust DTC
order in the non-Markovian regime even for aperiodic $\kappa(t)$, realized by
introducing a random noise in the dissipator coefficient.

We note that the results presented here shows the possibility of using DTC as
a probe for non-Markovian dynamics at long times, for intermediate strengths
of the dissipator coefficient $\kappa_{0}$. In the absence of an external
modulation, the steady-state of a dynamics may be identical for the Markovian
and the non-Markovian dynamics [31, 35]. Consequently, in general a probe for
differentiating between the Markovian and non-Markovian regimes is applicable
only for short times, before the system reaches the steady-state, which can be
highly challenging in cases of short thermalization times. In contrast, as
discussed here, the long-time state reached by a system can change drastically
in the presence of a periodic modulation; for intermediate values of
$\kappa_{0}$ we get a TISS in case of Markovian dynamics, which changes to a
DTC phase for some ranges of $\epsilon$ in the non-Markovian regime (see Figs.
3). This may provide us with a novel way of estimating the nature of the
dynamics even at long times, when conventional probes acting in absence of
periodic modulation will fail.

We expect the dissipative DTCs studied here can be realized experimentally in
currently existing setups. For example, the presence of DTCs has been verified
experimentally in different platforms, both in closed quantum systems, such as
in ion traps [7], and in the presence of dissipation, for example in optical
cavities [41]. Experimental studies of DTCs in the presence of different forms
of dissipation would require control over bath spectral functions, which for
example can be achieved through the introduction of filters [42].

Our results show that non-Markovian dynamics can be highly relevant for time-
translational symmetry breaking in many-body open quantum systems, and can
show features distinct from those seen in the Markovian regime. This also
raises open questions regarding the fate of DTCs in more generic forms of non-
Markovian dynamics, and regarding the possible role played by information
backflow in the formation of DTCs.

## ACKNOWLEDGEMENTS

B.D. acknowledges support from Prime Minister Research Fellowship(PMRF). V.M.
acknowledges support from Science and Engineering Research Board (SERB)
through MATRICS (Project No. MTR/2021/000055) and Seed Grant from IISER
Berhampur. V.M. also acknowledges Rosario Fazio and Manabendra Bera for
helpful discussions.

Appendix

*

##  Appendix A Evolution equations and steady states

As discussed in Ref. [30], one can use Eq. (1) to arrive at the following
time-evolution equations for the scaled variables
$x=\langle\hat{a}+\hat{a}^{\dagger}\rangle/\sqrt{2N\omega}$,
$p=i\langle\hat{a}-\hat{a}^{\dagger}\rangle/\sqrt{2N\omega}$, and
$\mathbb{j}=(j_{x},j_{y},j_{z})$ with $j_{\mu}=\langle\hat{J_{\mu}}\rangle/N$,

$\displaystyle\frac{d\textbf{j}}{dt}$ $\displaystyle=$
$\displaystyle(-\omega_{0}\textbf{e}_{z}+2\lambda_{t}\sqrt{2\omega}x\textbf{e}_{x})\times\textbf{j}$
(6) $\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle
p-\frac{\kappa(t)}{2}x$ $\displaystyle\frac{dp}{dt}$ $\displaystyle=$
$\displaystyle-\omega^{2}x-\frac{\kappa(t)}{2}p-2\lambda_{t}\sqrt{2\omega}xj_{x}.$

in thermodynamic limit of $N\gg 1$.

For a time-independent $\lambda_{t}=\lambda_{0}$, the setup considered here
has two-symmetry broken steady states $\rho_{ss}$ and $\rho_{ss}^{\prime}$
[38, 30]

$\displaystyle j_{x}$ $\displaystyle=$
$\displaystyle\pm\sqrt{1-\frac{\lambda_{c}^{4}}{\lambda^{4}}},~{}~{}j_{y}=0,~{}~{}j_{z}=\frac{-\lambda_{c}^{2}}{\lambda^{2}}$
$\displaystyle x$ $\displaystyle=$
$\displaystyle\mp\lambda\frac{\sqrt{2\omega(1-\frac{\lambda_{c}^{4}}{\lambda^{4}})}}{\omega^{2}+\kappa_{0}^{2}/4},$
$\displaystyle p$ $\displaystyle=$
$\displaystyle\mp\kappa_{0}/2\frac{\sqrt{2\omega(1-\frac{\lambda_{c}^{4}}{\lambda^{4}})}}{\omega^{2}+\kappa_{0}^{2}/4}.$
(7)

##  Appendix B Dynamics for large $\kappa_{max}$

In case of large $\kappa_{max}$ dissipation dominates the dynamics, such that
the system reaches a TISS, as shown in Fig.5 a) or a thermal phase, as shown
in Fig. 5 b).

Figure 5: Stroboscopic dynamics for $\kappa_{max}=10$, $\epsilon=0.02$ , a)
$m=\kappa_{0}$ and b) $m=\frac{\kappa_{0}}{4}$.

## References

  * [1] Anatoli Polkovnikov.  Universal adiabatic dynamics in the vicinity of a quantum critical point.  Phys. Rev. B, 72:161201, Oct 2005.
  * [2] Wojciech H. Zurek, Uwe Dorner, and Peter Zoller.  Dynamics of a quantum phase transition.  Phys. Rev. Lett., 95:105701, Sep 2005.
  * [3] Analabha Roy and Arnab Das.  Fate of dynamical many-body localization in the presence of disorder.  Phys. Rev. B, 91:121106, Mar 2015.
  * [4] Rozhin Yousefjani, Sougato Bose, and Abolfazl Bayat.  Floquet-induced localization in long-range many-body systems.  Phys. Rev. Res., 5:013094, Feb 2023.
  * [5] Krzysztof Sacha.  Modeling spontaneous breaking of time-translation symmetry.  Phys. Rev. A, 91:033617, Mar 2015.
  * [6] Michael P. Zaletel, Mikhail Lukin, Christopher Monroe, Chetan Nayak, Frank Wilczek, and Norman Y. Yao.  Colloquium: Quantum and classical discrete time crystals.  Rev. Mod. Phys., 95:031001, Jul 2023.
  * [7] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe.  Observation of a discrete time crystal.  Nature, 543(7644):217–220, Mar 2017.
  * [8] Andrew D. King, Sei Suzuki, Jack Raymond, Alex Zucca, Trevor Lanting, Fabio Altomare, Andrew J. Berkley, Sara Ejtemaee, Emile Hoskinson, Shuiyuan Huang, Eric Ladizinsky, Allison J. R. MacDonald, Gaelen Marsden, Travis Oh, Gabriel Poulin-Lamarre, Mauricio Reis, Chris Rich, Yuki Sato, Jed D. Whittaker, Jason Yao, Richard Harris, Daniel A. Lidar, Hidetoshi Nishimori, and Mohammad H. Amin.  Coherent quantum annealing in a programmable 2,000 qubit ising chain.  Nature Physics, 18(11):1324–1328, Nov 2022.
  * [9] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte, and R. Fazio.  Boundary time crystals.  Phys. Rev. Lett., 121:035301, Jul 2018.
  * [10] Victor Montenegro, Marco G. Genoni, Abolfazl Bayat, and Matteo G. A. Paris.  Quantum metrology with boundary time crystals.  Communications Physics, 6(1):304, Oct 2023.
  * [11] Phatthamon Kongkhambut, Jim Skulte, Ludwig Mathey, Jayson G. Cosme, Andreas Hemmerich, and Hans Keßler.  Observation of a continuous time crystal.  Science, 377(6606):670–673, August 2022.
  * [12] Midhun Krishna, Parvinder Solanki, Michal Hajdušek, and Sai Vinjanampathy.  Measurement-induced continuous time crystals.  Phys. Rev. Lett., 130:150401, Apr 2023.
  * [13] Tongjun Liu, Jun-Yu Ou, Kevin F MacDonald, and Nikolay I Zheludev.  Photonic metamaterial analogue of a continuous time crystal.  Nature Physics, 19(7):986–991, 2023.
  * [14] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath.  Discrete time crystals: Rigidity, criticality, and realizations.  Phys. Rev. Lett., 118:030401, Jan 2017.
  * [15] Michael P. Zaletel, Mikhail Lukin, Christopher Monroe, Chetan Nayak, Frank Wilczek, and Norman Y. Yao.  Colloquium: Quantum and classical discrete time crystals.  Rev. Mod. Phys., 95:031001, Jul 2023.
  * [16] F. M. Gambetta, F. Carollo, M. Marcuzzi, J. P. Garrahan, and I. Lesanovsky.  Discrete time crystals in the absence of manifest symmetries or disorder in open quantum systems.  Phys. Rev. Lett., 122:015701, Jan 2019.
  * [17] Jiehang Zhang, Paul W Hess, A Kyprianidis, Petra Becker, A Lee, J Smith, Gaetano Pagano, I-D Potirniche, Andrew C Potter, Ashvin Vishwanath, et al.  Observation of a discrete time crystal.  Nature, 543(7644):217–220, 2017.
  * [18] Soonwon Choi, Joonhee Choi, Renate Landig, Georg Kucsko, Hengyun Zhou, Junichi Isoya, Fedor Jelezko, Shinobu Onoda, Hitoshi Sumiya, Vedika Khemani, et al.  Observation of discrete time-crystalline order in a disordered dipolar many-body system.  Nature, 543(7644):221–225, 2017.
  * [19] Xiao Mi, Matteo Ippoliti, Chris Quintana, Ami Greene, Zijun Chen, Jonathan Gross, Frank Arute, Kunal Arya, Juan Atalaya, Ryan Babbush, et al.  Time-crystalline eigenstate order on a quantum processor.  Nature, 601(7894):531–536, 2022.
  * [20] Philipp Frey and Stephan Rachel.  Realization of a discrete time crystal on 57 qubits of a quantum computer.  Science advances, 8(9):eabm7652, 2022.
  * [21] C. W. von Keyserlingk, Vedika Khemani, and S. L. Sondhi.  Absolute stability and spatiotemporal long-range order in floquet systems.  Phys. Rev. B, 94:085112, Aug 2016.
  * [22] Norman Y. Yao and Chetan Nayak.  Time crystals in periodically driven systems.  Physics Today, 71(9):40–47, 09 2018.
  * [23] Dominic V. Else, Christopher Monroe, Chetan Nayak, and Norman Y. Yao.  Discrete time crystals.  Annual Review of Condensed Matter Physics, 11(Volume 11, 2020):467–499, 2020.
  * [24] Biao Huang, Ying-Hai Wu, and W. Vincent Liu.  Clean floquet time crystals: Models and realizations in cold atoms.  Phys. Rev. Lett., 120:110603, Mar 2018.
  * [25] Angelo Russomanno, Fernando Iemini, Marcello Dalmonte, and Rosario Fazio.  Floquet time crystal in the lipkin-meshkov-glick model.  Phys. Rev. B, 95:214307, Jun 2017.
  * [26] Achilleas Lazarides and Roderich Moessner.  Fate of a discrete time crystal in an open system.  Phys. Rev. B, 95:195135, May 2017.
  * [27] Bihui Zhu, Jamir Marino, Norman Y Yao, Mikhail D Lukin, and Eugene A Demler.  Dicke time crystals in driven-dissipative quantum many-body systems.  New Journal of Physics, 21(7):073028, jul 2019.
  * [28] Achilleas Lazarides, Sthitadhi Roy, Francesco Piazza, and Roderich Moessner.  Time crystallinity in dissipative floquet systems.  Phys. Rev. Res., 2:022002, Apr 2020.
  * [29] Federico Carollo, Kay Brandner, and Igor Lesanovsky.  Nonequilibrium many-body quantum engine driven by time-translation symmetry breaking.  Phys. Rev. Lett., 125:240602, Dec 2020.
  * [30] Zongping Gong, Ryusuke Hamazaki, and Masahito Ueda.  Discrete time-crystalline order in cavity and circuit qed systems.  Phys. Rev. Lett., 120:040404, Jan 2018.
  * [31] Heinz-Peter Breuer, Elsi-Mari Laine, and Jyrki Piilo.  Measure for the degree of non-markovian behavior of quantum processes in open systems.  Phys. Rev. Lett., 103:210401, Nov 2009.
  * [32] Dariusz Chruściński and Andrzej Kossakowski.  Non-markovian quantum dynamics: Local versus nonlocal.  Phys. Rev. Lett., 104:070406, Feb 2010.
  * [33] Ángel Rivas, Susana F Huelga, and Martin B Plenio.  Quantum non-markovianity: characterization, quantification and detection.  Reports on Progress in Physics, 77(9):094001, aug 2014.
  * [34] Heinz-Peter Breuer, Elsi-Mari Laine, Jyrki Piilo, and Bassano Vacchini.  Colloquium: Non-markovian dynamics in open quantum systems.  Rev. Mod. Phys., 88:021002, Apr 2016.
  * [35] V Mukherjee et al.  Efﬁciency of quantum controlled non-markovian thermalization.  New J. Phys., 2015:240602, June 2015.
  * [36] R. H. Dicke.  Coherence in spontaneous radiation processes.  Phys. Rev., 93:99–110, Jan 1954.
  * [37] Klaus Hepp and Elliott H Lieb.  On the superradiant phase transition for molecules in a quantized radiation field: the dicke maser model.  Annals of Physics, 76(2):360–404, 1973.
  * [38] F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael.  Proposed realization of the dicke-model quantum phase transition in an optical cavity qed system.  Phys. Rev. A, 75:013804, Jan 2007.
  * [39] Bogna Bylicka, Markus Johansson, and Antonio Acín.  Constructive method for detecting the information backflow of non-markovian dynamics.  Phys. Rev. Lett., 118:120501, Mar 2017.
  * [40] Yun-Yi Hsieh, Zheng-Yao Su, and Hsi-Sheng Goan.  Non-markovianity, information backflow, and system-environment correlation for open-quantum-system processes.  Phys. Rev. A, 100:012120, Jul 2019.
  * [41] Hossein Taheri, Andrey B. Matsko, Lute Maleki, and Krzysztof Sacha.  All-optical dissipative discrete time crystals.  Nature Communications, 13(1):848, Feb 2022.
  * [42] M. Tahir Naseem, Avijit Misra, Özgür E. Müstecaplioğlu, and Gershon Kurizki.  Minimal quantum heat manager boosted by bath spectral filtering.  Phys. Rev. Res., 2:033285, Aug 2020.


# Resolvent Analysis of an Under-expanded Planar Supersonic Impinging Jet

Qiong Liu111Present address: Department of Aerospace Engineering, UIUC.
<EMAIL_ADDRESS>The Ohio State University, OH 43210, USA  Chitrarth
Prasad222Postdoctoral Researcher, Department of Mechanical and Aerospace
Engineering and Datta V. Gaitonde333John Glenn Chair Professor, Department of
Mechanical and Aerospace Engineering The Ohio State University, OH 43210, USA

###### Abstract

This investigation aims to assess the effect of different types of actuator
forcing on the feedback loop of an under-expanded Mach 1.27 planar impinging
jet using a resolvent framework. To this end, we employ a Large Eddy
Simulation database as a truth model. The time and spanwise-averaged mean flow
is taken as an input to global stability and resolvent analyses with the
purpose of examining both the intrinsic instability and input-output
characteristics. The results show that the inherent instability and primary
energy amplification are attributed to the Kelvin-Helmholtz (K-H) instability.
Moreover, the K-H response modes obtained from the resolvent analysis are in
reasonable agreement with Spectral Proper Orthogonal Decomposition (SPOD)
modes from the unsteady LES data. Insights into noise control are obtained by
localizing the actuator forcing to the nozzle lip and the ground plate by
imposing component-wise forcing to mimic different notional actuators. It is
observed that energy amplification obtained for the localized component-wise
forcing is different from the global resolvent analysis and dependent on the
type of actuator. This provides insights into the type, wavenumber and
frequency of actuators for active flow control.

## Nomenclature

$\rho$ = | Density
---|---
$u$ = | Streamwise velocity
$T$ = | Temperature
$a_{\infty}$ = | speed of sound
$c$ = | Local speed of sound
$D$ = | Distance between nozzle walls
$\bm{B}$ = | Hydrodynamic FT Mode
$\psi^{\prime}$ = | Irrotational Scalar Potential
$\psi_{a}^{\prime}$ = | Acoustic Scalar Potential
$\psi_{t}^{\prime}$ = | Thermal Scalar Potential
$S$ = | Entropy
$St$ = | Oscillation frequency
$\beta$ = | Spanwise wavenumber
$\mathbb{R}$ = | Resolvent operator
$\mathbb{Q}=[\hat{\bm{q}}_{1},\hat{\bm{q}}_{2},\dots,\hat{\bm{q}}_{n}]$ = | response modes
$\mathbb{F}=[\hat{\bm{f}}_{1},\hat{\bm{f}}_{2},\dots,\hat{\bm{f}}_{n}]$ = | forcing modes
Subscripts
$j$ = | Jet exit quantities
$\infty$ = | Ambient quantities

##  1 Introduction

The study of a supersonic jet impinging on a flat plate is of great practical
relevance for the design of V/STOL aircraft and launch vehicles. The
interaction of the jet with the impingement surface produces a very complex
flow-field, which can be divided into a free jet region, an impingement zone
and a wall jet region. In addition to the acoustic radiation from the free-jet
and its subsequent reflection from the impingement surface, the extra noise
sources due to the complex flow-field at the impinging zone and the ensuing
wall jet often result in the flow-field being dominated by a self-reinforcing
aeroacoustic feedback loop [1]. This feedback process has been the subject of
several investigations, a comprehensive review of which can be found in Ref.
[2]. The impingement of the jet on the plate produces upstream travelling
acoustic waves that propagate to the nozzle exit and provide a periodic
forcing to the convectively unstable thin shear layer. These perturbations
grow into large scale structures that convect downstream and impinge on the
surface resulting in acoustic waves at a suitable phase and frequency which
complete the feedback loop [3]. The acoustic tones produced by this feedback
phenomenon result in significant fluctuating loads that can be up to 10dB
louder than their free-jet counterparts [4].

Both active and passive control techniques have been investigated to target
the individual elements of this feedback loop. Passive control techniques
focus on modifying the impingement geometry itself so as to avoid these
acoustic waves from perturbing the shear layer. These designs include, but are
not limited to, the use of oblique impingement [5, 6], the installation of
front-covers within the vertical flame path [7] and the introduction of a wall
with a circular cut-out to block the upstream traveling waves from the
impingement area towards the nozzle lip [8]. Active noise reduction methods on
the other hand, often focus on fluid injection at certain strategic locations
around the nozzle exit and beneath the impingement plate [9, 10, 11, 12].
Although these approaches have enjoyed some measure of success, a more
systematic study of the sensitivity of an impinging jet flow-field to
different types of actuators forcing is lacking.

Active control techniques are more attractive than their passive counterparts,
as they can be operated on-demand. The most general parameters in the design
of any active noise control methodology is the type and spatial distribution
(wavenumber) of the actuators, and the amplitude and frequency of actuation.
The goal of the present study is assess the effect of choosing different types
of actuator forcing and location on the feedback loop of an under-expanded
Mach 1.27 planar impinging jet using a resolvent framework. This is motivated
by the promising success of resolvent analysis in controlling the feedback
process in cavity flows [13, 14]. Although the mean flow needed for resolvent
analysis could be obtained from a Reynolds Averaged Navier-Stokes (RANS)
solution; we employ a carefully computed Mach 1.27 LES database to serve as a
truth model and therefore use the same database to obtain the mean flow for
analysis purposes. A detailed description of the LES database is provided in
Section 2. Section 3 provides a brief discussion of the resolvent framework
used in this study. Some preliminary results are presented in Section 4. A
summary of ongoing work and concluding remarks are given in Section 5.

##  2 LES Database

The LES database used in the present study consists of an unheated planar
under-expanded Mach 1.27 jet impinging on a flat plate placed at 4 jet
diameters from the nozzle exit. This impingement height is chosen based on
observed resonant mode sensitivity to nearby heights in previous single
impinging jet studies [15, 16, 17]. The database is generated by solving the
full 3D compressible Navier-Stokes equations in generalized curvilinear
coordinates using a well-validated inhouse solver SAFF. A third-order upwind
biased Roe scheme [18] in conjunction with the van Leer harmonic limiter [19]
is employed for the spatial discretization of the inviscid fluxes. The viscous
fluxes are computed using a second-order central differencing scheme, whereas
the time integration is achieved by using an implicit second-order
diagonalized [20] Beam Warming scheme [21].

A schematic of the computational domain is presented in Fig. 1.

Figure 1: A schematic of the computational domain for LES calculations. The
geometry consists of two parallel walls at distance of 1 inch with a sonic
inflow.

The geometry consists of two parallel walls separated at a distance $D$ of 1
inch covering the entire spanwise extent of the domain. A cartesian coordinate
system is chosen. with the jet exhausting in the $x-$direction. The nozzle
exit is located at $x=0.5D$, whereas the plate is located at $x=4.5D$. The
domain extends up to $15D$ in the $y-$direction on either side of the jet
axis, whereas the spanwise extent of the domain ($z-$direction) is $10D$. The
calculations are performed in a structured mesh containing 521, 917 and 81
points in the $x$, $y$ and $z$ directions respectively. Non-reflecting
freestream boundary conditions are applied at all the outer boundaries. Iso-
thermal wall conditions with $T_{w}=1.12T_{\infty}$ are applied at all solid
boundaries including the nozzle walls and the impingement plate. The inflow
consists of a laminar boundary layer profile with thickness $\delta=0.1D$
following the approach described by Bogey and Bailly [22] to match the onset
of shear layer turbulence observed in experiments [23]. This manner of inflow
specification has proven effective in replicating the experimentally observed
shear layer development in both single and twin impinging jets [24].

All flow quantities are used in their dimensionless form. The jet exit
conditions, $U_{j}=343.7\mathrm{m/s},T_{j}=293\mathrm{K}\text{ and
}\rho_{j}=1.2061\mathrm{kg/m^{3}}$ are used as the reference quantities for
velocity, temperature and density respectively. The nozzle exit diameter,
$D=0.0254m$ is used as the reference length, whereas the pressure is
normalized by $\rho_{j}{U_{j}}^{2}$. The ambient temperature is
$T_{\infty}=293\mathrm{K}$.

Figure 2: Density gradient magnitude contours with background dilatation
contours in a 2D streamwise plane.

Figure 2 shows the density gradient magnitude contours (color) superimposed on
velocity dilatation contours (gray-scale) at an arbitrary time step in a 2D
streamwise plane. The density gradient contours illustrate the turbulent
structures and the shock cells in the jet plume, whereas the background
dilatation contours give a qualitative picture of the acoustic waves generated
due to this mixing. Two distinct types of sound waves can be identified: the
broadband shock associated noise in the sideline and upstream direction at
locations where the shock cells interact with the shear layer and the upstream
travelling waves due to the impinging zone and the ensuing supersonic wall
jet.

Figure 3 shows the power spectral density of pressure as a function of
Strouhal number ($St=fU_{j}/D$) at the point marked by yellow circle in Fig.
2. The high-amplitude discrete frequency feedback tones due to the feedback
loop between the nozzle exit and the ground plate are clearly visible in the
spectrum. In order to validate these impinging tone frequencies, we make use
of Powell’s feedback model. Following Powell [25], the fundamental impingement
tone frequency $F$, and its harmonics $n$ are commonly expressed as

$\frac{n}{F}=\frac{H}{U}+\frac{H}{a}+p,$ (1)

where $H$ represents the height of impingement, $U$ is the average shear layer
convection speed, and $a$ is the speed of sound. The time lag in the
receptivity and sound generation components is accounted for in the phase lag
term p in order to recover the observed acoustic tones. In the present
investigation, the speed of the downstream convecting coherent structures and
upstream propagating acoustic waves are simultaneously measured by a line of
simulated probes along the jet shear layer from the nozzle lip to the
impingement plate (not shown). The pressure fluctuations at every point along
the shear-layer line of probes is cross-correlated with that at the nozzle
exit location to calculate the shear layer convection speed. The speed of
sound is assumed to be $343.7$ m/s, resulting in an equivalent fundamental
frequency equal to $St=0.10$. This approach sets the phase lag term $p$ in
equation 1 to zero by assuming it is accounted for in the computed average
shear-layer convection speed. Using this key result the impinging tones
observed in Fig. 3 can be identified as harmonics $n=2$ to $6$ of the
fundamental frequency. This provides confidence in the use of this database in
the present work. A detailed calculation of the individual elements of the
feedback loop will be presented in the final manuscript.

Figure 3: Power Spectral Density of near-field pressure. The impinging tone
harmonics $n=2$ to $6$ are identified.

The LES flow-field is averaged in both time and spanwise direction over
$300D/U_{j}$ time snapshots to provide an input for the global stability and
resolvent analyses. A brief summary of these techniques is presented next.

##  3 Modal analysis

###  3.1 Global stability analysis

An adoption of Reynolds averaging in the turbulent flow separates the mean
Reynolds stresses and fluctuating Reynolds stress. The mean Reynolds stresses
are implicitly included in the turbulent mean flow that are obtained from the
LES. In comparison, the fluctuating Reynolds stresses are amplified/diminished
at certain frequencies, which is assumed as a nonlinear forcing acting on the
large-scale structures in the turbulent flow. This point of view is similar to
the one used in Refs [26, 27] for modal analyses of high Reynolds number
turbulent flows and has been shown to provide significant insights in
understanding flow dynamics and informing flow control. In the present study,
we adopt this methodology by first neglecting the nonlinear forcing and
assessing the unforced behavior of the linear modes based on a time averaged
impinging jet flow field. This information is then utilized to assess the
response of the flow to non-linear forcing as shown later.

The instantaneous flow variables are Reynolds decomposed into a time-averaged
flow state $\bar{\bm{q}}\equiv[\bar{\rho},\bar{u},\bar{v},\bar{w},\bar{T}]$
and a statistically stationary fluctuating components
$\bm{q}^{\prime}\equiv[\rho^{\prime},u^{\prime},v^{\prime},w^{\prime},T^{\prime}]$.
By substituting the Reynolds decomposed state variables into the Navier–Stokes
(NS) equations, the governing equation for the fluctuation $\bm{q}^{\prime}$
becomes

$\frac{\partial\bm{q}^{\prime}}{\partial
t}=\mathbb{L}(\bar{\bm{q}})\bm{q}^{\prime},$ (2)

where $\mathbb{L}(\bar{\bm{q}})$ is the Navier-Stokes operator linearized
about the time-averaged flow state $\bar{\bm{q}}$.

Since the present planar impinging flow is periodic in the spanwise direction,
the fluctuation term $\bm{q}^{\prime}$ can be expressed as Fourier modes with
real a spanwise wavenumber $\beta$ and a complex frequency $\lambda$ as
follows

$\bm{q}^{\prime}(x,y,z,t)={\hat{\bm{q}}}(x,y)e^{\text{-i}(\beta z+\lambda
t)}.$ (3)

Here the spanwise wavenumber is normalized by $D$ and the complex value of
$\lambda$ is non-dimensionalized using $D/2\pi U_{j}$. The real part
$\lambda_{r}$ represents the oscillation frequencies and imaginary part
$\lambda$ denotes the growth ($\lambda>0$) or damping rates ($\lambda<0$) of
linear modes. Substituting the modal expression of $\bm{q}^{\prime}(x,y,z,t)$
from eqn. 3 into eqn. 2 yields an eigenvalue problem of the linear system

$\mathbb{L}(\bar{\bm{q}},\beta)\hat{\bm{q}}=\lambda\hat{\bm{q}}.$ (4)

Given the time and spanwise averaged flow with prescribed spanwise
wavenumbers, this large eigenvalue problem can be solved to obtain linear
modes in terms of their eigenvalues and eigenvectors $(\lambda,\hat{\bm{q}})$.
$\mathbb{L}$ represents the system of the linearized compressible continuity,
Navier-Stokes, and energy equations for an ideal gas in Cartesian coordinates.

For computational efficiency, we construct the linear operator $\mathbb{L}$
using a smaller computational domain and reducing the grid relative to the one
used in LES. This computational domain extents from $-9D\leq y\leq 9D$ and is
sufficient to illustrate the primary flow physics of interest that occur
between the jet nozzle and impingement surface. At inflow, Dirichlet boundary
condition is used on all perturbation variables. At the outflow, sponge zones
with width of $3D$ are prescribed. The spatial discretion uses spectral
element methods as described in Refs. [28, 29] with element size of $1327$ and
polynomial order of $5$. The eigenvalue problem is solved using an efficient
shift-and-invert Arnoldi algorithm with six solved eigensolution at Krylov
space dimension of 16 and a residual tolerance of $10^{-7}$ at each shift
value. The results converge to at least seven significant figures and are
verified to be independent of the domain size and mesh resolution.

###  3.2 Resolvent analysis

The input-output characteristics of the impinging jet flow are obtained from a
resolvent analysis. We now retain the effect of the nonlinear forcing and
consider its effect on the flow through the resolvent operator. The resulting
fluctuation NS equations can be expressed as an input-output system [30, 27,
31] given by,

$\frac{\partial\bm{q}^{\prime}}{\partial
t}=\mathbb{L}(\bar{\bm{q}})\bm{q}^{\prime}+\mathbb{M}\bm{f}^{\prime}$ (5)

where $\mathbb{L}(\bar{q})$ is the Navier-Stokes operator linearized about the
base state $\bar{\bm{q}}$, the finite-amplitude nonlinear terms are
incorporated in ${\bm{f}}^{\prime}$, and $\mathbb{M}$ is a coupling matrix and
is discussed later.

We use of a similar modal ansatz with spanwise wavenumber $\beta$ and temporal
frequency $\lambda$ for both the state vector fluctuations and the forcing as
shown below

$\begin{array}[]{c}\bm{q}^{\prime}(x,y,z,t)=\hat{\bm{q}}(x,y)e^{-\text{i}(\lambda
t+\beta z)},\\\ \bm{f}^{\prime}(x,y,z,t)=\hat{\bm{f}}(x,y)e^{-\text{i}(\lambda
t+\beta z)}.\end{array}$ (6)

Inserting these into eqn. 5, we obtain

$\hat{{\bm{q}}}=[-\text{i}\lambda\mathbb{I}-\mathbb{L}({\bm{\bar{q}}},\beta)]^{-1}\mathbb{M}\hat{{\bm{f}}},$
(7)

where
$\mathbb{R}=[\text{i}\lambda\mathbb{I}-\mathbb{L}({\bm{\bar{q}}},\beta)]^{-1}\mathbb{M}$
is the resolvent operator, which serves as a transfer function between the
input $\hat{\bm{f}}$ and the corresponding output $\hat{\bm{q}}$ for a given
flow state ($\bm{\bar{q}}$) and modal parameters ($\beta$ and $\lambda$).

The energy amplification of the system is evaluated from the ratio of output
to input energy $\frac{||\hat{\bm{q}}||_{E}}{||\hat{\bm{f}}||_{E}}$, where
$||\cdot||_{E}$ is an energy norm. A proper selection of the energy norm is
crucial in order to prevent divergence in non-modal results. Here we use the
compressible energy _Chu_ -norm [32] given by

$E=\int_{S}\left[\frac{\bar{a}^{2}\rho^{2}}{\gamma\bar{\rho}}+\bar{\rho}({u}^{2}+{v}^{2}+{w}^{2})+\frac{\bar{\rho}C_{v}T^{2}}{\bar{T}}\right]{\text{d}}s$

where $S$ is the domain of interest. A singular value decomposition (SVD) of
the resolvent operator facilitates a ranking of the energy amplification ratio
in descending order. This yields

$W^{\frac{1}{2}}\mathbb{R}W^{-\frac{1}{2}}=\mathbb{Q}\Sigma\mathbb{F}^{\ast},$
(8)

in terms of the weight matrix, $W$, based on the compressible energy norm $E$.
The matrix
$\mathbb{Q}=[\hat{\bm{q}}_{1},\hat{\bm{q}}_{2},\dots,\hat{\bm{q}}_{n}]$ holds
the set of optimal response directions and
$\mathbb{F}=[\hat{\bm{f}}_{1},\hat{\bm{f}}_{2},\dots,\hat{\bm{f}}_{n}]$
contains the corresponding forcing directions, where
$\hat{\bm{q}}_{i}=(\hat{\rho}_{r},\hat{u}_{r},\hat{v}_{r},\hat{w}_{r},\hat{T}_{r},)$
and
$\hat{\bm{f}}_{i}=(\hat{\rho}_{f},\hat{u}_{f},\hat{v}_{f},\hat{w}_{f},\hat{T}_{f},)$
with $n$ is number of solved singular values. The superscript $\ast$ denotes
the Hermitian transpose. The singular values
$\Sigma=\text{diag}(\sigma_{1},\sigma_{2},\dots,\sigma_{n})$ represent the
energy amplification (gain) between response and forcing modes.

##  4 Results and discussions

###  4.1 Global linear modes

Figure 4: Eigen-spectra of the time and spanwise-mean base flow at $\beta=0$,
$\pi$ and 2$\pi$. The blue and black circles are used to indicate the
eigenvalues of global mode shown in figure 5.  Figure 5: Real component of
transverse velocity for global modes at $\beta=0$, $\pi$ and $2\pi$.

The time and spanwise averaged impinging jet flow is unstable. We perform
linear stability analysis of time-averaged flow from $\beta=0$ to
$\beta=10\pi$ as described in Section 3.1. For each spanwise wavenumber
$\beta$, we shift the search region sixty positions along with the frequency
range of $0\leq St\leq 1.25$ and seven values along with the growth rate range
of $-0.06\leq\lambda_{r}/2\pi U_{r}ef\leq 0.14$ to cover the eigenvalues of
interest. Figure 4 shows the eigenspectra at $\beta=0$, $\pi$ and $2\pi$. The
eigenvalues have positive growth rates, which reveal that the spanwise and
time-averaged mean flow is unstable. Figure 5 shows the modal structure of
these unstable modes at $St=0.41$ and $0.73$ as a function of $\beta$ using
the real component of the transverse velocity $\hat{v}$. These modes present
distinctively strong structures along the shear layer in the free jet region.
The structures of the modes gradually vanish at the impinging region as well
as the wall jet region. This is unsurprising as these modes resemble the
Kelvin-Helmholtz instability, which is well known to be the intrinsic
instability in the impinging jet flow. As the frequency increases, the scale
of the mode structures changes reversely. At a similar frequency range, the
global modes exhibit symmetric or antisymmetric structures with different
growth rates. For example, at $\beta=0$ and frequency $St\sim 0.41$, the K-H
instability is symmetric for $\lambda_{r}D/2\pi a_{\infty}=0.1$, while it is
anti-symmetric for $\lambda_{r}D/2\pi a_{\infty}=0.08$.

The existence of both symmetric and anti-symmetric modes however, is not
surprising as the previous studies on circular jets [25] have shown that the
shear layer can support four resonant modes labeled as axisymmetric mode A
(divided into A1 and A2 ), flapping mode B, helical mode C, and sinuous mode
D. For the present case, the instantaneous LES flow field shows a strong
flapping motion in the jet plume at $H/D=4$, indicating the dominant role of
the anti-symmetric modes. This dominance however, is likely to vary with the
impingement height as seen previously by Davis _et al._ [16], where a change
in height from $H/D=4$ to $4.5$ resulted in different mode shapes. Similar
analysis for different heights is being pursued by the authors and will be
presented in the final manuscript.

###  4.2 Discounted resolvent analysis

Since the time-averaged flow is temporally unstable, it becomes crucial to
highlight amplifications that occur on a shorter time scale than those
associated with the asymptotic behavior observed with classical instability
theory. This consideration aids in achieving the main objective of finding
preferred energy transfer mechanisms from the mean flow to the fluctuation
field, which are necessary to provide physical insight into potential flow
control strategies. To achieve this objective, the discounting technique [33,
34, 14] is employed to obtain forcing and response modes. The method
introduces a free parameter, denoted the discounting parameter $\kappa$, the
choice of which is predicated on information about the most unstable growth
rate as obtained from the stability analysis. Other methods to address the
unstable linear operator may be found in Ref. [35]. Consequently, we
substitute the resolvent operator $\mathbb{R}$ by discounted resolvent
operator $\mathbb{R^{\prime}}$ in eqn.7 as follows

$\mathbb{R^{\prime}}=[\text{i}(\lambda+\text{i}\kappa)\mathbb{I}-\mathbb{L}({\bm{\bar{q}}},\beta)]^{-1}\mathbb{M}.$
(9)

For the present test case, the discounting parameter is chosen as
$\kappa=0.159$ based on the global stability analysis results.

Figure 6: The leading energy amplification $\sigma_{1}$ over $\beta$ and $St$
for a discount parameter of $\kappa=0.14$

Figure 6 shows the leading singular value as a function of $\beta$ and $St$,
revealing the optimal energy amplification to harmonic forcing. A large
amplification around $0.2<St<0.6$ and $\pi<\beta<3\pi$ is observed from the
figure. In addition, since the shear layer modes primarily stem from two-
dimensional shear-layer instabilities, the gain distribution with $\beta=0$
exhibits discrete peaks at the impingement tones.

Figure 7: Leading response mode at $\beta=\pi$ at frequency $St=0.4$, $0.6$,
$0.75$ and $0.95$, illustrating with the real component of $\hat{v}$.

The representative response modes are presented in Fig. 7 for $\beta=\pi$ at
$St=0.4$, 0.6, 0.75 and 0.95 using the real transverse velocity component.
Similar to the previous stability analysis, these modes resemble a Kelvin-
Helmholtz type of mode structure. As expected, with an increase in frequency,
the mode structures become smaller and exhibit a smaller wavelength. Both
symmetric and anti-symmetric response modes are observed. As the spanwise
wavenumber increases, the response mode structures are confined to the shear-
layer region with smaller-scale structures (not shown).

In order to assess the accuracy by which resolvent modes capture the mode
shapes, we utilize Spectral Proper Orthogonal Decomposition (SPOD) of the
unsteady LES data. SPOD is a data-driven technique that decomposes the flow-
field into a set of orthogonal modes in the frequency domain that optimally
capture the flow’s energy based on a user-defined norm. It is well known that
SPOD modes are identical to the resolvent modes when the resolvent mode
expansion coefficients are uncorrelated, i.e., white-noise forcing [36].

SPOD modes aredetermined by solving the following Fredholm integral eigenvalue
problem:

$\int
R\left(\textbf{x},\textbf{x'},\textit{St}\right)W(\textbf{x'})\phi_{n}\left(\textbf{x'}\right)\text{d{x'}}=\lambda_{n}\left(\textit{St}\right)\phi_{n}\left(\textbf{x},\textit{St}\right).$
(10)

Here $\phi_{n}\left(\textbf{x},\textit{St}\right)$ represents the
$n^{\text{th}}$ SPOD mode at a particular _St_ value and $\lambda_{n}$ is the
energy associated with that mode. $W$ and $R$ represent the weight and the
cross-spectral density matrices respectively. In order to perform the Fourier
transform, the flow-field is arranged in 46 blocks of 256 snapshots with a
$75\%$ overlap. A periodic Hann window is employed to minimize spectral
leakage. Special care is taken to ensure that the conclusions drawn from the
SPOD modes are insensitive to these choices.

Figure 8 shows the leading pressure SPOD modes along the jet center plane at
$St=0.4$ and $0.6$. The mode shapes exhibit a symmetric and anti-symmetric
behaviour similar to the resolvent modes at these frequencies as shown
previously in Fig. 7 for transverse velocity. A more comprehensive one-to-one
comparison will be presented in the final manuscript. Nevertheless, this
demonstrates the capability of resolvent analysis to accurately represent the
flow’s response to harmonic forcing. This is leveraged next to provide
insights into flow control.

(a)

(b)

Figure 8: Leading SPOD pressure modes at (a) $St=0.4$ and (b) $St=0.6$.

###  4.3 Flow response to localized component-wise forcing

A control strategy may be designed based on the information of the response
output and forcing input. The primary forcing modes (not shown) exhibit white
noise forcing in the shear layer near the nozzle lip. As such, forcing around
the nozzle lip is a natural choice for flow control, with the potential to
significantly modify the behavior of shear layer structures in the free jet
region. Inspired by the location control actuation, we consider localization
of the input forcing to two locations: (i) the nozzle walls near the nozzle
exit $(x/D\leq 0.5)$ and (ii) impinging ground surface $(x/D>4.1)$.

The component-wise forcing input is selected by prescribing constraints in the
coupling matrix $\mathbb{M}$ in eqn. 7. The forcing contains five components
expressed as
$\hat{\bm{f}}=[\hat{\rho}_{f},\hat{u}_{f},\hat{v}_{f},\hat{w}_{f},\hat{T}_{f}]$
in eqn. (2), of which the spanwise velocity component forcing $\hat{w}_{f}$ is
neglected. In the coupling matrix $\mathbb{M}$, the elements at the location
of interest are set to the chosen component forcing, which localizes the input
term $\mathbb{M}\hat{\bm{f}}$. The matrix $\mathbb{M}$ thus serves two
purposes, spatial restriction, and imposition of a component-wise forcing
filter. The results identify the influence of different types of forcing at
these two locations on energy amplification and response structures.

Figure 9 show the energy amplification results with the localized component-
wise forcings corresponding to location (i). The streamwise velocity forcing
results in the strongest energy amplification, followed by the $\hat{v}_{f}$ ,
$\hat{\rho}_{f}$ and $\hat{T}_{f}$ forcings in that order. The localization of
$\hat{u}_{f}$ forcing dramatically reduces the energy amplification to about
one-third that with the global forcing shown in Fig. 6. Nonetheless, compared
to the global forcing structures, the much smaller localized forcing at a
position where an actuator may actually be placed attains substantial energy
amplification. A slightly different trend of energy amplification is observed
for the three-dimensional cases ($\beta>0$) with another gain peak around
$0.6<St<0.8$.

Figure 9: Optimal gain of localized-componentwise forcing versus spanwise
wavenumbers and frequencies.

For location (ii), a strong energy amplification is detected in the streamwise
velocity forcing followed by the $\hat{\rho}_{f}$, $\hat{v}_{f}$, and
$\hat{T}_{f}$ forcings in that order. Contrary to the gain maps at location
(i), the gain obtained from a three-dimensional forcing around
$2\pi<\beta<4\pi$ and $0.2<St<0.6$ is significantly higher than gains obtained
at $\beta=0$. These gain contours provide insights into the type, wavenumber
and the frequency of actuators that can be employed for noise control. A more
complete analysis in underway and will be presented in the final manuscript.

##  5 Ongoing work and Conclusion

We perform a large-eddy simulation over a supersonic planar impinging jet at
$M=1.27$ and Reynolds number based on nozzle size of 527,000. The frequency
tones established between nozzle exit and ground surface agree with the
prediction from Powell’s formula.

We further adopt global stability and resolvent analysis to identify intrinsic
instability. From the global stability analysis, we learn that the time and
spanwise averaged flow is unstable. The mode structure indicates the Kelvin-
Helmholtz instability as the unique instability existing in the impinging jet
flow. Both symmetric and anti-symmetric modes are observed from the global
stability analysis. Resolvent analysis is performed over frequency $(St=fU/D)$
and spanwise wavenumber $\beta$ ranges of $0.1\leq St\leq 1.05$ and
$0\leq\beta\leq 10\pi$, respectively. The gain-peak resides in the $0.2\leq
St\leq 0.5$ and $\pi\leq\beta\leq 3\pi$ ranges. We observe that the energy
amplification mechanics is dominated by a Kelvin–Helmholtz type of response.
At $St\sim 0.4$, this KH response exhibits a symmetric structure at $\beta=0$,
$\pi$, $3\pi$, $4\pi$, $5\pi$ and $6\pi$ with higher energy amplification
compared to its anti-symmetric counterparts at other $\beta$-values. The
efforts focus on localizing the forcing near the nozzle to consider notional
actuators that can separately introduce individual velocity components,
density, and temperature forcing, respectively, to target noise mitigation.
The result showed a slightly different energy amplification by localized
component-wise forcing.

## Acknowledgments

This work was performed in part under the sponsorship of the Office of Naval
Research (Contract N00014-18-1-2506) with Dr. D. Gonzalez serving as Project
Monitor. The views and conclusions contained herein are those of the authors
and do not represent the opinion of the Office of Naval Research or the U.S.
government. Computational resources were provided by the DoD High Performance
Computing Modernization Program as well as the Ohio Supercomputer Center.
Several figures were made using FieldView software with licenses obtained from
the Intelligent Light University Partnership Program.

## References

  * Powell [1988] Powell, A., “The sound-producing oscillations of round underexpanded jets impinging on normal plates,” _The Journal of the Acoustical Society of America_ , Vol. 83, No. 2, 1988, pp. 515–533.
  * Edgington-Mitchell [2019] Edgington-Mitchell, D., “Aeroacoustic resonance and self-excitation in screeching and impinging supersonic jets–a review,” _International Journal of Aeroacoustics_ , Vol. 18, No. 2-3, 2019, pp. 118–188.
  * Prasad et al. [2021a] Prasad, C., Stahl, S., and Gaitonde, D., “Exchange Mechanisms between Hydrodynamic and Acoustic Components of an Under-expanded Supersonic Impinging Jet,” _AIAA AVIATION 2021 FORUM_ , 2021a, p. 2118.
  * Krothapalli et al. [1999] Krothapalli, A., Rajkuperan, E., Alvi, F., and Lourenco, L., “Flow field and noise characteristics of a supersonic impinging jet,” _4th AIAA/CEAS aeroacoustics conference_ , 1999, p. 2239.
  * Worden et al. [2013-2187] Worden, T., Gustavsson, J., Shih, C., and Alvi, F. S., “Acoustic measurements of high-temperature supersonic impinging jets in multiple configurations,” _19th AIAA/CEAS Aeroacoustics Conference_ , 2013-2187.
  * Brehm et al. [2016] Brehm, C., Housman, J. A., and Kiris, C. C., “Noise generation mechanisms for a supersonic jet impinging on an inclined plate,” _Journal of Fluid Mechanics_ , Vol. 797, 2016, pp. 802–850.
  * Tsutsumi et al. [2015] Tsutsumi, S., Ishii, T., Ui, K., Tokudome, S., and Wada, K., “Study on acoustic prediction and reduction of Epsilon launch vehicle at liftoff,” _Journal of Spacecraft and Rockets_ , 2015.
  * Prasad et al. [2021b] Prasad, C., Yenigelen, E., and Morris, P. J., “Effect of Launchpad Modification on the Hydrodynamic and Acoustic Modes of an Impinging Jet,” _AIAA Scitech 2021 Forum_ , 2021b, p. 1417.
  * Alvi et al. [2003] Alvi, F. S., Shih, C., Elavarasan, R., Garg, G., and Krothapalli, A., “Control of supersonic impinging jet flows using supersonic microjets,” _AIAA journal_ , Vol. 41, No. 7, 2003, pp. 1347–1355.
  * Norum [2004] Norum, T., “Reductions in multi-component jet noise by water injection,” _10th AIAA/CEAS Aeroacoustics Conference_ , 2004, p. 2976.
  * Fukuda et al. [2011-2814] Fukuda, K., Tsutsumi, S., Shimizu, T., Takaki, R., and Ui, K., “Examination of sound suppression by water injection at lift-off of launch vehicles,” _17th AIAA/CEAS Aeroacoustics Conference (32nd AIAA Aeroacoustics Conference)_ , 2011-2814.
  * Salehian et al. [2018-0519] Salehian, S., Kourbatski, K., Golubev, V. V., and Mankbadi, R. R., “Numerical Aspects of Rocket Lift-off Noise with Launch-Pad Aqueous Injection,” _2018 AIAA Aerospace Sciences Meeting_ , 2018-0519.
  * Liu et al. [2018] Liu, Q., Sun, Y., Cattafesta, L. N., Ukeiley, L. S., and Taira, K., “Resolvent Analysis of Compressible Flow over a Long Rectangular Cavity,” AIAA Paper 2018-0588, 2018.
  * Liu et al. [2021] Liu, Q., Sun, Y., Yeh, C.-A., Ukeiley, L. S., Cattafesta, L. N., and Taira, K., “Unsteady control of supersonic turbulent cavity flow based on resolvent analysis,” _Journal of Fluid Mechanics_ , Vol. 925, 2021.
  * Bhargav et al. [2021] Bhargav, V. N., Song, M., Sellappan, P., Alvi, F. S., and Kumar, R., “Experimental Characterization of Supersonic Single- and Dual-Impinging Jets,” _AIAA Journal_ , 2021, pp. 1–17.  10.2514/1.j059687.
  * Davis et al. [2015] Davis, T., Edstrand, A., Alvi, F., Cattafesta, L., Yorita, D., and Asai, K., “Investigation of impinging jet resonant modes using unsteady pressure-sensitive paint measurements,” _Experiments in Fluids_ , Vol. 56, No. 5, 2015, pp. 1–13.  10.1007/s00348-015-1976-9.
  * Kumar et al. [2013] Kumar, R., Wiley, A., Venkatakrishnan, L., and Alvi, F., “Role of coherent structures in supersonic impinging jets,” _Physics of Fluids_ , Vol. 25, No. 7, 2013.  10.1063/1.4811401.
  * Roe [1981] Roe, P. L., “Approximate Riemann solvers, parameter vectors, and difference schemes,” _Journal of computational physics_ , Vol. 43, No. 2, 1981, pp. 357–372.
  * Van Leer [1979] Van Leer, B., “Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method,” _Journal of computational Physics_ , Vol. 32, No. 1, 1979, pp. 101–136.
  * Pulliam and Chaussee [1981] Pulliam, T. H., and Chaussee, D., “A diagonal form of an implicit approximate-factorization algorithm,” _Journal of Computational Physics_ , Vol. 39, No. 2, 1981, pp. 347–363.
  * Beam and Warming [1978] Beam, R. M., and Warming, R., “An implicit factored scheme for the compressible Navier-Stokes equations,” _AIAA journal_ , Vol. 16, No. 4, 1978, pp. 393–402.
  * Bogey and Bailly [2010] Bogey, C., and Bailly, C., “Influence of nozzle-exit boundary-layer conditions on the flow and acoustic fields of initially laminar jets,” _Journal of Fluid Mechanics_ , Vol. 663, 2010, p. 507–538.
  * Nataraj Bhargav et al. [2020] Nataraj Bhargav, V., Song, M., Sellappan, P., Alvi, F. S., and Kumar, R., “Unsteady Characteristics of Resonant Supersonic Dual Impinging Jets,” _AIAA Scitech 2020 Forum_ , 2020.
  * Stahl et al. [2021] Stahl, S. L., Prasad, C., and Gaitonde, D. V., “Distinctions between single and twin impinging jet dynamics,” _The Journal of the Acoustical Society of America_ , Vol. 150, No. 2, 2021, pp. 734–744.
  * Powell [1953] Powell, A., “On edge tones and associated phenomena,” _Acustica_ , Vol. 3, 1953, pp. 233–43.
  * Farrell and Ioannou [1993] Farrell, B. F., and Ioannou, P. J., “Stochastic forcing of the linearized Navier–Stokes equations,” _Physics of Fluids A: Fluid Dynamics_ , Vol. 5, No. 11, 1993, pp. 2600–2609.
  * McKeon and Sharma [2010] McKeon, B. J., and Sharma, A. S., “A critical-layer framework for turbulent pipe flow,” _J. Fluid Mech._ , Vol. 658, 2010, pp. 336–382.
  * Liu et al. [2016] Liu, Q., Gómez, F., and Theofilis, V., “Linear instability analysis of low-${R}e$ incompressible flow over a long rectangular finite-span open cavity,” _Journal of Fluid Mechanics_ , Vol. 799, 2016.
  * Liu and Gaitonde [2021] Liu, Q., and Gaitonde, D., “Acoustic response of turbulent cavity flow using resolvent analysis,” _Physics of Fluids_ , Vol. 33, No. 5, 2021, p. 056102\.
  * Jovanović and Bamieh [2005] Jovanović, M. R., and Bamieh, B., “Componentwise energy amplification in channel flows,” _J. Fluid Mech._ , Vol. 534, 2005, pp. 145–183.
  * Schmid and Henningson [2012] Schmid, P. J., and Henningson, D. S., _Stability and transition in shear flows_ , Springer, 2012.
  * Chu [1965] Chu, B.-T., “On the energy transfer to small disturbances in fluid flow (Part I),” _Acta Mechanica_ , Vol. 1, No. 3, 1965, pp. 215–234.
  * Jovanović [2004] Jovanović, M. R., “Modeling, analysis, and control of spatially distributed systems,” Ph.D. thesis, University of California, Santa Barbara, 2004\.
  * Yeh and Taira [2019] Yeh, C.-A., and Taira, K., “Resolvent-analysis-based design of airfoil separation control,” _J. Fluid Mech._ , Vol. 867, 2019, pp. 572–610.
  * Pickering et al. [2019-2454] Pickering, E. M., Rigas, G., Sipp, D., Schmidt, O. T., and Colonius, T., “Eddy viscosity for resolvent-based jet noise models,” _AIAA Paper No._ , 2019-2454.
  * Towne et al. [2018] Towne, A., Schmidt, O. T., and Colonius, T., “Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis,” _Journal of Fluid Mechanics_ , Vol. 847, 2018, pp. 821–867.


# Explicit, time-reversible and symplectic integrator for Hamiltonians in
isotropic uniformly curved geometries

Ana Silva  Eitan Ben Av  Efi Efrati<EMAIL_ADDRESS>Department of
Physics of Complex Systems, Weizmann Institute of Science, Rehovot 76100,
Israel

###### Abstract

The kinetic term of the $N$-body Hamiltonian system defined on the surface of
the sphere is non-separable. As a result, standard explicit symplectic
integrators are inapplicable. We exploit an underlying hierarchy in the
structure of the kinetic term to construct an explicit time-reversible
symplectic scheme of second order. We use iterative applications of the method
to construct a fourth order scheme and demonstrate its efficiency.

##  I Introduction

Most Hamiltonian systems cannot be solved analytically and their long time
dynamics are thus commonly studied using numerical schemes that approximate
the solutions of Hamilton’s equations. Symplectic methods conserve the
symplectic structure of the evolution equations and are therefore considered
particularly suited and stable for numerically integrating autonomous
Hamiltonian systems Yoshida (1993); Gladman et al. (1991); Channell and Scovel
(1990); Reich (1996). Simulating the long time dynamics of such systems
requires not only stability and accuracy, but also efficiency, thus favoring
explicit symplectic methods.

Explicit symplectic methods are abundant for separable Hamiltonians where the
kinetic term depends solely on the momenta and the potential term depends
solely on the coordinates (Yoshida, 1993; Haruo, 1990; Forest and Ruth, 1990;
Leimkuhler and Reich, 2004a; Hairer et al., 2006; Yoshida, 1992; Preto and
Tremaine, 1999). Additional explicit symplectic schemes were devised for non-
separable Hamiltonians of specific forms Chin (2009); Tao (2016a, b); Benettin
et al. (2001). However, there is no explicit symplectic scheme for general
Hamiltonians. The Hamiltonian of an autonomous $N$-body system in uniformly
curved space, such as the surface of a sphere, is non-separable, and is not
amenable to known existing explicit schemes. The obstruction to separability
arises from the non-trivial metric coefficients in the kinetic term. We
identify a hierarchy in the structure of the metric of uniform and isotropic
spaces that allows the construction of an explicit, second order time
reversible symplectic scheme. We use iterative applications of this method to
obtain a fourth order integrator and evaluate its performance in comparison to
other methods.

##  II Constructing an explicit symplectic and time-reversible method

The goal of the present section is to construct an explicit, symplectic, and
time-reversible scheme that approximates the exact flow map of a Hamiltonian
with the following structure:

$H(\mathbf{z})=H_{1}(p_{\theta})+H_{2}(\theta,p_{\varphi})+H_{3}(\theta,\varphi)\;,$
(1)

with
$\mathbf{z}=(\mathbf{q},\mathbf{p})^{T}=(\theta,\varphi,p_{\theta},p_{\varphi})^{T}$.
The Hamiltonian $H(\mathbf{z})$ admits the splitting into the sum of several
terms, but it is nonetheless non-separable. The non-separability arises due to
the kinetic term $H_{2}(\theta,p_{\varphi})$, which prevents the complete
splitting between terms that only depend on the position coordinates, and
those that solely depend on the momentum coordinates. Such kinetic terms are
typical for Hamiltonians defined on isotropic and uniform curved geometries
(see SM). While the lack of full separability prevents the naïve
implementation of standard splitting techniques, e.g. Haruo (1990), the
hierarchical structure of the kinetic terms guides the design of a numerical
scheme with the desired features along similar lines. To see this, let us
start by examining Hamilton’s equations:

$\frac{d\mathbf{z}(t)}{dt}=\\{\mathbf{z},H(\mathbf{z})\\}=\hat{D}_{H}\mathbf{z},$
(2)

where the curly braces represent the Poisson bracket, given by

$\\{\mathbf{z},H(\mathbf{z})\\}=\sum^{2}_{l=1}\frac{\partial\mathbf{z}}{\partial
q_{l}}\frac{\partial H(\mathbf{z})}{\partial
p_{l}}-\frac{\partial\mathbf{z}}{\partial p_{l}}\frac{\partial
H(\mathbf{z})}{\partial q_{l}}\;.$ (3)

In eq.(2), the differential operator $\hat{D}_{H}$ is defined as
$\hat{D}_{H}f=\\{f,H(\mathbf{z})\\}$, with $f$ an arbitrary function.
Rewriting Hamilton’s equations in operator form allows the introduction of a
formal solution to the equations of motion Haruo (1990); Donnelly and Rogers
(2005), namely

$\mathbf{z}(\tau)=e^{\tau\hat{D}_{H}}\mathbf{z}(0)\;.$ (4)

From a computational stand point, the solution in eq. (4) is often
impractical. Its implementation requires knowing the form of the flow map
$e^{\tau\hat{D}_{H}}$, and therefore, requires knowing exactly how to solve
the equations of motion. It, nevertheless, motivates the use of Yoshida’s
method Haruo (1990) in the construction of an approximate scheme for the flow
map $e^{\tau\hat{D}_{H}}$. Recall that $H$ is given by the sum of three terms.
As a consequence, the operator $\hat{D}_{H}$ also allows for a similar
splitting, $\hat{D}_{H}=\hat{D}_{H_{1}}+\hat{D}_{H_{2}}+\hat{D}_{H_{3}}$
resulting in
$e^{\tau\hat{D}_{H}}=e^{\tau(\hat{D}_{H_{1}}+\hat{D}_{H_{2}}+\hat{D}_{H_{3}})}$.
We can thus aim at approximating the exact flow map $e^{\tau\hat{D}_{H}}$ as
the successive concatenation of the individual flow maps
$e^{\tau\lambda_{j}\hat{D}_{H_{i}}}$, with the number of concatenations $k$
and the value of the set of coefficients $\lambda_{j}=\\{c_{j},d_{j},f_{j}\\}$
determining the order of the integrator, $n$:

$\begin{split}&e^{\tau(\hat{D}_{H_{1}}+\hat{D}_{H_{2}}+\hat{D}_{H_{3}})}=\\\
&=\prod^{k}_{j=1}e^{\tau c_{j}\hat{D}_{H_{1}}}e^{\tau
d_{j}\hat{D}_{H_{2}}}e^{\tau
f_{j}\hat{D}_{H_{3}}}+\mathcal{O}(\tau^{n+1})\;.\end{split}$ (5)

Whenever the individual Hamiltonian terms $H_{i}$ are integrable, each flow
map $e^{\tau\hat{D}_{H_{j}}}$ is known exactly, and the method above can
automatically be interpreted as the composition of symplectic maps, giving
rise to a final method that is also symplectic Benettin et al. (2001);
Leimkuhler and Reich (2004a). In what follows, we show how to use Yoshida’s
framework to build explicit, symplectic and time-reversible methods of second
and forth order.

A second order method is constructed from eq.(5) by setting $k=3$,
$c_{1}=d_{1}=d_{2}=c_{3}=\frac{1}{2}$, $f_{1}=1$ and
$c_{2}=f_{2}=d_{3}=f_{3}=0$,

$\begin{split}&e^{\tau(\hat{D}_{H_{1}}+\hat{D}_{H_{2}}+\hat{D}_{H_{3}})}=\\\
&=e^{\frac{\tau}{2}\hat{D}_{H_{1}}}e^{\frac{\tau}{2}\hat{D}_{H_{2}}}e^{\tau\hat{D}_{H_{3}}}e^{\frac{\tau}{2}\hat{D}_{H_{2}}}e^{\frac{\tau}{2}\hat{D}_{H_{1}}}+\mathcal{O}(\tau^{3})\;.\end{split}$
(6)

To better understand the origins of eq.(6), we label
$\hat{X}=\frac{\tau}{2}\hat{D}_{H_{2}}$,
$\hat{Y}=\frac{\tau}{2}\hat{D}_{H_{1}}$ and $\hat{U}=\tau\hat{D}_{H_{3}}$,
which allows the right hand side of eq.(6) to be expressed as the product
$e^{\hat{Y}}e^{\hat{X}}e^{\hat{U}}e^{\hat{X}}e^{\hat{Y}}$. This product can be
further simplified to $e^{\hat{Z^{\prime}}}e^{\hat{T}}$, by setting
$e^{\hat{Z^{\prime}}}=e^{\hat{Y}}e^{\hat{X}}$ and
$e^{\hat{T}}=e^{\hat{U}}e^{\hat{Z}}=e^{\hat{U}}e^{\hat{X}}e^{\hat{Y}}$. We are
then left with solving the equation
$e^{\hat{Z^{\prime}}}e^{\hat{T}}=e^{\hat{W}}$, which immediately suggests the
use of the Baker-Campbell-Hausdorff (BCH) formula to obtain a solution for
$\hat{W}$. Given that the operators $\hat{X}$, $\hat{Y}$ and $\hat{U}$ are
proportional to $\tau$, we only need to take the terms in the BCH formula that
involve at most the product of two operators. Therefore, up to the desired
accuracy, we get the following solution for $\hat{W}$:

$\begin{split}\hat{W}&=\hat{Z^{\prime}}+\hat{T}+\frac{1}{2}[\hat{Z^{\prime}},\hat{T}]+\mathcal{O}(\tau^{3})\\\
&=2\hat{X}+2\hat{Y}+\hat{U}+\mathcal{O}(\tau^{3})\;.\end{split}$ (7)

All the commutators arising from the term $[\hat{Z^{\prime}},\hat{T}]$ cancel
each other exactly, and the next contributions stemming from the BCH formula
are already of the order $\tau^{3}$. We, thus, arrive at eq.(6). For this
result to be useful, nonetheless, we need to be able to write down the
individual flow maps $e^{\tau\hat{D}_{H_{i}}}$. We consider these to be given
by the elementary one-step Euler method Leimkuhler and Reich (2004a):

$e^{\tau\hat{D}_{H_{i}}}\mathbf{z}^{(j)}=\mathbf{z}^{(j+1)}=\mathbf{z}^{(j)}+\tau\hat{D}_{H_{i}}\mathbf{z}^{(j)}\;,$
(8)

with the action of the operator $\hat{D}_{H_{i}}$ as given in eq.(2). By
acting successively with the individual flow maps $e^{\tau\hat{D}_{H_{i}}}$ on
$\mathbf{z}^{(j)}$, as described in eq.(6), we obtain the following
integrator:

$\begin{split}&\theta^{(j+\frac{1}{2})}=\theta^{(j)}+\frac{\tau}{2}\frac{\partial
H_{1}}{\partial p_{\theta}}\big{(}p^{(j)}_{\theta}\big{)}\\\
&\varphi^{(j+\frac{1}{2})}=\varphi^{(j)}+\frac{\tau}{2}\frac{\partial
H_{2}}{\partial
p_{\varphi}}\big{(}\theta^{(j+\frac{1}{2})},p^{(j)}_{\varphi}\big{)}\\\
&p^{(j+1)}_{\varphi}=p^{(j)}_{\varphi}-\tau\frac{\partial
H_{3}}{\partial\varphi}\big{(}\theta^{(j+\frac{1}{2})},\varphi^{(j+\frac{1}{2})}\big{)}\\\
&p^{(j+1)}_{\theta}=p^{(j)}_{\theta}-\tau\Bigg{(}\frac{\partial
H_{3}}{\partial\theta}\big{(}\theta^{(j+\frac{1}{2})},\varphi^{(j+\frac{1}{2})}\big{)}+\\\
&+\frac{1}{2}\frac{\partial
H_{2}}{\partial\theta}\big{(}\theta^{(j+\frac{1}{2})},p^{(j)}_{\varphi}\big{)}+\frac{1}{2}\frac{\partial
H_{2}}{\partial\theta}\big{(}\theta^{(j+\frac{1}{2})},p^{(j+1)}_{\varphi}\big{)}\Bigg{)}\\\
&\theta^{(j+1)}=\theta^{(j+\frac{1}{2})}+\frac{\tau}{2}\frac{\partial
H_{1}}{\partial p_{\theta}}\big{(}p^{(j+1)}_{\theta}\big{)}\\\
&\varphi^{(j+1)}=\varphi^{(j+\frac{1}{2})}+\frac{\tau}{2}\frac{\partial
H_{2}}{\partial
p_{\varphi}}\big{(}\theta^{(j+\frac{1}{2})},p^{(j+1)}_{\varphi}\big{)}\;.\end{split}$
(9)

It then follows that the proposed numerical integrator, when executed in the
order of instructions listed in eq.(9), is entirely composed from explicit
steps. This greatly simplifies its numerical implementation.

We may write the integrator in eq.(9) in a more compact form, namely as the
set of instructions

$\begin{split}&\mathbf{z}^{(j+\frac{1}{2})}=\mathbf{z}^{(j)}+\frac{\tau}{2}J\nabla_{\mathbf{z}}H\big{(}\mathbf{q}^{(j+\frac{1}{2})},\mathbf{p}^{(j)}\big{)}\\\
&\mathbf{z}^{(j+1)}=\mathbf{z}^{(n+\frac{1}{2})}+\frac{\tau}{2}J\nabla_{\mathbf{z}}H\big{(}\mathbf{q}^{(j+\frac{1}{2})},\mathbf{p}^{(j+1)}\big{)}\;,\end{split}$
(10)

where $J$ is the canonical structure matrix, given by
$J=\begin{pmatrix}\mathbf{0}&\mathbf{1}\\\
-\mathbf{1}&\mathbf{0}\end{pmatrix}$, with $\mathbf{0}$ representing a
$2\times 2$ matrix filled with zeros, and $\mathbf{1}$ the $2\times 2$
identity matrix. Although the more compact notation used in eq.(10) hides the
fact that integrator can be made explicit, it allows us to deduce the
symplectic and time-reversible properties of the method. Written in the
compact notation of eq.(10), the method proposed in eq.(6) can be identified
with the generalised Störmer-Verlet/leapfrog method, which is known to be a
second order symplectic and time-reversible method Leimkuhler and Reich
(2004a); Hairer et al. (2006, 2003). Concatenating three such steps with
appropriately chosen time differences, namely

$\tilde{\tau}_{1}=\frac{\tau}{2-2^{1/3}},\quad\tilde{\tau}_{2}=-\frac{2^{1/3}\tau}{2-2^{1/3}},\quad\tilde{\tau}_{3}=\frac{\tau}{2-2^{1/3}},$

results in a fourth order integration scheme advancing the system by the time
step $\tau$ Haruo (1990).

##  III Comparison of numerical code

Figure 1: Comparison of the numerical integration schemes: Simulating the time
evolution of the Hamiltonian in equation 11, for special symmetrical initial
conditions. The exact solution, given in the SMV, displays a closed trajectory
whose projection onto the XY plane yields a 6-petaled shape. The schemes
compared are: A commercially available explicit Runge-Kutta (4,5)Dormand and
Prince (1980); Shampine and Reichelt (1997) (first column), an implicit
symplectic midpoint second order scheme (second column) Leimkuhler and Reich
(2004b), and the $2^{nd}$ and $4^{th}$ order symplectic explicit spherical
integrator (SESI) introduced in this work (third and fourth columns,
respectively). (a),(b),(c),(d) - Depict the projection of the trajectory of
one of the masses on the XY plane. The results of the numerical integrators
shown in (green and blue) dots overlaid by the exact solution, plotted in a
continuous (black) line. The shading (color) of the dots vary with time from
dark (green) to light (blue). (e),(f),(g),(h) - Display energy deviation,
where $\Delta E\equiv E-E_{0}$. (i),(j),(j),(k) - Present the numerical
angular momentum drift on $x,y,z$ axes, in blue, red, and yellow respectively.

To benchmark our integrator, we simulate the following Hamiltonian

$\begin{split}H=\sum_{i}\frac{1}{2m}\left(p_{i,\theta}^{2}+\left(\frac{p_{i,\phi}}{\sin\left(\theta_{i}\right)}\right)^{2}\right)\\\
+\sum_{i,j}\left(\frac{\sqrt{1-\left(\frac{2\sin\left(\frac{L_{ij}}{2R}\right)}{\sqrt{3}}\right)^{2}}-\cos\left(\theta_{0}\right)}{\frac{2\sin\left(\frac{L_{ij}}{2R}\right)}{\sqrt{3}}}\right)^{2},\end{split}$
(11)

which allows an explicit solution (See SMV for details). By fine tuning the
initial conditions we obtain a closed trajectory in real space, which enables
an easier visualization of the deviations from the exact solution in long time
trajectories. Figure 1 shows the comparison of the Symplectic Explicit
Spherical Integrator (SESI) derived here for both second order and fourth
order schemes, with the implicit midpoint symplectic second order scheme
Leimkuhler and Reich (2004b) (see SMV for details), as well as with a
commercially available explicit Runge-Kutta (RK) method based on the Dormand-
Prince (4,5) pair Dormand and Prince (1980); Shampine and Reichelt (1997). We
observe that both the $2^{nd}$ and the $4^{th}$ order SESI schemes demonstrate
a bounded error in energy, while the non symplectic RK(4,5) displays an
unbounded error in energy that grows with time. The error in angular momentum
is unbounded, yet better in the symplectic integrators than the RK(4,5).
Moreover, we observe that the $4^{th}$ order symplectic integrator conserves
the closed shape of the solution better than the other solvers, for the given
time window. The run time of the implicit symplectic solver was two orders of
magnitude longer than the comparable $2^{nd}$ order SESI.

##  IV discussion

Efficient and reliable long time numerical integration of Hamiltonian systems
favor explicit symplectic schemes. The kinetic term in many-body Hamiltonians
in uniformly curved space is non-seperable, precluding the usage of existing
explicit symplectic methods. We exploit the geometric hierarchy in the
structure of the kinetic term on the surface of the sphere, to produce an
explicit symplectic scheme. We note that while the resulting integrator is
symplectic and energy deviations remain bounded, other conserved quantities
arising from the symmetry of the metric, such as the angular momentum, are not
identically conserved.

These results may be generalized to higher dimensions as well as additional
Hamiltonians with a certain structure. The metrics of all uniformly curved
spaces of $N$-dimensions share a similar hierarchy to the one we exploited
above; namely that the kinetic term for each of the masses can be written as
the sum of $N$ terms $K_{i}$, where each $K_{i}$ depends on $p_{i}$ and on the
preceding coordinates $q_{1},q_{2},\dots,q_{i-1}$:

$K=\sum_{i=1}^{N}K_{i}(q_{1},q_{2},\dots,q_{i-1},p_{i}).$

Whenever the non-separable kinetic term displays this form, we can utilize a
similar maneuver to the one employed in this work to produce an explicit
symplectic scheme as expressed in equation (10): Advancing the $i^{th}$
coordinate, $q_{i}$, from time $t$ to time $t+\tau/2$ is performed using the
knowledge of $q_{k<i}$ at time $t+\tau/2$, and of $p_{i}$ at time $t$. Due to
the hierarchical structure described above, this could be carried out
explicitly. Advancing the corresponding momenta to time $t+\tau/2$ can then be
carried out explicitly as they too depend only on the momenta at time $t$ and
on the already calculated coordinates at time $t+\tau/2$. To obtain time-
reversibility we next perform the mirrored version of the above step: We note
that the hierarchy above implies that $\partial H/\partial q_{i}$ may only
depend on momenta terms $p_{k}$ such that $i<k$. Thus, starting with the
momenta of largest index, we may successively advance the momenta from time
$t+\tau/2$ to time $t+\tau$ using the knowledge of the coordinates at time
$t+\tau/2$ and the already calculated momenta at time $t+\tau$. We conclude by
advancing the coordinates from time $t+\tau/2$ to time $t+\tau$ by using the
coordinates at time $t+\tau/2$ and the already calculated momenta at time
$t+\tau$. We thus obtain an explicit, time-reversible, symplectic,
integration-scheme. The resulting second order scheme may be improved to yield
higher even orders using standard procedures Haruo (1990).

## References

  * Yoshida (1993) H. Yoshida, Qualitative and Quantitative Behaviour of Planetary Systems pp. 27–43 (1993).
  * Gladman et al. (1991) B. Gladman, M. Duncan, and J. Candy, Celestial Mechanics and Dynamical Astronomy 52, 221 (1991).
  * Channell and Scovel (1990) P. J. Channell and C. Scovel, Nonlinearity 3, 231 (1990).
  * Reich (1996) S. Reich, SIAM journal on numerical analysis 33, 475 (1996).
  * Haruo (1990) Y. Haruo, Physics letters A 150, 262 (1990).
  * Forest and Ruth (1990) E. Forest and R. D. Ruth, Physica D: Nonlinear Phenomena 43, 105 (1990).
  * Leimkuhler and Reich (2004a) B. Leimkuhler and S. Reich, _Simulating hamiltonian dynamics_ , 14 (Cambridge university press, 2004a).
  * Hairer et al. (2006) E. Hairer, C. Lubich, and G. Wanner, _Geometric numerical integration: structure-preserving algorithms for ordinary differential equations_ , vol. 31 (Springer Science & Business Media, 2006).
  * Yoshida (1992) H. Yoshida, in _Symposium-International Astronomical Union_ (Cambridge University Press, 1992), vol. 152, pp. 407–411.
  * Preto and Tremaine (1999) M. Preto and S. Tremaine, The Astronomical Journal 118, 2532 (1999).
  * Chin (2009) S. A. Chin, Physical Review E 80, 037701 (2009).
  * Tao (2016a) M. Tao, Journal of Computational Physics 327, 245 (2016a).
  * Tao (2016b) M. Tao, Physical Review E 94, 043303 (2016b).
  * Benettin et al. (2001) G. Benettin, A. M. Cherubini, and F. Fassò, SIAM Journal on Scientific Computing 23, 1189 (2001).
  * Donnelly and Rogers (2005) D. Donnelly and E. Rogers, American Journal of Physics 73, 938 (2005).
  * Hairer et al. (2003) E. Hairer, C. Lubich, and G. Wanner, Acta numerica 12, 399 (2003).
  * Dormand and Prince (1980) J. R. Dormand and P. J. Prince, Journal of computational and applied mathematics 6, 19 (1980).
  * Shampine and Reichelt (1997) L. F. Shampine and M. W. Reichelt, SIAM journal on scientific computing 18, 1 (1997).
  * Leimkuhler and Reich (2004b) B. Leimkuhler and S. Reich, _Simulating hamiltonian dynamics_ , 14 (Cambridge university press, 2004b).

##  V  Supplementary Material:
Explicit, time-reversible and symplectic integrator for Hamiltonians in
isotropic uniformly curved geometries

##  VI Order of the integrator

For the purpose of verifying the order of the integration scheme we run the
simulation to a fixed final time $t=10$ using different time steps $\Delta t$,
and examine the difference between the value of the final position $x(t=10)$
for the different time steps. Figure 2 shows the difference in the final
position $x(t=10)$ between two runs with nearby time steps, as a function of
the time step magnitude indicating that the SESI 2 and SESI 4 schemes are
indeed second and fourth order respectively.

Figure 2: Numerical schemes convergence rate.
$\vec{x}_{\Delta_{i}}\equiv\left(\theta(t=10),\phi(t=10),p_{\theta}(t=10),p_{\phi}(t=10)\right)$
calculated using a time step of $\Delta_{i}$. Filled dots correspond to SESI
2; hollow diamonds correspond to SESI 4.

##  VII Finding an explicit solution

To benchmark our numerical schemes’ accuracy we seek to find a potential that
will allow an explicit analytic solution. We seek this solution as a highly
symmetric quasi-one-dimensional solution for a specially designed interaction.

We start by considering the three mass system on the sphere. Each of the
masses position is given by

$\mathbf{r}_{i}=(\sin(\theta_{i})\cos(\phi_{i}),\sin(\theta_{i})\sin(\phi_{i}),\cos(\theta_{i})).$

We assume the sought solution is associated with the three fold symmetry for
the problem, i.e.

$\theta_{i}=\theta,\quad\phi_{1}=\phi,\quad\phi_{2}=\phi+2\pi/3,\quad\phi_{3}=\phi+4\pi/3.$

The geodesic (or any other) distance between every two particles becomes in
this highly symmetric case only a function of $\theta$. The real space
distance simply measures $\sqrt{3}\sin(\theta)$. The kinetic terms are
proportional to $\dot{\theta}^{2}+\sin(\theta)^{2}\dot{\phi}^{2}$. The
Lagrangian of the system thus reduces (after dividing by the constant
$3m/2R^{2}$):

$\mathcal{L}=\dot{\theta}^{2}+\sin(\theta)^{2}\dot{\phi}^{2}-V(\theta).$

The conserved angular momentum of the system is given by

$L=\sin(\theta)^{2}\dot{\phi}.$

The Euler Lagrange equation thus yields

$2\sin(\theta)^{-2}\cot(\theta)L^{2}-V^{\prime}(\theta)-2\ddot{\theta}=0,$

which could be multiplied by $\dot{\theta}$ and explicitly integrated to give

$V(\theta)+L^{2}\cot(\theta)^{2}+\dot{\theta}^{2}=c_{1}.$

We may now identify the system’s conserved energy and rewrite the equation as

$V(\theta)+L^{2}(\sin(\theta)^{-2}-1)+\dot{\theta}^{2}=c_{1}=E_{0}-L^{2}.$

Setting $\psi=\cos(\theta)$ we have

$\dot{\psi}^{2}=E_{0}-L^{2}-E_{0}\psi^{2}-U(\psi)$

where $U(\psi)=U(\cos(\theta))=V(\theta)\sin(\theta)^{2}$. We are thus
required to integrate

$\int\frac{d\psi}{\sqrt{E_{0}-L^{2}-E_{0}\psi^{2}-U(\psi)}}=t.$

What remains is to choose $V(\theta)$ and the corresponding $U(\psi)$ to
render the equations solvable. One very convenient choice for the interactions
reads

$V(\theta)=(\cos(\theta)-\cos(\theta_{0}))^{2}/\sin(\theta)^{2}.$

While this potential diverges near the origin, it gives a nice parabolic
behavior in the vicinity of $\theta_{0}$, which is where we make use of it.
Moreover, this potential yields $U(\psi)=(\psi-\psi_{0})^{2}$. It is now
straightforward to explicitly calculate the full time evolution for $\psi$:

$\psi=\frac{\psi_{0}}{1+E_{0}}+\sqrt{\frac{E_{0}-L^{2}-\frac{E_{0}\psi_{0}^{2}}{1+E_{0}}}{1+E_{0}}}\sin(\sqrt{1+E_{0}}t).$

Writing the three dimensional distance between every two masses as
$L^{3D}=\sqrt{3}R\sin(\theta)$ and its relation to the corresponding geodesic
distance on the sphere through

$L^{S}=2R\arcsin(L^{3D}/2R),$

yields the potential in the form (11).

##  VIII Simulation details

We set $\theta_{0}=\frac{\pi}{4}$. In order to get a closed shape trajectory
in real space, we choose $p_{\theta}=\frac{1}{4}$ and
$p_{\phi}=0.1727174029854043$ such that
$\int_{0}^{\frac{2\pi}{\sqrt{E_{0}+1}}}\dot{\varphi}dt=\frac{2\pi}{6}$. The
integration time step is $\tau=0.1$ and the total integration time is
$t=800=8000\tau$.

The implicit midpoint symplectic second order integration scheme is [19]

$\begin{split}p_{n+\frac{1}{2}}=p_{n}-\frac{\tau}{2}\frac{\partial H}{\partial
q}\left(q_{n+\frac{1}{2}},p_{n+\frac{1}{2}}\right)\\\
q_{n+\frac{1}{2}}=q_{n}+\frac{\tau}{2}\frac{\partial H}{\partial
p}\left(q_{n+\frac{1}{2}},p_{n+\frac{1}{2}}\right)\\\
q_{n+1}=q_{n+\frac{1}{2}}+\frac{\tau}{2}\frac{\partial H}{\partial
p}\left(q_{n+\frac{1}{2}},p_{n+\frac{1}{2}}\right)\\\
p_{n+1}=p_{n+\frac{1}{2}}-\frac{\tau}{2}\frac{\partial H}{\partial
q}\left(q_{n+\frac{1}{2}},p_{n+\frac{1}{2}}\right)\end{split}$ (12)