# Supplementary information: Data-driven classification of single cells by
their non-Markovian motion

###### Contents

  1. I Video captions
  2. II Effective kernel follows uniquely from correlation functions
  3. III Green’s function is given in terms of positional two-point correlation function
  4. IV Long-time MSDs
  5. V Relation between the oscillation periods of the friction kernel and the VACF
  6. VI Derivation of the analytical expression for the MSD
  7. VII Localization noise fit
  8. VIII Cluster analysis in lower dimensions
  9. IX Comparing friction kernel expressions eq. (7) and eq. (26)
  10. X Exemplary non-equilibrium model describing CR motion
  11. XI Information on the cell orientation is contained in the cell-center trajectories
  12. XII Effects of smoothing trajectory data
  13. XIII Markovian embedding of the friction kernel eq. (26)

##  I Video captions

Video 1.: A representative synchro. High-speed video microscopy at 500 frames
per second obtained by phase-contrast imaging of a synchro CR cell showing the
planar and synchronous breaststroke motion of the flagella. Between t $\sim
330-390\,ms$ the synchronous beat of the flagella exhibits a phase slip,
meaning the synchronicity of the flagella is disturbed in that short time
interval.
Video 2.: A representative wobbler. High-speed video microscopy at 500 frames
per second obtained by phase-contrast imaging of a CR cell which paddles the
flagella in an asynchronous and irregular manner resulting in the wobbling
motion of the cell body.

##  II Effective kernel follows uniquely from correlation functions

We can write the GLE with an arbitrary potential $U(x(t))$ as

$\ddot{x}(t)=-\nabla
U(x(t))-\int_{t_{0}}^{t}\Gamma_{v}(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}+F_{R}(t)\,.$
(S1)

Multiplying eq. (S1) by $v(t_{0})$, inserting
$\Gamma(|t|)=\Gamma_{R}(t)=\Gamma_{v}(|t|)$ and averaging the entire equation
leads to the Volterra equation

$\dfrac{d}{dt}C_{vv}(t)=-C_{v\nabla
U}(t)-\int_{0}^{t}\Gamma(t-s)C_{vv}(s)ds\,,$ (S2)

using that the random force is not correlated with the initial velocity at the
projection time, i.e. $\langle v(t_{0})F_{R}(t)\rangle=0$ and setting
$t_{0}=0$. Volterra equations are well-behaved and for given functions
$\Gamma(t)$ and $C_{v\nabla U}(t)$ one can find a solution of eq. (S2) in
terms of $C_{vv}(t)$. Inversely, for given correlation functions $C_{vv}(t)$
and $C_{v\nabla U}(t)$ one can solve eq. (S2) in terms of the friction kernel
$\Gamma(t)$ by Laplace transformation. With the definition of the Laplace
transform of a function $f(t)$

$\hat{f}(q)=\int_{0}^{\infty}f(t)e^{-qt}dt,$ (S3)

the unique solution for $\Gamma(t)$ in Laplace space is given by is given by

$\hat{\Gamma}(q)=-\frac{1}{\hat{C}_{vv}(q)}\left(q\hat{C}_{vv}(q)-C_{vv}(0)+\hat{C}_{v\nabla
U}(q)\right)\,.$ (S4)

The existence of a unique friction kernel for any given input correlation
functions $C_{vv}(t)$ and $C_{v\nabla U}(t)$ assures that one can always find
an effective friction kernel $\Gamma(t)$ that, when employed in the GLE and
using $\Gamma(|t|)=\Gamma_{R}(t)=\Gamma_{v}(|t|)$, reproduces the two-point
correlation functions. Thus, every non-equilibrium model described by eq. (1)
with $\Gamma_{R}(t)\neq\Gamma_{v}(|t|)$ can be mapped on an effective model
with $\Gamma(|t|)=\Gamma_{R}(t)=\Gamma_{v}(|t|)$, because the Green’s
function, which completely describes a Gaussian process, is solely given in
terms of the positional two-point correlation function, see SI Sec. III.

##  III Green’s function is given in terms of positional two-point correlation
function

In the following we map a non-equilibrium process described by the GLE eq.
(S1) with the random force fluctuations given by eq. (2) and a harmonic
potential $U(x)=Kx^{2}/2$ on an effective model with an effective friction
kernel $\Gamma(t)$ fulfilling $\langle
F_{R}(0)F_{R}(t)\rangle=k_{B}T\Gamma(|t|)$ and an effective harmonic potential
${U_{\rm{eff}}(x)=K_{\rm{eff}}x^{2}/2}$. For this we show that the Green’s
function for the non-equilibrium system eqs. (S1), (2) is given in terms of
the positional autocorrelation function.

We define the average over the random force by the path integral

$\langle
X\rangle=\mathcal{N}\int\mathcal{D}\tilde{F}_{R}(\cdot)X\exp\left(-\frac{1}{4\pi
B}\int_{-\infty}^{\infty}d\omega\int_{-\infty}^{\infty}d\omega^{\prime}\frac{\tilde{F}_{R}(\omega)\delta(\omega+\omega^{\prime})\tilde{F}_{R}(\omega^{\prime})}{\tilde{\Gamma}_{R}(\omega)}\right)\,.$
(S5)

We also define the generating functional as

$\displaystyle\mathcal{Z}[\tilde{h}]=$
$\displaystyle\mathcal{N}\int\mathcal{D}\tilde{F}_{R}(\cdot)\exp\left(-\frac{1}{4\pi
B}\int_{-\infty}^{\infty}d\omega\frac{\tilde{F}_{R}(\omega)\tilde{F}_{R}(-\omega)}{\tilde{\Gamma}_{R}(\omega)}+\int_{-\infty}^{\infty}d\omega\tilde{h}(-\omega)\tilde{F}_{R}(\omega)\right)$
$\displaystyle=$ $\displaystyle\exp\left(2\pi
B\int_{-\infty}^{\infty}d\omega\tilde{h}(\omega)\tilde{h}(-\omega)\tilde{\Gamma}_{R}(\omega)\right)$
(S6)

where $\tilde{h}(\omega)$ is a generating field. Equation (S6) has the
functional derivatives

$\displaystyle\dfrac{\delta\mathcal{Z}[\tilde{h}]}{\delta\tilde{h}(\omega)}\mid_{\tilde{h}=0}=\langle\tilde{F}_{R}(\omega)\rangle=0$
(S7)
$\displaystyle\dfrac{\delta^{2}\mathcal{Z}[\tilde{h}]}{\delta\tilde{h}(\omega)\delta\tilde{h}(\omega^{\prime})}\mid_{\tilde{h}=0}=\langle\tilde{F}_{R}(\omega)\tilde{F}_{R}(\omega^{\prime})\rangle=2\pi
B\delta(\omega+\omega^{\prime})\tilde{\Gamma}_{R}(\omega)\,,$ (S8)

where $\tilde{\Gamma}_{R}(\omega)=\tilde{\Gamma}_{R}(-\omega)$. For harmonic
potential ${U(x)=Kx^{2}/2}$ we can Fourier-transform the GLE eq. (S1) by which
we obtain the form of eq. (S28) with the response function given by

$\tilde{\chi}(\omega)=\left(K-\omega^{2}+i\omega\tilde{\Gamma}_{v}^{+}(\omega)\right)^{-1}\,.$
(S9)

The Green’s function can be written as

$G(x,t\mid x^{\prime},t^{\prime})=\frac{P(x,t\mathchar
24635\relax\;x^{\prime},t^{\prime})}{P(x^{\prime},t^{\prime})},$ (S10)

which is the ratio of the joint two-point probability distribution to find $x$
at $t$ and $x^{\prime}$ at $t^{\prime}$ and the base probability to find
$x^{\prime}$ at $t^{\prime}$. Those probabilities can be expressed as

$P(x^{\prime},t^{\prime})=\langle\delta(x^{\prime}-x(t^{\prime}))\rangle\,,\,P(x,t\mathchar
24635\relax\;x^{\prime},t^{\prime})=\langle\delta(x^{\prime}-x(t^{\prime}))\delta(x-x(t))\rangle\,.$
(S11)

From eq. (S28) we can write the solution of the GLE eq. (S1) using the linear
response function $\chi(t)$ as

$x(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{i\omega
t}\tilde{\chi}(\omega)\tilde{F}_{R}(\omega).$ (S12)

Inserting this result into the probability distributions eq. (S11) we obtain

$P(x^{\prime},t^{\prime})=\int_{-\infty}^{\infty}\frac{dp}{2\pi}e^{-ipx^{\prime}}\Biggl{\langle}\exp\left(ip\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\tilde{\chi}(\omega)\tilde{F}_{R}(\omega)e^{i\omega
t^{\prime}}\right)\Biggl{\rangle}$ (S13)

and

$\displaystyle P(x,t\mathchar
24635\relax\;x^{\prime},t^{\prime})=\int_{-\infty}^{\infty}\frac{dp}{2\pi}\int_{-\infty}^{\infty}\frac{dq}{2\pi}e^{-iqx-
ipx^{\prime}}$
$\displaystyle\Biggl{\langle}\exp\Biggl{(}iq\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\tilde{\chi}(\omega)\tilde{F}_{R}(\omega)e^{i\omega
t}$
$\displaystyle+ip\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\tilde{\chi}(\omega)\tilde{F}_{R}(\omega)e^{i\omega
t^{\prime}}\Biggl{)}\Biggl{\rangle}.$ (S14)

By comparing eqs. (S13), (S14) to eq. (S6) we find the generating field
$\tilde{h}(\omega)$ for the base probability $P(x^{\prime},t^{\prime})$

$\tilde{h}_{b}(-\omega)=\frac{ip}{2\pi}\tilde{\chi}(\omega)e^{i\omega
t^{\prime}}$ (S15)

and for the joint probability $P(x,t\mathchar
24635\relax\;x^{\prime},t^{\prime})$ we find

$\tilde{h}_{j}(-\omega)=\frac{i}{2\pi}\tilde{\chi}(\omega)\left(qe^{i\omega
t}+pe^{i\omega t^{\prime}}\right)\,,$ (S16)

where the subscript $b$ stands for the base probability and $j$ for the joint
probability. Using that the position correlation function can be written as

$C_{xx}(t-t^{\prime})=B\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\tilde{\chi}(\omega)\tilde{\Gamma}_{R}(\omega)\tilde{\chi}(-\omega)e^{i\omega(t-t^{\prime})},$
(S17)

we arrive at the Gaussian form

$\displaystyle P(x^{\prime},t^{\prime})=$
$\displaystyle\int_{-\infty}^{\infty}\frac{dp}{2\pi}\exp\left(-C_{xx}(0)p^{2}-ipx^{\prime}\right)$
$\displaystyle=$ $\displaystyle\sqrt{\frac{1}{2\pi
C_{xx}(0)}}\exp\left(-\frac{x^{\prime 2}}{2C_{xx}(0)}\right)$ (S18)

and

$\displaystyle P(x,t\mathchar
24635\relax\;x^{\prime},0)=\int_{-\infty}^{\infty}\frac{dp}{2\pi}\int_{-\infty}^{\infty}\frac{dq}{2\pi}\exp\left(-\frac{q^{2}+p^{2}}{2}C_{xx}(0)-qpC_{xx}(t)-iqx-
ipx^{\prime}\right)$ $\displaystyle=$ $\displaystyle
2\sqrt{\frac{1}{2(C_{xx}(0)^{2}-C_{xx}(t)^{2})}}\exp\left(-\frac{(x^{\prime}+x)^{2}}{4(C_{xx}(0)+C_{xx}(t))}-\frac{(x^{\prime}-x)^{2}}{4(C_{xx}(0)-C_{xx}(t))}\right)\,.$
(S19)

For simplicity we set $t^{\prime}=0$ as the position correlation function only
depends on the time difference $t-t^{\prime}$. Inserting eqs. (S18) (S19) into
eq. (S10), we finally obtain the position Green’s function as a Gaussian

$G(x,t\mid
x^{\prime},0)=\frac{\exp(-\frac{(x-\mu)^{2}}{2\sigma^{2}})}{\sqrt{2\pi\sigma^{2}}}$
(S20)

with mean and standard deviation

$\mu=\mu(x^{\prime},t)=\frac{C_{xx}(t)}{C_{xx}(0)}x^{\prime}\,,\,\sigma^{2}=\sigma^{2}(t)=C_{xx}(0)\left(1-\frac{C_{xx}(t)^{2}}{C_{xx}(0)^{2}}\right)\,.$
(S21)

Thus the Green’s function is entirely described in terms of the two-point
positional correlation function $C_{xx}(t)$, note that
$C_{vv}(t)=-\dfrac{d^{2}}{dt^{2}}C_{xx}(t)$.

We next calculate the effective harmonic potential strength $K_{\rm{eff}}$.
The condition to find an effective description of the GLE eq. (S1) arises from
setting eq. (S17) equal to its effective version

$\tilde{\chi}(\omega)\tilde{\Gamma}_{R}(\omega)\tilde{\chi}(-\omega)=\tilde{\chi}_{\rm{eff}}(\omega)\tilde{\Gamma}(\omega)\tilde{\chi}_{\rm{eff}}(-\omega)\,,$
(S22)

which assures that the effective model has the same Green’s function as the
non-equilibrium model, as shown by eqs. (S20) and (S21). Here,
$\tilde{\chi}_{\rm{eff}}(\omega)$ is the effective response function and
$\Gamma(t)$ is the effective kernel. Inserting eq. (S9) into eq. (S22) leads
to

$\frac{\tilde{\Gamma}_{R}(\omega)}{(K-\omega^{2}+i\omega\tilde{\Gamma}_{v}^{+}(\omega))(K-\omega^{2}-i\omega\tilde{\Gamma}_{v}^{+}(-\omega))}=\frac{\tilde{\Gamma}(\omega)}{(K_{\rm{eff}}-\omega^{2}+i\omega\tilde{\Gamma}_{+}(\omega))(K_{\rm{eff}}-\omega^{2}-i\omega\tilde{\Gamma}_{+}(-\omega))}\,.$
(S23)

Defining the potential of mean force for the effective system according to

$U_{\rm{eff}}(x)=-B\ln(p(x))=K_{\rm{eff}}x^{2}/2,$ (S24)

with $p(x)$ being the probability to observe $x$, and comparing eq. (S24) to
the long time limit of the non-equilibrium distribution eq. (S20), namely

$p(x)=\lim_{t\to\infty}G(x,t\mid
x^{\prime},0)=\frac{\exp(-\frac{x^{2}}{2C_{xx}(0)})}{\sqrt{2\pi
C_{xx}(0)}}\,,$ (S25)

we find the effective harmonic potential strength to be

$K_{\rm{eff}}=\frac{B}{C_{xx}(0)}\,.$ (S26)

Note that the effective potential will in general deviate from the potential
in the non-equilibrium GLE, i.e. $K_{\rm{eff}}\neq K$. The present calculation
can also be done for a general multi-point distribution roland_neq_2023b .

##  IV Long-time MSDs

The trajectories analyzed in the main text (examples shown in Figs. 2a,b) are
recorded with a 40 fold magnification at 500 frames per second, which leads to
trajectories of lengths up to $14\,s$ (Fig. 2h). Here, we analyze a different
data set to resolve the long time behavior of the MSD of the CR cells, which
we obtain by using a lower magnification (tenfold) and longer time step of
$0.02\,s$. A lower magnification results in a larger observation window, which
in turn allows recording longer trajectories. The decreased localization
precision and the longer time step do not allow to observe the details of the
oscillatory features of the motion. Hence, we do not apply the cluster
analysis relying on the extraction of the oscillating friction kernels, but
rather use an approximate distinction by the mean squared velocities $B$. We
assign cells with $B>110\,\mu m^{2}/s^{2}$ to be putative wobblers and cells
with $B<35\,\mu m^{2}/s^{2}$ to be putative synchros. In this way, we avoid
the assignment of cells with intermediate mean squared velocities $B$, for
which this approximate classification scheme is unreliable. For the slow cells
we additionally exclude cells for which the MSD stays below $1\,\mu m^{2}$
during the entire trajectory. These cells are most likely stuck or have some
defects prohibiting them to move. This distinction into putative wobblers and
synchros is biased by our choice of the upper and lower bounds for $B$, which
we choose according to the high resolution fastest synchro and slowest wobbler
when using only every tenth step of the high resolution data in order to
obtain comparable velocities with a time step of $0.02\,s$. Nevertheless, in
Fig. S1 we demonstrate that the present approximate classification leads to
extended MSDs (black lines) that match closely the MSDs shown in the man text
based on the high-resolution videos (gray lines). The extended MSDs for
synchros exhibit a transition from the ballistic to the long time diffusive
regime at a time of about $t\approx 2$s, whereas for the wobblers the
ballistic regime extends over more than ten seconds. The prefactors of the MSD
in the diffusive regimes, i.e. the diffusivity $D$, is for the synchros shown
in Fig. S1.

Figure S1: MSDs from low-resolution data with 10x magnification and a time
step of $\Delta=0.02\,s$ of (a) wobblers and (b) synchros according to an
approximate classification scheme explained in the text. The average of the
colored single-cell MSDs is shown in black, the grey line represents the
average MSD from Fig. 3 with higher temporal and spatial resolution of 40x and
$\Delta=0.002\,s$.

##  V Relation between the oscillation periods of the friction kernel and the
VACF

In Fig. 4 in the main text one can note a subtle difference between the
oscillation period of the friction kernel and the oscillation period of the
VACF when comparing Fig. 4b to 4c or Fig. 4e to 4f. Since we derive an
analytical expression for the VACF in eq. (S52) for the friction kernel given
in eq. (7), we can evaluate the dependence of the VACF oscillation frequency,
denoted by $\omega_{vv}$, on the kernel frequency $\Omega$. From eq. (S52) we
know that there are three complex frequencies in the VACF that can lead to
oscillations, where the frequency $\omega_{i}$ with the smallest real part
dictates $\omega_{vv}$. In Fig. S2a we show the dependence of $\omega_{vv}$ on
$\Omega$ for different $b$ and find the VACF frequency $\omega_{vv}$ to be
constant for small kernel frequencies $\Omega$ followed by a transition to the
two frequencies being the same at a threshold value
$\Omega_{t}=\omega_{vv}(\Omega\rightarrow 0)$, which is defined by the
constant value $\omega_{vv}(\Omega\rightarrow 0)$ for small kernel frequencies
intersecting with $\omega_{vv}=\Omega$. This threshold value depends on the
kernel parameters $a$, $b$, $\tau$. It can be described by

$\Omega_{t}=\begin{cases}\ \sqrt{b}&,\text{for}\
b>\max(a^{2}/4,\tau^{-2}/4)\\\ \ a/2&,\text{for}\ b<a^{2}/4\land\tau>a^{-1}\\\
\tau^{-1}/2&,\text{for}\ b<a^{2}/4\land\tau<a^{-1}\,,\end{cases}$ (S27)

which follows from Figs. S2a-c.

Figure S2: (a) The dependence of the VACF frequency $\omega_{vv}$ on the
kernel frequency $\Omega$ of the kernel eq. (7) is shown for different kernel
oscillation amplitudes $b$ and fixed $\delta$-peak amplitude $a=100\tau^{-1}$.
The black line indicates the linear behavior $\omega_{vv}=\Omega$. (b) The
dependence of the VACF frequency $\omega_{vv}$ on the kernel oscillation
amplitude $b$ is shown for varying kernel $\delta$-peak amplitude $a$ and
$\Omega=\tau^{-1}$, which represents the low kernel frequency plateau value
shown in (a) for small $\Omega$. The horizontal lines represent
$\omega_{vv}=a/2$ for the respective values of $a$, the diagonal line
represents $\omega_{vv}=\sqrt{b}$ and the vertical lines represent the
intersection of the other lines, which mark a minimum in the parameter $b$
positioned at $b=a^{2}/4$. (c) The dependence of the VACF frequency
$\omega_{vv}$ on the kernel decay time $\tau$ shown for different
$\delta$-peak amplitudes $a$ and fixed $b=10^{5}\Omega^{2}$. The regime shown
for $\tau/\Omega^{-1}$ represents the low kernel frequency plateau value shown
in (a) for small $\Omega\tau<100$. The black line indicates the scaling
$\omega_{vv}=\tau^{-1}/2$.

Since the median kernel frequency of the CR cells is of the order ${\Omega\sim
100\,\tau^{-1}}$ and the oscillation amplitude is of the order $b\sim
10^{4}\tau^{-2}$, the data lies in the transition regime between the constant
VACF frequency to the linear regime $\omega_{vv}=\Omega$, as can be seen in
Fig. S2a. From Fig. S2a we know that the VACF frequency is always greater or
equal to the kernel frequency, i.e. $\omega_{vv}\geq\Omega$. Moreover, the
oscillation amplitudes of CR usually fulfill $b>a^{2}/4$ and the inverse decay
time $\tau^{-1}$ is small compared to the other kernel parameters, such that
$a>\tau^{-1}$ and $b>\tau^{-2}/4$ are always fulfilled. Therefore, the cells
exhibit parameters for which $b>\max(a^{2}/4,\tau^{-2}/4)$ holds. Eq. (S27)
thus tells us that the threshold frequency is given by $\Omega_{t}=\sqrt{b}$.
Combining this knowledge with $\omega_{vv}\geq\Omega$ leads to
$\omega_{vv}\geq\sqrt{b}$. Since the CR data lies in the transition regime
between $\omega_{vv}=\sqrt{b}$ and $\omega_{vv}=\Omega$, the VACF frequency
$\omega_{vv}$ is strongly influenced by the oscillation amplitude $b$. In
order to further estimate the exhibited parameter range, the reader is
referred to Fig. 5 of the main text, which includes all parameters of all
cells.

The fact that the VACF frequency $\omega_{vv}$ increases with increasing
kernel oscillation amplitude $b$ can be interpreted in the context of the
Markovian embedding model in Sec. XIII, where the coupling strength of the
harmonically coupled parts is proportional to the oscillation amplitude $b$,
as shown in eqs. (S75), (S77), which leads to faster oscillations for higher
$b$. Since every oscillation leads to a net forward motion of the cell, it
leads to fast net forward motion for cells if they exhibit a large kernel
oscillation amplitude $b$ compared to the squared inverse decay time
$\tau^{-2}$ and the squared $\delta$-amplitude $a^{2}$, which is the case in
our data set (see Fig. 5). By realizing that the decay time $\tau$ is related
to how long a cell ’remembers’ its past trajectory, one can conclude that a
fast decay of the oscillation would lead to a cell not being able to maintain
its direction for a long time. Thus, it is useful to achieve fast net forward
motion if the decay time is longer than the oscillation period, i.e. if
$\tau\Omega>1$. By net forward motion we mean directed motion on relevant time
scales, since of course in the long-time limit the cells exhibit diffusion as
shown in Fig. S1. Both criteria, the oscillation amplitude being large, i.e.
$b>\max(\tau^{-2},a^{2})$ and the decay time being long compared to the
oscillation period, i.e. $\tau\Omega>1$, are fulfilled by the CR cell
parameters shown in Fig. 5, which means that CR cells exhibit a range of
kernel parameters that is suited to achieve fast net forward motion.

##  VI Derivation of the analytical expression for the MSD

In Fourier space the GLE eq. (1) can be written as

$\tilde{x}(\omega)=\tilde{\chi}(\omega)\tilde{F}_{R}(\omega)\,,$ (S28)

where $\tilde{\chi}(\omega)$ is the Fourier transform of the position response
function, which takes the form

$\tilde{\chi}(\omega)=\left(-\omega^{2}+i\omega\tilde{\Gamma}_{v}^{+}(\omega)\right)^{-1}$
(S29)

with $\tilde{\Gamma}_{v}^{+}(\omega)$ being the half sided Fourier transform
of the friction kernel defined via

$\tilde{\Gamma}_{v}^{+}(\omega)=\int_{0}^{\infty}e^{-i\omega
t}\Gamma_{v}(t)dt\,.$ (S30)

The Fourier transform of the position correlation function $C_{xx}(t)=\langle
x(0)x(t)\rangle$ can be written as

$\tilde{C}_{xx}(\omega)=B\tilde{\chi}(\omega)\tilde{\Gamma}_{R}(\omega)\tilde{\chi}(-\omega)\,,$
(S31)

where we made use of

$\begin{split}&\tilde{C}_{xx}(\omega)=\int_{-\infty}^{\infty}dte^{-i\omega
t}\langle x(0)x(t)\rangle\\\ &=\int_{-\infty}^{\infty}e^{-i\omega
t}dt\int_{\infty}^{\infty}e^{i\omega
t}\frac{d\omega}{2\pi}\int_{-\infty}^{\infty}\frac{d\omega^{\prime}}{2\pi}\langle\tilde{x}(\omega)\tilde{x}(\omega^{\prime})\rangle\end{split}$
(S32)

and

$\langle\tilde{F}_{R}(\omega)\tilde{F}_{R}(\omega^{\prime})\rangle=2\pi
B\delta(\omega+\omega^{\prime})\tilde{\Gamma}_{R}(\omega^{\prime})\,.$ (S33)

Equation S31 can be rewritten as

$\tilde{C}_{xx}(\omega)=\frac{B}{i\omega}\left(\tilde{\chi}(\omega)-\tilde{\chi}(-\omega)\right)\frac{\tilde{\Gamma}_{R}(\omega)}{\tilde{\Gamma}_{v}(\omega)}\,,$
(S34)

which leads to the MSD

$\begin{split}C_{\textrm{MSD}}(t)&=2(C_{xx}(0)-C_{xx}(t))\\\
&=B\int_{-\infty}^{\infty}\frac{d\omega}{\pi}\frac{e^{i\omega
t}-1}{i\omega}\left(\tilde{\chi}(\omega)-\tilde{\chi}(-\omega)\right)\frac{\tilde{\Gamma}_{R}(\omega)}{\tilde{\Gamma}_{v}(\omega)}\end{split}\,.$
(S35)

For the decaying oscillation model friction kernel of eq. (7) and using the
effective form $\Gamma(|t|)=\Gamma_{R}(t)=\Gamma_{v}(|t|)$ one obtains the
half sided Fourier transform of the kernel as

$\tilde{\Gamma}_{+}(\omega)=\int_{0}^{\infty}\left(2a\delta(t)+be^{-t/\tau}\cos(\Omega
t)\right)e^{-i\omega
t}dt=a+b\left(\frac{i\omega+\frac{1}{\tau}}{\Omega^{2}+(i\omega+\frac{1}{\tau})^{2}}\right)\,.$
(S36)

Inserting the result of eq. (S36) into eq. (S29) we find the response function

$\tilde{\chi}(\omega)=\left(-\omega^{2}+i\omega a+i\omega
b\left(\frac{i\omega+\frac{1}{\tau}}{\Omega^{2}+(i\omega+\frac{1}{\tau})^{2}}\right)\right)^{-1}\,,$
(S37)

which by inserting into eq. (S35) leads to

$C_{\textrm{MSD}}(t)=B\int_{-\infty}^{\infty}\frac{d\omega}{\pi}\frac{(e^{i\omega
t}-1)(k_{1}+k_{2}\omega^{2}+k_{3}\omega^{4})}{\omega^{2}(c_{1}+c_{2}\omega^{2}+c_{3}\omega^{4}+\tau^{4}\omega^{6})}\,,$
(S38)

where the constants $c_{i}$ and $k_{i}$ are given by

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle(a+b\tau+a\Omega^{2}\tau^{2})^{2}$ (S39) $\displaystyle c_{2}$
$\displaystyle=$
$\displaystyle(1+2(a^{2}-b+\Omega^{2})\tau^{2}+2ab\tau^{3}+(-2a^{2}\Omega^{2}+(b+\Omega^{2})^{2})\tau^{4})$
(S40) $\displaystyle c_{3}$ $\displaystyle=$
$\displaystyle\tau^{2}(2+(a^{2}-2(b+\Omega^{2}))\tau^{2})$ (S41)
$\displaystyle k_{1}$ $\displaystyle=$
$\displaystyle-2(1+\Omega^{2}\tau^{2})(a+b\tau+a\Omega^{2}\tau^{2})$ (S42)
$\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle-2\tau^{2}(b\tau+a(2-2\Omega^{2}\tau^{2}))$ (S43) $\displaystyle
k_{3}$ $\displaystyle=$ $\displaystyle-2a\tau^{4}\,.$ (S44)

Interpreting the integrand in eq. (S38) as a sum of the three terms each
proportional to $k_{i}$, we see that all terms have the same poles, where the
term proportional to $k_{1}$ has an additional double pole at $\omega=0$. The
solutions of $\omega^{2}$ to the equation

$c_{1}+c_{2}\omega^{2}+c_{3}\omega^{4}+\tau^{4}\omega^{6}=\tau^{4}(\omega^{2}-\omega_{1}^{2})(\omega^{2}-\omega_{2}^{2})(\omega^{2}-\omega_{3}^{2})=0\,,$
(S45)

define the remaining poles, which we denote by $\pm\sqrt{\omega_{i}^{2}}$ and
which we compute by numerically solving eq. (S45). Next we use the partial
fraction decompositions

$\displaystyle\frac{1}{\omega^{2}\prod_{i=1}^{3}(\omega^{2}-\omega_{i}^{2})}$
$\displaystyle=$
$\displaystyle-\frac{1}{\omega^{2}\prod_{i=1}^{3}\omega_{i}^{2}}+\sum_{i=1}^{3}\frac{1}{\omega_{i}^{2}(\omega^{2}-\omega_{i}^{2})\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}$ (S46)
$\displaystyle\frac{1}{\prod_{i=1}^{3}(\omega^{2}-\omega_{i}^{2})}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{3}\frac{1}{(\omega^{2}-\omega_{i}^{2})\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}$ (S47)
$\displaystyle\frac{\omega^{2}}{\prod_{i=1}^{3}(\omega^{2}-\omega_{i}^{2})}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{3}\frac{\omega_{i}^{2}}{(\omega^{2}-\omega_{i}^{2})\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}$ (S48)

to rewrite the fraction of eq. (S38) as a sum of terms proportional to
$(\omega^{2}-\omega_{i}^{2})^{-1}$ and one term proportional to $\omega^{-2}$.
Using the solution of the integrals

$\displaystyle\int_{-\infty}^{\infty}\frac{e^{i\omega
t}-1}{\omega^{2}-\omega_{i}^{2}}d\omega$ $\displaystyle=$
$\displaystyle\frac{\pi(e^{-\sqrt{-\omega_{i}^{2}}t}-1)}{\sqrt{-\omega_{i}^{2}}}$
(S49) $\displaystyle\int_{-\infty}^{\infty}\frac{e^{i\omega
t}-1}{\omega^{2}}d\omega$ $\displaystyle=$ $\displaystyle-\pi t$ (S50)

for $t>0$ with the condition
$\text{Re}(\omega_{i}^{2})<0\lor\text{Im}(\omega_{i}^{2})\neq
0\land\text{Im}(\sqrt{\omega_{i}^{2}})\neq 0$, we can rewrite the integral of
the MSD eq. (S38) as

$C_{\textrm{MSD}}(t)=\frac{B}{\tau^{4}}\left(\frac{k_{1}t}{\omega_{1}^{2}\omega_{2}^{2}\omega_{3}^{2}}+\sum_{i=1}^{3}\frac{e^{-\sqrt{-\omega_{i}^{2}}t}-1}{\sqrt{-\omega_{i}^{2}}\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}\left[\frac{k_{1}}{\omega_{i}^{2}}+k_{2}+k_{3}\omega_{i}^{2}\right]\right)\,.$
(S51)

The VACF can be computed from the MSD by using eq. (23) as

$C_{vv}(t)=\frac{B}{2\tau^{4}}\left(\sum_{i=1}^{3}\frac{\sqrt{-\omega_{i}^{2}}e^{-\sqrt{-\omega_{i}^{2}}t}}{\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}\left[\frac{k_{1}}{\omega_{i}^{2}}+k_{2}+k_{3}\omega_{i}^{2}\right]\right)\,.$
(S52)

The calculation of the MSD for the friction kernel of eq. (26) proceeds in the
same way as shown for the kernel eq. (7). The response function then takes the
slightly different form

$\tilde{\chi}(\omega)=\left(-\omega^{2}+i\omega a+i\omega
b\left(\frac{i\omega+\frac{2}{\tau}}{\Omega^{2}+(i\omega+\frac{1}{\tau})^{2}}\right)\right)^{-1}\,.$
(S53)

The integral eq. (S35) resulting in the MSD can be written in the same form as
before shown in eq. (S38), where the constants $c_{1}$, $c_{2}$, $k_{1}$,
$k_{2}$ take different values of

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle(a+2b\tau+a\Omega^{2}\tau^{2})^{2}$ (S54) $\displaystyle c_{2}$
$\displaystyle=$
$\displaystyle(1+2(a^{2}-3b+\Omega^{2})\tau^{2}+(-2a^{2}\Omega^{2}+(b+\Omega^{2})^{2})\tau^{4})$
(S55) $\displaystyle k_{1}$ $\displaystyle=$
$\displaystyle-2(1+\Omega^{2}\tau^{2})(a+2b\tau+a\Omega^{2}\tau^{2})$ (S56)
$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle
4a\tau^{2}(\Omega^{2}\tau^{2}-1)\,.$ (S57)

Thus, the MSD of the GLE eq. (1) with the friction kernel eq. (26) can be
described by the same form eq. (S51) as for the friction kernel of eq. (7),
where only the constants $c_{i}$ and $k_{i}$ have slightly different values as
seen by comparing eqs. (S39)-(S44) to eqs. (S54)-(S57).

##  VII Localization noise fit

The fitting procedure including finite time step and localization noise
described in Sec. 4.7 of the main text does sometimes not converge to results
that agree well with the VACF of the data. In Fig. S3 we show two examples for
which the fit including the localization noise (Methods eq. (24)) does not
agree well with the VACF data. For 14 out of the 59 cells, the fits of the
VACF using eq. (25) do not converge to a stable set of parameters. Still the
direct fit of the friction kernel model eq. (7) agrees perfectly with the
extracted friction kernels, which is why we use this direct fit to extract
parameters, that are later used in the cluster analysis.

A typical size of the localization noise width resulting from the fits is
${\sigma_{\textrm{loc}}\approx 0.02\,\mu m}$. Considering the resolution of
the microscope of roughly ${\sim 0.5\,\mu m}$ and approximating the area of a
cell as roughly ${25\pi\,\mu m^{2}}$, one can estimate the number of pixels
per cell as $100\pi$. Thus, the error of the mean position of a cell evaluated
by all pixels is estimated by ${0.5\mu m/\sqrt{100\pi}\approx 0.03\,\mu m}$,
which is very close to the fitted localization noise width of
${\sigma_{\textrm{loc}}\approx 0.02\,\mu m}$. For a typical velocity around
$v=100\,\mu m/s$, the displacement during one time step $\Delta=0.002\,s$ is
$0.2\,\mu m$. Therefore, the localization noise accounts for $\sim 10\%$ of
the cell displacement. This relatively small localization noise explains why
the direct fit of the model eq. (7) to the extracted friction kernels works so
well and the inclusion of localization noise effects is not necessary for our
data set.

Figure S3: Two example cells for which the localization-noise fit to VACF does
not agree well with the data. The MSD $C_{\textrm{MSD}}(t)$ (eq. (6)), VACF
$C_{vv}(t)$ (eq. (3)) and friction kernel $\Gamma(t)$ extracted by eq. (15) of
a single (a)-(c) wobbler and (d)-(f) synchro are displayed as orange dots.
Blue crosses in the VACF figures represent the results from the localization
noise fit according to eq. (24) and the blue crosses in the MSD figures result
from inserting the fit parameters of the VACF into eq. (22). The black dashed
line in (c) and (f) is the direct fit of eq. (7) to the extracted kernels,
which is used to obtain parameters for the cell classification, see Methods
Sec. 4.6. The dashed black lines in the MSD figures result from inserting the
fitting result of the friction kernel (dashed lines in (c), (f) respectively)
and the mean squared velocity $B=C_{vv}(0)$ into the analytical solution of
the MSD eq. (S51).

##  VIII Cluster analysis in lower dimensions

Since the wobblers are found to exhibit higher average speeds than the
synchros, an intuitive and simple approach would be to distinguish the cells
solely by their mean squared velocity $B$. Applying the cluster analysis
described in Sec. 4.8 of the Methods in one dimension for $B$, results in a
assignment agreeing only to $69\,\%$ with the full classification result, see
Fig. S4a. Especially slow wobblers are wrongly assigned as synchrose, as
indicated in Fig. S4a by the black crosses.

As discussed in the main text, high-dimensional data sets are often analyzed
with principal component analysis (PCA), which determines the directions
explaining most of the data variance. Applying a PCA to the five dimensional
parameter space of the extracted friction kernel parameters eq. (7) and $B$
shown in Fig. 5, we find the first two components of the PCA to explain
$83\,\%$ of the total parameter variance. Here, we use the PCA tool
implemented in sklearn for python. The normalized vectors of the shown two PCA
components in the order ${(a,\hskip 2.84526ptb,\hskip 2.84526pt\tau,\hskip
2.84526pt\Omega,\hskip 2.84526ptB)}$ are given by ${(-0.49,\hskip
2.84526pt-0.29,\hskip 2.84526pt0.16,\hskip 2.84526pt-0.01,\hskip
2.84526pt0.81)}$ and ${(0.63,\hskip 2.84526pt0.41,\hskip 2.84526pt-0.5,\hskip
2.84526pt0.05,\hskip 2.84526pt0.43)}$, respectively. This indicates, that no
parameter alone can describe the complete variance and therefore no parameter
alone can explain all the differences between synchros and wobblers. The
parameters most important for the variance and discriminability of the single
CR cells are $a$, $b$ and $B$. However, the distinction into the two groups of
wobblers and synchros by our cluster analysis (Sec. 4.8 in the Methods) works
the best using the complete five dimensional set of parameters. Applying the
cluster analysis to the first two PCA components results in an accuracy of
$90\,\%$, as shown in Fig. S4b. Six wobblers are assigned to belong to the
cluster of synchros (defined by the classification of mondal_strong_2021  that
agrees with our cluster analysis of the complete parameter set Fig. 5) as they
lie in the transition area between the two clusters and are indicated by the
black crosses in Fig. S4b.

Figure S4: (a) Mean squared velocities $B$ for all cells; a cluster analysis
solely based on B leads to the wrong assignment of 18 wobblers as synchros,
indicated by black crosses. (b) Projection of the five dimensional extracted
parameters shown in Fig. 5 on the first two components of the PCA, which
explain $61\,\%$ and $22\,\%$ of the total parameter variance, respectively.
The black crosses indicate six wobbler cells that are wrongly assigned as
synchros when applying the cluster analysis only to the shown two PCA
components.

##  IX Comparing friction kernel expressions eq. (7) and eq. (26)

The friction kernel of eq. (7) in the main text used to extract parameters
from the data is very similar to the friction kernel eq. (26), which is
derived in Sec. XIII from a system of harmonically coupled particles. The term
$1/(\tau\Omega)$ in front of the sine in eq. (26) becomes small when the decay
time $\tau$ is longer than the oscillation period $1/\Omega$. Several
oscillations occur before the friction kernel has decayed for the CR data in
Fig. 3 in the main text and Fig. S5, thus, the term $1/(\tau\Omega)$ is indeed
small and the two friction kernel models have very similar shapes, as we show
in Fig. S5. Here, we insert the extracted parameters from the friction kernel
eq. (7) of two example cells into eq. (26), the deviations are seen to be very
small.

Figure S5: Friction kernel $\Gamma(t)$ extracted by eq. (15) of (a) a single
wobbler and (b) single synchro shown as green data points. The black line
shows the direct fit of eq. (7) to the data points and the red dotted line
results from inserting the fitted parameters from the black dashed line into
eq. (26) derived from the Hamiltonian eq. (27) in Sec. XIII. The agreement of
the black dashed line and the red dotted line means that in the parameter
range exhibited by the CR cells, the two kernels eq. (7) and eq. (26) describe
the data equally well.

##  X Exemplary non-equilibrium model describing CR motion

The effective friction kernel $\Gamma(t)$ given by eq. (26) leads to the
Fourier transform

$\displaystyle\tilde{\Gamma}_{+}(\omega)$
$\displaystyle=\int_{0}^{\infty}\left(2a\delta(t)+be^{-t/\tau}\left(\cos(\Omega
t)+\frac{1}{\tau\Omega}\sin(\Omega t)\right)\right)e^{-i\omega t}dt$
$\displaystyle=a+b\left(\frac{i\omega+\frac{2}{\tau}}{\Omega^{2}+(i\omega+\frac{1}{\tau})^{2}}\right)\,.$
(S58)

Inserting eq. (S58) into eq. (11) yields the Fourier-transformed VACF

$\begin{split}&\frac{\tilde{C}_{vv}(\omega)}{B}=\bigg{(}2a+2b\tau+4a\Omega^{2}\tau^{2}+2b\Omega^{2}\tau^{3}+2a\Omega^{4}\tau^{4}+4ib\tau^{2}\omega\\\
&+(4a\tau^{2}+2b\tau^{3}-4a\Omega^{2}\tau^{4})\omega^{2}+2a\tau^{4}\omega^{4}\bigg{)}/\\\
&\bigg{(}a+2b\tau+a\Omega^{2}\tau^{2})^{2}+(1+2(a^{2}-3b+\Omega^{2})\tau^{2}+(-2a^{2}\Omega^{2}+(b+\Omega^{2})^{2})\tau^{4})\omega^{2}\\\
&+\tau^{2}(2+(a^{2}-2(b+\Omega^{2}))\tau^{2})\omega^{4}+\tau^{4}\omega^{6}\bigg{)}\\\
\end{split}\,.$ (S59)

In order to find a non-equilibrium model, that leads to the same VACF as the
kernel in eq. (26), we assume a specific form of the velocity friction kernel

$\Gamma_{v}(t)=2a_{v}\delta(t)$ (S60)

which leads to the Fourier transform ${\tilde{\Gamma}_{v}^{+}(\omega)=a_{v}}$.
We insert this together with the result from eq. (S59) into eq. (12) and solve
for $\tilde{\Gamma}_{R}(\omega)$. By Fourier back transformation the result
can be written in the form

$\Gamma_{R}(t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{i\omega
t}\frac{k_{0}+k_{1}\omega+k_{2}\omega^{2}+k_{3}\omega^{3}+k_{4}\omega^{4}+k_{6}\omega^{6}}{c_{0}+c_{1}\omega^{2}+c_{2}\omega^{4}+c_{3}\omega^{6}}$
(S61)

with the constants $c_{i}$ and $k_{i}$ being

$\displaystyle c_{0}$ $\displaystyle=$
$\displaystyle(a+2b\tau+a\Omega^{2}\tau^{2})^{2}$ $\displaystyle c_{1}$
$\displaystyle=$ $\displaystyle
1+2(a^{2}-3b+\Omega^{2})\tau^{2}+(-2a^{2}\Omega^{2}+(b+\Omega^{2})^{2})\tau^{4}$
$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle\tau^{2}(2+(a^{2}-2(b+\Omega^{2}))\tau^{2})$ $\displaystyle
c_{3}$ $\displaystyle=$ $\displaystyle\tau^{4}$ $\displaystyle k_{0}$
$\displaystyle=$ $\displaystyle
2aa_{v}^{2}+2a_{v}^{2}b\tau+4aa_{v}^{2}\Omega^{2}\tau^{2}+2a_{v}^{2}b\Omega^{2}\tau^{3}+2aa_{v}^{2}\Omega^{4}\tau^{4}$
$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle 4a_{v}^{2}b\tau^{2}$
$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle
2(b\tau(1+a_{v}^{2}\tau^{2}+\Omega^{2}\tau^{2})+a((1+\Omega^{2}\tau^{2})^{2}+a_{v}^{2}(2\tau^{2}-2\Omega^{2}\tau^{4})))$
$\displaystyle k_{3}$ $\displaystyle=$ $\displaystyle 4b\tau^{2}$
$\displaystyle k_{4}$ $\displaystyle=$ $\displaystyle
2\tau^{2}(b\tau+a(2+a_{v}^{2}\tau^{2}-2\Omega^{2}\tau^{2}))$ $\displaystyle
k_{6}$ $\displaystyle=$ $\displaystyle 2a\tau^{4}\,.$

We solve this Fourier transform by performing the same steps as for the
calculation of the integral (S38) derived in Sec. VI, where $\omega_{i}^{2}$
are the roots of the denominator of eq. (S61) in $\omega^{2}$, which we obtain
numerically, and we additionally use

$\int_{-\infty}^{\infty}\frac{\omega e^{i\omega
t}}{\omega^{2}-\omega_{i}^{2}}d\omega=i\pi e^{-\sqrt{-\omega_{i}^{2}}t}\,.$
(S62)

The friction kernel in time domain is then retrieved as

$\begin{split}\Gamma_{R}(t)=2a\delta(t)+\frac{1}{2c_{3}}\sum_{i=1}^{3}\bigg{(}\frac{e^{-\sqrt{-\omega_{i}^{2}}t}}{\sqrt{-\omega_{i}^{2}}\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}\left[k_{0}+k_{2}\omega_{i}^{2}+k_{4}(\omega_{i}^{2})^{2}+k_{6}(\omega_{i}^{2})^{3}\right]\\\
-\frac{e^{-\sqrt{-\omega_{i}^{2}}t}}{\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}\left[k_{1}+k_{3}\omega_{i}^{2}\right]\bigg{)}.\end{split}\,$
(S63)

This exemplary mapping shows, that a simple Markovian friction kernel given by
eq. (S60) combined with a random force correlation function that contains
additional oscillating components,given by eq. (S63), lead to the same
correlation function as described by the GLE with the effective friction
kernel given by eq. (26). We have thus derived one possible non-equilibrium
model that describes CR cell motion.

Furthermore, this non-equilibrium model with $\Gamma_{R}(t)$ given by eq.
(S63) and $\Gamma_{v}(t)$ given by eq. (S60) corresponds to the coupled system
of differential equations

$\displaystyle\dot{x}(t)$ $\displaystyle=v(t)$ (S64) $\displaystyle\dot{v}(t)$
$\displaystyle=-a_{v}v(t)+\xi_{x}(t)+\sum_{i=1}^{3}F_{i}(t)$ (S65)
$\displaystyle\dot{F}_{i}(t)$
$\displaystyle=-\frac{1}{\tau_{i}}\left(F_{i}(t)-\xi^{R}_{i}(t)\right)\,,$
(S66)

with $\langle\xi_{x}(0)\xi_{x}(t)\rangle=2aB\delta(t)$ and
$\langle\xi^{R}_{i}(0)\xi^{R}_{j}(t)\rangle=2a^{R}_{i}B\delta_{ij}\delta(t)$.
Solving eq. (S66) for $F_{i}(t)$ leads to

$F_{i}(t)=-\frac{1}{\tau_{i}}\int_{0}^{t}e^{-(t-t^{\prime})/\tau_{i}}\xi^{R}_{i}(t^{\prime})dt^{\prime}\,,$
(S67)

which yields $\langle
F_{i}(0)F_{i}(t)\rangle=B\frac{a^{R}_{i}}{\tau_{i}}e^{-t/\tau_{i}}$. Now
defining the random force as $F_{R}(t)=\xi_{x}(t)+\sum_{i=1}^{3}F_{i}(t)$, we
can write eq. (S65) as the GLE eq. (1) with

$\Gamma_{R}(t)=2a\delta(t)+\sum_{i=1}^{3}\frac{a^{R}_{i}}{\tau_{i}}e^{-t/\tau_{i}}$
(S68)

and $\Gamma_{v}(t)=2a_{v}\delta(t)$. Thus, the system of eqs. (S64)-(S66) is
equivalent to the GLE eq. (1) with $\Gamma_{R}(t)$ given by eq. (S63) and
$\Gamma_{v}(t)$ given by eq. (S60), when the parameters are given by

$\displaystyle\tau_{i}$ $\displaystyle=\frac{1}{\sqrt{-\omega_{i}^{2}}}$ (S69)
$\displaystyle a^{R}_{i}$
$\displaystyle=\frac{1}{-2c_{3}\omega_{i}^{2}\prod_{j\neq
i}(\omega_{i}^{2}-\omega_{j}^{2})}\left[k_{0}+k_{2}\omega_{i}^{2}+k_{4}(\omega_{i}^{2})^{2}+k_{6}(\omega_{i}^{2})^{3}-\sqrt{-\omega_{i}^{2}}(k_{1}+k_{3}\omega_{i}^{2})\right]\,.$
(S70)

##  XI Information on the cell orientation is contained in the cell-center
trajectories

The mean direction of motion, exponentially weighted over the past trajectory,
is defined by

$\vec{V}(t)=\frac{1}{T}\int_{-\infty}^{t}e^{-(t-t^{\prime})/T}\vec{v}(t^{\prime})dt^{\prime}\hskip
8.53581pt,$ (S71)

where $T$ is the exponential decay time and $\vec{v}(t)=(v_{x}(t),v_{y}(t))$
is the two-dimensional velocity li2011dicty . This running-average velocity
vector $\vec{V}(t)$ points in the direction in which a cell has moved in the
past, the running-average position follows as

$\vec{R}(t)=\int_{0}^{t}\vec{V}(t^{\prime})dt^{\prime}\,.$ (S72)

In Fig. S6a we show the trajectory of a synchro and its running average
position $\vec{R}(t)$ for $T=0.1\,s$, $\vec{R}(t)$ resolves the mean direction
of motion but not the back-and-forth motion on time scales below $T$. From the
circular path of the trajectory it follows that the mean direction of motion
changes, this is reflected in Fig. S6b by the change of the angle of
$\vec{R}(t)$ with the x-axis (red data).

The orientation of cells is determined from the videos by fitting an ellipsoid
to the cell perimeter and given by the main axis of the ellipsoid. Since CR
cells move with their flagella in front, the orientation vector is defined to
point towards the side where the flagella are anchored.

The cell orientation vector coincides very well with the direction of the
running average velocity $\vec{V}(t)$, as shown in Fig. S6b. Therefore, the
information on the orientation of the CR cells is included in the cell
trajectory $(v_{x}(t),v_{y}(t))$ via eq. (S71). Our description of the CR
motion by the GLE eq. (1) captures the correct behavior of the cell position
and velocity over time, which is demonstrated by Figs. 4a-f. Since the cell
orientational information is included in the trajectory, we conclude that the
GLE model also correctly describes the CR orientational dynamics.

Figure S6: (a) Trajectory of a synchro shown in black, the running average
position defined by eqs. (S71), (S72) with a decay time $T=0.1\,s$ is shown in
red. The inset in the middle shows a zoom into the first $0.2\,s$ of the
trajectory. (b) The angle between the x-axis and the cell orientation directly
extracted from video data is shown for the trajectory in (a) as black crosses.
The angle between the x-axis and the running-average velocity $\vec{V}(t)$
determined by eq. (S71) is shown in red.

##  XII Effects of smoothing trajectory data

Since the localization noise modifies the VACF and thereby the friction kernel
that is extracted from the data, it is important to smooth the data at a level
that reduces the noise without loosing detail. For every data set, there is an
optimal level of smoothing.

A natural way of smoothing data is by sequentially averaging over neighboring
data points, which conserves the discretization time step. Here, we use an
iterative averaging of the position data over two consecutive points according
to

$x_{i}^{n+1}=\frac{x_{i}^{n}+x_{i+1}^{n}}{2}\,,$ (S73)

with $n$ being the number of smoothing iterations. We show the effect of the
smoothing level on the VACF and the localization noise fit in Fig. S7 and the
effect on the friction kernel in Fig. S8. For high smoothing iterations, the
cell speed is underestimated, which is reflected by a decreasing mean squared
velocity $B$ for increasing $n$ in Fig. S7. At the same time, the fitted
localization noise strength $\sigma_{\textrm{loc}}$ decreases with increasing
$n$. This is reflected by the smaller difference between the first two data
points of the VACF.

Figure S7: The VACF $C_{vv}(t)$ of a wobbler (left) and a synchro (right) is
shown for different smoothing iterations $n$ described by eq. (S73), where
$n=0$ denotes the original non-smoothed data. The mean-squared velocities $B$
in units of $\mu m^{2}/s^{2}$ and the fitted localization noise width
$\sigma_{\textrm{loc}}$ in units of $\mu m$ are both given in the color of the
according smoothing iteration. The fitting result of the VACF including
localization noise, described in Sec. VII, is shown as solid lines for $n=0$
and $n=9$ in the respective color. Figure S8: The friction kernel $\Gamma(t)$
extracted according to eq. (15) from data at different smoothing iterations
described by eq. (S73) is given for a wobbler on the left and for a synchro on
the right. The solid lines connecting the data points are guides to the eye.

The estimate of the localization noise strength for high smoothing iterations
$n$ is not reliable anymore because it becomes very small, this explains why
it is not decreasing monotonically with growing $n$, see Fig. S7. Since the
first smoothing iteration $n=1$ simply adds two consecutive points according
to eq. (S73), one in principle expects the localization noise width
$\sigma_{\textrm{loc}}$ to be halved compared to the non-smoothed data. This
indeed one can see for the synchro shown in Fig. S7.

Since the smoothing of the data decreases the amplitude of the VACF, as seen
in Fig. S7, it consequently decreases the amplitudes of the friction kernel
parameters $a$ and $b$ and therefore at the same time it reduces the kernel
frequency $\Omega$, as $\Omega$ exhibits a complex relation with the frequency
of the VACF, which is explained in detail in Sec. V. Moreover, we find the dip
in the friction kernel at the first time step to decrease for higher smoothing
iterations shown in Fig. (S8), which indicates that the dip in the friction
kernel originates from the localization noise. This is explained by eq. (15),
since the first point of the VACF is overestimated and the second point of the
VACF is underestimated due to the localization noise, which then propagates
into the friction kernel and leads to a dip at short times.

The VACF of the smoothed data is less influenced by the localization noise,
nevertheless, with every smoothing iteration, the deviation of the smoothed
VACF from the non-smoothed VACF becomes larger, as shown in Fig. S7. Thus, we
choose $n=1$ smoothing iterations for our data as a compromise between
minimizing the effect of noise and keeping a good resolution that accurately
describes the actual cell velocities.

##  XIII Markovian embedding of the friction kernel eq. (26)

Here we derive the equivalency of the system of coupled equations of motion

$\displaystyle\dot{x}(t)$ $\displaystyle=$ $\displaystyle v(t)$ (S74)
$\displaystyle m\dot{v}(t)$ $\displaystyle=$
$\displaystyle-\gamma_{x}v+bm(y(t)-x(t))+F_{Rx}(t)$ (S75)
$\displaystyle\dot{y}(t)$ $\displaystyle=$ $\displaystyle v_{y}(t)$ (S76)
$\displaystyle m_{y}\dot{v}_{y}(t)$ $\displaystyle=$
$\displaystyle-\gamma_{y}v_{y}+bm(x(t)-y(t))+F_{Ry}(t)\,,$ (S77)

describing two harmonically coupled masses in a heat bath based on the
Hamiltonian eq. (27), to the GLE eq. (1) with the effective friction kernel of
eq. (26) brunig2022time . Here, the random forces have zero mean $\langle
F_{i}(t)\rangle=0$ and are delta correlated $\langle
F_{i}(0)F_{j}(t)\rangle=2\gamma_{i}k_{B}T\delta_{ij}\delta(t)$. We start by
finding the solution for $\vec{z}(t)=(y,\dot{y})^{T}$ to the subsystem eqs.
(S76) and (S77), which is the solution to

$\dot{\vec{z}}(t)=A\vec{z}(t)+\begin{pmatrix}0\\\
F_{Ry}(t)/m_{y}+\frac{bm}{m_{y}}x(t)\end{pmatrix}$ (S78)

with

$A=\begin{pmatrix}0&1\\\
-\frac{bm}{m_{y}}&-\frac{\gamma_{y}}{m_{y}}\end{pmatrix}$ (S79)

and consequently

$A^{-1}=\begin{pmatrix}-\frac{\gamma_{y}}{bm}&-\frac{m_{y}}{bm}\\\
1&0\end{pmatrix}\,.$ (S80)

One finds the solution to eq. (S78) as

$\vec{z}(t)=e^{A[t-t_{0}]}\vec{z}(t_{0})+\int_{t_{0}}^{t}dt^{\prime}e^{A[t-t^{\prime}]}\begin{pmatrix}0\\\
F_{Ry}(t^{\prime})/m_{y}+\frac{bm}{m_{y}}x(t^{\prime})\end{pmatrix}\,.$ (S81)

Performing a partial integration on the integral in eq. (S81) for the term
including $x(t^{\prime})$ leads to

$\begin{split}\vec{z}(t)=e^{A[t-t_{0}]}\vec{z}(t_{0})+\int_{t_{0}}^{t}dt^{\prime}A^{-1}e^{A[t-t^{\prime}]}\begin{pmatrix}0\\\
\frac{bm}{m_{y}}v(t^{\prime})\end{pmatrix}\\\
-\left[A^{-1}e^{A[t-t^{\prime}]}\begin{pmatrix}0\\\
\frac{bm}{m_{y}}x(t^{\prime})\end{pmatrix}\right]_{t_{0}}^{t}\\\
+\int_{t_{0}}^{t}dt^{\prime}e^{A[t-t^{\prime}]}\begin{pmatrix}0\\\
F_{Ry}(t^{\prime})/m_{y}\end{pmatrix},\end{split}\,,$ (S82)

where the matrix exponential $e^{At}$ can be expressed in terms of the
Eigenvalues of $A$

$\lambda_{1,2}=-\frac{\gamma_{y}}{2m_{y}}\pm\omega_{0}$ (S83)

with

$\omega_{0}=\sqrt{\left(\frac{\gamma_{y}}{2m_{y}}\right)^{2}-\frac{bm}{m_{y}}}$
(S84)

as

$\displaystyle e^{At}$
$\displaystyle=\frac{1}{\lambda_{2}-\lambda_{1}}\begin{pmatrix}\lambda_{2}e^{\lambda_{1}t}-\lambda_{1}e^{\lambda_{2}t}&e^{\lambda_{2}t}-e^{\lambda_{1}t}\\\
\lambda_{1}\lambda_{2}(e^{\lambda_{1}t}-e^{\lambda_{2}t})&\lambda_{2}e^{\lambda_{2}t}-\lambda_{1}e^{\lambda_{1}t}\end{pmatrix}$
(S85)
$\displaystyle=e^{-t\frac{\gamma_{y}}{2m_{y}}}\begin{pmatrix}\cosh(\omega_{0}t)+\frac{\sinh(\omega_{0}t)\gamma_{y}}{2m_{y}\omega_{0}}&\sinh(\omega_{0}t)/\omega_{0}\\\
-\frac{\sinh(\omega_{0}t)bm}{\omega_{0}m_{y}}&\cosh(\omega_{0}t)-\frac{\sinh(\omega_{0}t)\gamma_{y}}{2m_{y}\omega_{0}}\end{pmatrix}$

consequently leading to

$A^{-1}e^{At}=e^{-t\frac{\gamma_{y}}{2m_{y}}}\begin{pmatrix}-\frac{\cosh(\omega_{0}t)\gamma_{y}}{bm}+\frac{\sinh(\omega_{0}t)}{\omega_{0}}(1-\frac{\gamma_{y}^{2}}{2bmm_{y}})&-\frac{\sinh(\omega_{0}t)\gamma_{y}}{2bm\omega_{0}}-\frac{\cosh(\omega_{0}t)m_{y}}{bm}\\\
\cosh(\omega_{0}t)+\frac{\sinh(\omega_{0}t)\gamma_{y}}{2m_{y}\omega_{0}}&\frac{\sinh(\omega_{0}t)}{\omega_{0}}\end{pmatrix}\,.$
(S86)

The solution for $y(t)$ can now be found as the first element of the vector
$\vec{z}(t)$ from eqs. (S82)-(S86) when assuming a stationary system, i.e.
$t_{0}\rightarrow-\infty$ as

$\begin{split}y(t)&=x(t)-\int_{-\infty}^{t}dt^{\prime}v(t^{\prime})e^{-\frac{(t-t^{\prime})\gamma_{y}}{2m_{y}}}\left(\cosh(\omega_{0}(t-t^{\prime}))+\frac{\gamma_{y}}{2m_{y}\omega_{0}}\sinh(\omega_{0}(t-t^{\prime}))\right)\\\
&+\int_{-\infty}^{t}dt^{\prime}F_{Ry}(t^{\prime})e^{-\frac{(t-t^{\prime})\gamma_{y}}{2m_{y}}}\frac{\sinh(\omega_{0}(t-t^{\prime}))}{m_{y}\omega_{0}}.\end{split}$
(S87)

Inserting eq. (S87) into eq. (S75), one recovers the GLE eq. (1) for $x(t)$

$m\ddot{x}(t)=-\int_{-\infty}^{t}dt^{\prime}\bar{\Gamma}(t-t^{\prime})+\bar{F}_{R}(t),$
(S88)

where the bars indicate multiplication by the mass
$\Gamma(t)=\bar{\Gamma}(t)/m$, $F_{R}=\bar{F}_{R}/m$ and the friction kernel
takes the oscillatory form

$\Gamma(t)=\frac{\gamma_{x}}{m}\delta(t)+be^{-t\frac{\gamma_{y}}{2m_{y}}}\left(\cosh(\omega_{0}t)+\frac{\gamma_{y}}{2m_{y}\omega_{0}}\sinh(\omega_{0}t)\right)$
(S89)

and the random force is explicitly given by

$F_{R}(t)=\bar{F}_{Rx}/m+b\int_{-\infty}^{t}dt^{\prime}\bar{F}_{Ry}(t^{\prime})e^{-\frac{(t-t^{\prime})\gamma_{y}}{2m_{y}}}\frac{\sinh(\omega_{0}(t-t^{\prime}))}{m_{y}\omega_{0}}\,.$
(S90)

Computing the random force correlation, one finds

$\langle F_{R}(0)F_{R}(t)\rangle=B\Gamma(t)\,.$ (S91)

The friction kernel of eq. (26) in the Methods is equivalent to the friction
kernel of eq. (S89) with the parameters given by

$\displaystyle a$ $\displaystyle=$ $\displaystyle\frac{\gamma_{x}}{m}$
$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\frac{2m_{y}}{\gamma_{y}}$
$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle i\omega_{0}\,.$

## References

  *   * (1) Netz, R.R.: Multi-point distribution for Gaussian non-equilibrium non-Markovian observables (2023). https://arxiv.org/abs/2310.08886
  * (2) Mondal, D., Prabhune, A.G., Ramaswamy, S., Sharma, P.: Strong confinement of active microalgae leads to inversion of vortex flow and enhanced mixing. eLife 10, 67663 (2021). https://doi.org/10.7554/eLife.67663. Publisher: eLife Sciences Publications, Ltd. Accessed 2022-04-27
  * (3) Li, L., Cox, E.C., Flyvbjerg, H.: ‘dicty dynamics’: Dictyostelium motility as persistent random motion. Physical biology 8(4), 046006 (2011)
  * (4) Brünig, F.N., Geburtig, O., Canal, A.v., Kappler, J., Netz, R.R.: Time-dependent friction effects on vibrational infrared frequencies and line shapes of liquid water. The Journal of Physical Chemistry B 126(7), 1579–1589 (2022)


Fix $w\in\mathbb{Y}$ and note that the events $\\{y^{i}_{t+1}=w\\}$ and
$\\{\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\\}$ are equivalent. We claim that

$\displaystyle{\rm
Pr}^{\bm{\pi}}_{\nu}\left[\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right]={\rm
Pr}^{\bm{\pi}}_{\nu}\left[\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}\in\varphi^{i\;-1}(\varphi^{i}(\textbf{s})),u^{i}_{t}=a^{i}\right]$
(6)

holds for any $\textbf{s}\in\textbf{X}$ and any $a^{i}\in\mathbb{U}$.

To see this, let
$\tilde{\textbf{s}}\in\varphi^{i\;-1}(\varphi^{i}(\textbf{s}))$; that is,
$s^{i}=\tilde{s}^{i}$ and
$\mu(\cdot|\tilde{\textbf{s}})=\mu(\cdot|\textbf{s})$. Note that there exists
a permutation $\sigma$ on $\mathcal{N}$ such that $\sigma(i)=i$ and
$\tilde{s}^{j}=s^{\sigma(j)}$ for all $j\in\mathcal{N}\setminus\\{i\\}$. (In
words, player $i$’s part remains unchanged while the local states of the
remaining agents are permuted.)

We overload notation so that the permutation $\sigma$ also acts on X and U:
for any $\textbf{a}\in\textbf{U}$, let $\sigma(\textbf{a})\in\textbf{U}$ be
the joint action such that $\sigma(a)^{j}=a^{\sigma(j)}$ for all
$j\in\mathcal{N}$; the permutation acts on X in an analogous manner. With this
overloaded notation, we have that $\tilde{\textbf{s}}=\sigma(\textbf{s})$.

Since the policies $\pi^{j}$ and $\pi^{\sigma(j)}$ are mean-field symmetric,
for any $\textbf{a}\in\textbf{U}$ and $\textbf{s}\in\textbf{X}$, we have

${\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{u}_{t}=\textbf{a}\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right)={\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{u}_{t}=\sigma(\textbf{a})\middle|\textbf{x}_{t}=\sigma(\textbf{s}),u^{i}_{t}=a^{i}\right).$
(7)

Then, since the local state variables $\\{x^{j}_{t+1}\\}_{j\in\mathcal{N}}$
are conditionally independent given $\textbf{x}_{t}$ and $\textbf{u}_{t}$, the
global state transition function has a product structure, which is invariant
under permutations of its terms:

${\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}=\textbf{s}^{\prime}\middle|\textbf{x}_{t}=\textbf{s},\textbf{u}_{t}=\textbf{a}\right)=\prod_{j\in\mathcal{N}}P_{{\rm
loc}}\left(s^{\prime j}\middle|s^{j},\mu(\cdot|\textbf{s}),a^{j}\right).$

Using the law of total probability, conditioning on $\textbf{u}_{t}$, and
recalling that $\tilde{\textbf{s}}=\sigma(\textbf{s})$, it follows that for
any $\textbf{s}^{\prime}\in\textbf{X}$,

$\displaystyle{\rm
Pr}^{\bm{\pi}}_{\nu}\left[\textbf{x}_{t+1}=\textbf{s}^{\prime}\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right]={\rm
Pr}^{\bm{\pi}}_{\nu}\left[\textbf{x}_{t+1}=\sigma(\textbf{s}^{\prime})\middle|\textbf{x}_{t}=\tilde{\textbf{s}},u^{i}_{t}=a^{i}\right].$
(8)

Noting that
$\mu(\cdot|\textbf{s}^{\prime})=\mu(\cdot|\sigma(\textbf{s}^{\prime}))$, we
have that
$\textbf{s}^{\prime}\in\varphi^{i\;-1}(w)\iff\sigma(\textbf{s}^{\prime})\in\varphi^{i\;-1}(w)$.
Thus,

$\varphi^{i\;-1}(w)=\bigcup_{\textbf{s}^{\prime}\in\varphi^{i\;-1}(w)}\\{\sigma(\textbf{s}^{\prime})\\}.$

It follows that

$\displaystyle{\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right)=$
$\displaystyle\sum_{\textbf{s}^{\prime}\in\varphi^{i\;-1}(w)}{\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}=\textbf{s}^{\prime}\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right)$
$\displaystyle=$
$\displaystyle\sum_{\textbf{s}^{\prime}\in\varphi^{i\;-1}(w)}{\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}=\sigma(\textbf{s}^{\prime})\middle|\textbf{x}_{t}=\sigma(\textbf{s}),u^{i}_{t}=a^{i}\right)$
$\displaystyle=$ $\displaystyle{\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}=\sigma(\textbf{s}),u^{i}_{t}=a^{i}\right).$

Since $\tilde{\textbf{s}}\in\varphi^{i\;-1}(\varphi^{i}(\textbf{s}))$ was
arbitrary, we see that

$\displaystyle{\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right)={\rm
Pr}^{\bm{\pi}}_{\nu}\left(\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}=\bar{\textbf{s}},u^{i}_{t}=a^{i}\right),$
(9)

for any $\bar{\textbf{s}}\in\varphi^{i\;-1}(\varphi^{i}(\textbf{s}))$. We
conclude that (6) holds by applying iterated expectations to its right-hand
side, conditioning on $\\{\textbf{x}_{t},u^{i}_{t}\\}$, and using (9) to
simplify the resulting summation.

Finally, since the event
$\\{\textbf{x}_{t}\in\varphi^{i\;-1}(\varphi^{i}(\textbf{s}))\\}$ is
equivalent to the event $\\{y^{i}_{t}=\varphi^{i}(\textbf{s})\\}$ and
$\\{\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\\}$ is equivalent to
$\\{y^{i}_{t+1}=w\\}$, we have that

$\displaystyle{\rm
Pr}^{\bm{\pi}}\left[\textbf{x}_{t+1}\in\varphi^{i\;-1}(w)\middle|\textbf{x}_{t}=\textbf{s},u^{i}_{t}=a^{i}\right]={\rm
Pr}^{\bm{\pi}}\left[y^{i}_{t+1}=w\middle|y^{i}_{t}=\varphi(\textbf{s}),u^{i}_{t}=a^{i}\right],$

for any $w\in\mathbb{Y}$, $\textbf{s}\in\textbf{X}$, and $a^{i}\in\mathbb{U}$.
The result follows. ∎

###  -B Proofs of Results in Section V

###### Proof:

Let $i,j\in\mathcal{N}$. For each $\textbf{s}\in\textbf{X}$, we define ${\rm
swap}_{ij}(\textbf{s})\in\textbf{X}$ to be the global state such that (i)
${\rm swap}_{ij}(\textbf{s})^{p}=s^{p}$ for each
$p\in\mathcal{N}\setminus\\{i,j\\}$, and (ii) we have ${\rm
swap}_{ij}(\textbf{s})^{j}=s^{i}$ and ${\rm swap}_{ij}(\textbf{s})^{i}=s^{j}$.

If the joint policy $\bm{\pi}$ is not soft, then by our definition of
$\mathcal{V}^{*}$ and $\mathcal{W}^{*}$, we have that $\pi^{i}\notin{\rm
Subj\text{-}BR}^{i}_{\epsilon}(\bm{\pi}^{-i},\mathcal{V}^{*},\mathcal{W}^{*})$
and $\pi^{j}\notin{\rm
Subj\text{-}BR}^{j}_{\epsilon}(\bm{\pi}^{-j},\mathcal{V}^{*},\mathcal{W}^{*})$.
We thus focus on joint policies $\bm{\pi}$ that are soft.

By Theorem 2, we have that
$V^{*i}_{\bm{\pi}}(\textbf{s})=J^{i}_{\bm{\pi}}(\textbf{s})$ for each
$\textbf{s}\in\textbf{X}$, $W^{*i}_{\bm{\pi}}=Q^{*i}_{\bm{\pi}^{-i}}$,
$V^{*j}_{\bm{\pi}}(\textbf{s})=J^{j}_{\bm{\pi}}(\textbf{s})$ for each
$\textbf{s}\in\textbf{X}$, and $W^{*j}_{\bm{\pi}}=Q^{*j}_{\bm{\pi}^{-j}}$.
Thus, we have that $\pi^{i}\in{\rm
Subj\text{-}BR}^{i}_{\epsilon}(\bm{\pi}^{-i},\mathcal{V}^{*},\mathcal{W}^{*})$
if and only if $\pi^{i}\in{\rm BR}^{i}_{\epsilon}(\bm{\pi}^{-i})$, and
analogously for $j$.

Toward obtaining a contradiction, assume that $\pi^{i}\in{\rm
Subj\text{-}BR}^{i}_{\epsilon}(\bm{\pi}^{-i},\mathcal{V}^{*},\mathcal{W}^{*})$
but $\pi^{j}\notin{\rm
Subj\text{-}BR}^{j}_{\epsilon}(\bm{\pi}^{-j},\mathcal{V}^{*},\mathcal{W}^{*})$.
Equivalently, assume that $\pi^{i}\in{\rm BR}^{i}_{\epsilon}(\pi^{-i})$ while
$\pi^{j}\in{\rm BR}^{i}_{\epsilon}(\pi^{-j})$. That is, we have that

$J^{i}_{\bm{\pi}}(\textbf{s})\leq\inf_{\tilde{\pi}^{i}\in\Gamma^{i}_{S}}J^{i}_{(\tilde{\pi}^{i},\bm{\pi}^{-i})}(\textbf{s})+\epsilon,\quad\forall\textbf{s}\in\textbf{X},$

while there exists $\textbf{s}^{*}\in\textbf{X}$ such that
$J^{j}_{\bm{\pi}}(\textbf{s}^{*})\geq\inf_{\tilde{\pi}^{j}\in\Gamma^{j}_{S}}J^{j}_{(\tilde{\pi}^{j},\bm{\pi}^{-j})}(\textbf{s}^{*})+\epsilon.$

Observe that if $\bm{\pi}\in\bm{\Gamma}_{S}$ is of the mean-field type and
$\tilde{\pi}^{i}\in\Gamma^{i}_{S}$ is any stationary policy for player $i$
(not necessarily of the mean-field type), then defining a policy
$\tilde{\pi}^{j}\in\Gamma^{j}_{S}$ state-wise by
$\tilde{\pi}^{j}(\cdot|\textbf{s})=\tilde{\pi}^{i}(\cdot|{\rm
swap}_{ij}(\textbf{s}))$ for all s, we have that

$J^{i}_{(\tilde{\pi}^{i},\bm{\pi}^{-i})}(\textbf{s})=J^{j}_{(\tilde{\pi}^{j},\bm{\pi}^{-j})}({\rm
swap}_{ij}(\textbf{s})),\quad\forall\textbf{s}\in\textbf{X}.$

It follows that, first, we have
$J^{i}_{\bm{\pi}}(\textbf{s})=J^{j}_{\bm{\pi}}({\rm swap}_{ij}(\textbf{s}))$
for each $\textbf{s}\in\textbf{X}$ and furthermore

$\inf_{\tilde{\pi}^{i}\in\Gamma^{i}_{S}}J^{i}_{(\tilde{\pi}^{i},\bm{\pi}^{-i})}(\textbf{s})=\inf_{\tilde{\pi}^{j}\in\Gamma^{j}_{S}}J^{j}_{(\tilde{\pi}^{j},\bm{\pi}^{-j})}({\rm
swap}_{ij}(\textbf{s})),\quad\forall\textbf{s}\in\textbf{X}.$

In particular, this holds for ${\rm swap}_{ij}(\textbf{s}^{*})$, which yields

$\displaystyle
J^{j}_{\bm{\pi}}(\textbf{s}^{*})>\inf_{\tilde{\pi}^{j}\in\Gamma^{j}_{S}}J^{j}_{(\tilde{\pi}^{j},\bm{\pi}^{-j})}(\textbf{s}^{*})+\epsilon=\inf_{\tilde{\pi}^{i}\in\Gamma^{i}_{S}}J^{i}_{(\tilde{\pi}^{i},\bm{\pi}^{-i})}({\rm
swap}_{ij}(\textbf{s}^{*}))+\epsilon\geq J^{i}_{\bm{\pi}}({\rm
swap}_{ij}(\textbf{s}^{*}))=J^{j}_{\bm{\pi}}(\textbf{s}^{*}),$

a contradiction, which completes the proof. ∎

###### Proof:

To prove this result, we focus on the learned value functions
$V^{*i}_{\bm{\pi}},W^{*i}_{\bm{\pi}}$, $V^{*j}_{\bm{\pi}}$, and
$W^{*j}_{\bm{\pi}}$. We argue that these are equal for players $i$ and $j$
whenever $\pi^{i}$ and $\pi^{j}$ are mean-field symmetric.

To rule out trivial cases, we note that if the policy $\bm{\pi}$ is not soft,
then we have $\pi^{i}\notin{\rm
Subj\text{-}BR}^{i}_{\epsilon}(\bm{\pi}^{-i},\mathcal{V}^{*},\mathcal{W}^{*})$
and also $\pi^{j}\notin{\rm
Subj\text{-}BR}^{j}_{\epsilon}(\bm{\pi}^{-j},\mathcal{V}^{*},\mathcal{W}^{*})$.
These hold by our definition of the naively learned subjective function family
$(\mathcal{V}^{*},\mathcal{W}^{*})$. We thus focus on the other case, where
$\bm{\pi}$ is soft.

In Section IV, we observed that, under Assumptions 2 and 5, the learned values
$V^{*i}_{\bm{\pi}}$ and $W^{*i}_{\bm{\pi}}$ were in fact the state value and
action value functions, respectively, for an MDP with state space
$\\{\varphi^{i}(\textbf{s}):\textbf{s}\in\textbf{X}\\}\subseteq\mathbb{Y}$. An
analogous statement holds for $V^{*j}_{\bm{\pi}}$ and $W^{*j}_{\bm{\pi}}$ and
player $j$.

We will argue that the (approximate belief) MDPs on
$\\{\varphi^{i}(\textbf{s}):\textbf{s}\in\textbf{X}\\}$ and on
$\\{\varphi^{j}(\textbf{s}):\textbf{s}\in\textbf{X}\\}$ are equivalent. Then,
since player $i$ and $j$’s policies are mean-field symmetric, they correspond
to the same stationary policy for this (approximate belief) MDP, and are
therefore both either $\epsilon$-optimal for that MDP or not; the result will
follow.

To see that the approximate belief MDP facing player $i$ is equivalent to that
facing player $j$, we observe that the construction of [49] depends only on
the stage cost function, the observation channel, and the unique invariant
distribution on the underlying state space. The stage cost and observation
channel are symmetric and shared by players $i$ and $j$, so it suffices to
show that the unique invariant distribution on X, say $\nu_{\bm{\pi}}$, is
symmetric in the following sense:

$\nu_{\bm{\pi}}({\rm
swap}_{ij}(\textbf{s}))=\nu_{\bm{\pi}}(\textbf{s}),\quad\forall\textbf{s}\in\textbf{X}.$

To see this, note that by Assumption 5, $\\{\textbf{x}_{t}\\}_{t\geq 0}$ is an
irreducible, aperiodic Markov chain on X under the soft policy $\bm{\pi}$.
Thus, we have for any $\tilde{\nu}\in\mathcal{P}(\textbf{X})$.

${\rm
Pr}^{\bm{\pi}}_{\tilde{\nu}}(\textbf{x}_{t}\in\cdot)\to\nu_{\bm{\pi}}(\cdot),$

as $t\to\infty$, where $\nu_{\bm{\pi}}$ is the (unique) invariant measure on X
induced by $\bm{\pi}$. In particular, putting the initial measure
$\tilde{\nu}=\text{Unif}(\textbf{X})$, we have that for each $t\geq 0$ and
each $\textbf{s}\in\textbf{X}$,

${\rm Pr}^{\bm{\pi}}_{\tilde{\nu}}(\textbf{x}_{t}=\textbf{s})={\rm
Pr}^{\bm{\pi}}_{\tilde{\nu}}(\textbf{x}_{t}={\rm swap}_{ij}(\textbf{s})).$

It follows that $\nu_{\bm{\pi}}(\textbf{s}):=\lim_{t\to\infty}{\rm
Pr}^{\bm{\pi}}_{\tilde{\nu}}(\textbf{x}_{t}=\textbf{s})=\lim_{t\to\infty}{\rm
Pr}^{\bm{\pi}}_{\tilde{\nu}}(\textbf{x}_{t}={\rm
swap}_{ij}(\textbf{s}))=:\nu_{\bm{\pi}}({\rm swap}_{ij}(\textbf{s}))$. From
this, it follows that the approximate belief MDP on $\mathbb{Y}$ faced by
player $i$ is the same as the approximate belief MDP on $\mathbb{Y}$ faced by
player $j$. ∎

###  -C Approximation Results on the Sequences of Learning Iterates

Remark: The contents of this and the next section closely resembles that of
[46, Appendix A]. The proof technique used here parallels the proof technique
of [46, Theorem 5.1].

In the next section, we prove that Algorithm 3 leads to the convergence of
joint policies as described in Theorems 6, 7, and 8. Since the evolution of
the policy process $\\{\bm{\pi}_{k}\\}_{k\geq 0}$ depends on the evolution of
the learning iterates $\\{\widehat{J}^{i}_{t},\widehat{Q}^{i}_{t}\\}_{t\geq
0,i\in\mathcal{N}}$, we begin by studying the convergence behaviour of these
iterates. We argue that when parameters are suitably selected, then these
learning iterates sampled at the end of each exploration phase will closely
approximate the subjective functions for that exploration phase, and
consequently the policy process of Algorithm 3 approximates the policy process
of the Markov chain resulting from Algorithm 2.

We note that when each agent $i\in\mathcal{N}$ uses Algorithm 3, it is
actually following a particular randomized, non-stationary policy. When all
agents use Algorithm 3 and the initial measure is
$\nu\in\mathcal{P}(\textbf{X})$, we use ${\rm Pr}_{\nu}$ (with no policy index
in the superscript) to denote the resulting probability measure on
trajectories of states and actions. For all other policies
$\tilde{\bm{\pi}}\in\bm{\Gamma}$, we continue to use ${\rm
Pr}^{\tilde{\bm{\pi}}}_{\nu}$, as before.

The policy process $\\{\bm{\pi}_{k}\\}_{k\geq 0}$ depends on the sequences
$\\{\widehat{J}^{i}_{t},\widehat{Q}^{i}_{t}\\}_{t\geq 0,i\in\mathcal{N}}$ only
through these sequences sampled at the end of exploration phases; that is, the
iterate sequences are relevant to the updating of policies only along the
subsequence of times $\\{t_{k}\\}_{k\geq 0}$. Recall that we used
$\\{\bar{Q}^{i}_{t}\\}_{t\geq 0}$ and $\\{\bar{J}^{i}_{t}\\}_{t\geq 0}$ to
denote the naively learned stochastic iterates obtained when player
$i\in\mathcal{N}$ employed Algorithm 1 and all players followed a soft
stationary policy. We now analyze the sequences
$\\{\widehat{Q}^{i}_{t}\\}_{t\geq 0}$ and $\\{\widehat{J}^{i}_{t}\\}_{t\geq
0}$ by comparison to the sequences $\\{\bar{Q}^{i}_{t}\\}_{t\geq 0}$ and
$\\{\bar{J}^{i}_{t}\\}_{t\geq 0}$.

The sequences $\\{\widehat{Q}^{i}_{t}\\}_{t\geq 0}$ and
$\\{\bar{Q}^{i}_{t}\\}_{t\geq 0}$ are related through the Q-factor update.
There are, however, two major differences. First, Algorithm 3 instructs player
$i$ to reset its counters at the end of the $k^{th}$ exploration phase (i.e.
after the update at time $t_{k+1}$, before the update at time $t_{k+1}+1$),
meaning the step sizes differ for the two iterate sequences
$\\{\bar{Q}^{i}_{t}\\}_{t\geq 0}$ and $\\{\widehat{Q}^{i}_{t}\\}_{t\geq 0}$
even when the state-action-cost trajectories observed by player $i$ are
identical. Second, Algorithm 3 instructs player $i$ to reset its Q-factors at
the end of the $k^{th}$ exploration phase, while Algorithm 1 does not involve
any resetting.

Consequently, one sees that the process $\\{\widehat{Q}^{i}_{t}\\}_{t\geq 0}$
depends on the initial condition
$\widehat{Q}^{i}_{0}=0\in\mathbb{R}^{\mathbb{Y}\times\mathbb{U}}$, the state-
action trajectory, and the resetting times $\\{t_{k}\\}_{k\geq 0}$. In
contrast, the process $\\{\bar{Q}^{i}_{t}\\}_{t\geq 0}$ depends only on the
initial value $\bar{Q}^{i}_{0}=0\in\mathbb{R}^{\mathbb{Y}\times\mathbb{U}}$
and the state-action trajectory. Analogous remarks hold relating
$\\{\widehat{J}^{i}_{t}\\}_{t\geq 0}$ and $\\{\bar{J}^{i}_{t}\\}_{t\geq 0}$.

Recall that the $k^{th}$ exploration phase begins with the stage game at time
$t_{k}$ and ends before the stage game at time $t_{k+1}=t_{k}+T_{k}$. During
the $k^{th}$ exploration phase, the sequences
$\\{\widehat{Q}^{i}_{t}\\}_{t=t_{k}}^{t_{k}+T_{k}}$ and
$\\{\widehat{J}^{i}_{t}\\}_{t=t_{k}}^{t_{k}+T_{k}}$ depend only on the the
state-action trajectory
$\textbf{x}_{t_{k}},\textbf{u}_{t_{k}},\cdots,\textbf{u}_{t_{k}+T_{k}-1},\textbf{x}_{t_{k}+T_{k}}$.
This leads to the following useful observation: for any
$(\textbf{s}_{0},\textbf{s}_{1},\cdots,\textbf{s}_{T_{k}})\in\mathbb{X}^{T_{k}+1}$
and $(\textbf{a}_{0},\cdots,\textbf{a}_{T_{k}-1})\in\textbf{U}^{T_{k}}$, we
have that

$\displaystyle{\rm
Pr}\left(\\{\textbf{x}_{t_{k}+T_{k}}=\textbf{s}_{T_{k}}\\}\bigcap_{j=0}^{T_{k}-1}\\{\textbf{x}_{t_{k}+j}=\textbf{s}_{j},\textbf{u}_{t_{k}+j}=\textbf{a}_{j}\\}\middle|\textbf{x}_{t_{k}}=\textbf{s},\bm{\pi}_{k}=\bm{\pi}\right)={\rm
Pr}^{\bm{\pi}}_{\textbf{s}}\left(\\{\textbf{x}_{T_{k}}=\textbf{s}_{T_{k}}\\}\bigcap_{j=0}^{T_{k}-1}\\{\textbf{x}_{j}=\textbf{s}_{j},\textbf{u}_{j}=\textbf{a}_{j}\\}\right).$

In words, once players following Algorithm 3 select a policy $\bm{\pi}$ for
the $k^{th}$ exploration phase, then the conditional probabilities of the
trajectories restricted to time indices in that exploration phase can be
described by ${\rm Pr}^{\bm{\pi}}$, with the indices of the events suitably
shifted to start at time 0. This leads to a series of useful lemmas, which we
include below for completeness.

###### Lemma 21

For any $i\in\mathcal{N},\bm{\pi}\in\bm{\Pi},k\geq 0$ and global state
$\textbf{s}\in\textbf{X}$, we have

${\rm
Pr}\left(\widehat{Q}^{i}_{t_{k+1}}\in\cdot\middle|\bm{\pi}_{k}=\bm{\pi},\textbf{x}_{t_{k}}=\textbf{s}\right)={\rm
Pr}^{\bm{\pi}}_{\textbf{s}}\left(\bar{Q}^{i}_{T_{k}}\in\cdot\right)$

and

${\rm
Pr}\left(\widehat{J}^{i}_{t_{k+1}}\in\cdot\middle|\bm{\pi}_{k}=\bm{\pi},\textbf{x}_{t_{k}}=\textbf{s}\right)={\rm
Pr}^{\bm{\pi}}_{\textbf{s}}\left(\bar{J}^{i}_{T_{k}}\in\cdot\right).$

Combining Lemma 21 with Theorem 2, we get the following result.

###### Lemma 22

For any joint policy $\bm{\pi}\in\bm{\Pi}$, global state
$\textbf{s}\in\textbf{X}$, and player $i\in\mathcal{N}$, we have

  1. 1.

${\rm
Pr}^{\bm{\pi}}_{\textbf{s}}\left(\lim_{t\to\infty}\bar{Q}^{i}_{t}=W^{*i}_{\bm{\pi}}\right)=1$;

  2. 2.

${\rm
Pr}^{\bm{\pi}}_{\textbf{s}}\left(\lim_{t\to\infty}\bar{J}^{i}_{t}=V^{*i}_{\bm{\pi}}\right)=1$;

  3. 3.

For any $\xi>0$, there exists $T=T(i,\bm{\pi},\xi)\in\mathbb{N}$ such that

${\rm Pr}^{\bm{\pi}}_{\textbf{s}}\left(\sup_{t\geq
T}\left\|\bar{Q}^{i}_{t}-W^{*}_{\bm{\pi}}\right\|<\xi\right)\geq 1-\xi,\text{
and }{\rm Pr}^{\bm{\pi}}_{\textbf{s}}\left(\sup_{t\geq
T}\left\|\bar{J}^{i}_{t}-V^{*}_{\bm{\pi}}\right\|<\xi\right)\geq 1-\xi$

Finally, we combine Lemmas 21 and 22 to obtain the following useful result on
conditional probabilities.

###### Lemma 23

Let $k,\ell\in\mathbb{Z}_{\geq 0}$ and suppose $k\leq\ell$. Let
$\mathcal{F}_{k}$ denote the $\sigma$-algebra generated by the random
variables $\bm{\pi}_{k}$ and $\textbf{x}_{t_{k}}$. For any $\xi>0$, there
exists $T=T(\xi)\in\mathbb{N}$ such that if $T_{\ell}\geq T$, then ${\rm
Pr}$-almost surely, we have

${\rm
Pr}\left(\bigcap_{i\in\mathcal{N}}\left\\{\left\|\widehat{Q}^{i}_{t_{\ell+1}}-W^{*i}_{\bm{\pi}_{\ell}}\right\|_{\infty}<\xi\right\\}\cap\left\\{\left\|\widehat{J}^{i}_{t_{\ell+1}}-V^{*i}_{\bm{\pi}_{\ell}}\right\|_{\infty}<\xi\right\\}\middle|\mathcal{F}_{k}\right)\geq
1-\xi.$

###  -D Proof of Theorem 6

###### Proof:

We begin by introducing the quantity $\bar{d}_{G}$, which will serve as the
upper bound for the tolerance parameters $d^{i},i\in\mathcal{N}$. The quantity
$\bar{d}_{G}$ depends on both $\epsilon>0$ and the subset of policies
$\bm{\Pi}\subset\bm{\Gamma}_{S}$ as follows:
$\bar{d}_{G}:=\min\mathcal{O}_{G},$ where
$\mathcal{O}_{G}=S_{G}\setminus\\{0\\}$ and $S_{G}$ is given by

$S_{G}:=\left\\{\left|\epsilon-\left(J^{i}_{\bm{\pi}}(\textbf{s})-\min_{a^{i}\in\mathbb{U}}W^{*i}_{\bm{\pi}}(\textbf{s},a^{i})\right)\right|:i\in\mathcal{N},\bm{\pi}\in\bm{\Pi},\textbf{s}\in\textbf{X}\right\\}.$

That is, $S_{G}$ is the collection of $\epsilon$-optimality gaps of the
various joint policies in $\bm{\Pi}$, and $\bar{d}_{G}$ is the minimum non-
zero separation between $\epsilon$ and the suboptimality gap of some player
$i$’s performance.

In Assumption 7, we assumed that each player’s $d^{i}$ parameter, which
represents tolerance for suboptimality in their policy evaluation step (Line
17 of Algorithm 3) and is included to account for noise in the learning
iterates, is positive and small: $d^{i}\in(0,\bar{d}_{G})$ for each
$i\in\mathcal{N}$.

We define
$\Xi:=\frac{1}{2}\min_{i\in\mathcal{N}}\\{d^{i},\bar{d}_{G}-d^{i}\\}$. The
quantity $\Xi$ will serve as a desirable upper bound on learning error: if
players jointly follow a joint policy $\bm{\pi}\in\bm{\Pi}$ and engage in
Q-learning and value function estimation, then convergence to within $\Xi$ of
the limiting values ensures that each player correctly assesses whether it is
$(\mathcal{V}^{*},\mathcal{W}^{*})$-subjectively $\epsilon$-best-responding by
using the comparison of Line 17 of Algorithm 3. We formalize this below.

For $k\geq 0$, we define ${\rm Event}(\Xi,k)$ as

${\rm
Event}(\Xi,k):=\left\\{\max\left\\{\|\widehat{J}^{i}_{t_{k+1}}-V^{*i}_{\bm{\pi}_{k}}\|_{\infty},\|\widehat{Q}^{i}_{t_{k+1}}-W^{*i}_{\bm{\pi}_{k}}\|_{\infty}:i\in\mathcal{N}\right\\}<\Xi\right\\}.$

Given ${\rm Event}(\Xi,k)$, we have that for any player $i\in\mathcal{N}$,

$\pi^{i}_{k}\in{\rm
Subj\text{-}BR}^{i}_{\epsilon}(\bm{\pi}^{-i}_{k},\mathcal{V}^{*},\mathcal{W}^{*})\iff\widehat{J}^{i}_{t_{k+1}}(\textbf{x})\leq\min_{a^{i}\in\mathbb{U}}\widehat{Q}^{i}_{t_{k+1}}(\textbf{x},a^{i})+\epsilon+d^{i}\quad\forall\textbf{x}\in\textbf{X}.$

For convenience, we also define the following intersection of events. For any
$k,\ell\in\mathbb{Z}_{\geq 0}$, let

$E_{k:k+\ell}:=\cap_{j=0}^{\ell}{\rm Event}(\Xi,k+j).$

That is, $E_{k:k+\ell}$ is the event where all agents obtain $\Xi$-accurate
learning estimates in each of the exploration phases $k,k+1,\cdots,k+\ell$.
For $k\in\mathbb{Z}_{\geq 0}$, we use
$G_{k}:=\\{\bm{\pi}_{k}\in\bm{\Pi}\cap{\rm
Subj}_{\epsilon}(\mathcal{V}^{*},\mathcal{W}^{*})\\}$ denote the event that
the policy $\bm{\pi}_{k}$ is an $\epsilon$-equilibrium.

From the preceding discussion on $\bar{d}_{G},\\{d^{i}\\}_{i\in\mathcal{N}}$
and the choice of $\Xi$, we have that for any $\ell\geq 0$,

${\rm Pr}\left(G_{k+\ell}\middle|G_{k}\cap E_{k:k+\ell}\right)=1.$ (10)

Recall the quantity $L:=\max\\{L_{\bm{\pi}_{0}}:\bm{\pi}_{0}\in\bm{\Pi}\\}$,
where for each $\bm{\pi}_{0}\in\bm{\Pi}$, the number $L_{\bm{\pi}_{0}}<\infty$
is defined as the shortest $(\mathcal{V}^{*},\mathcal{W}^{*})$-subjective
$\epsilon$-satisficing path within $\bm{\Pi}$ starting at $\bm{\pi}_{0}$ and
ending in $\bm{\Pi}\cap{\rm
Subj}_{\epsilon}(\mathcal{V}^{*},\mathcal{W}^{*})$. From our assumptions on
$\bm{\Pi}$, such a path exists for every $\bm{\pi}_{0}\in\bm{\Pi}$. If
$L_{\bm{\pi}_{0}}<L$ for a particular initial policy $\bm{\pi}_{0}$, we may
extend this path to have length $L$ by repeating the final term. Thus, for
every $\bm{\pi}\in\bm{\Pi}$, we obtain the inequality

${\rm Pr}\left(G_{k+L}\middle|\\{\bm{\pi}_{k}=\bm{\pi}\\}\cap
E_{k:k+L}\right)\geq p_{\min}>0,$ (11)

where
$p_{\min}:=\prod_{i\in\mathcal{N}}\left(\frac{e^{i}}{|\Pi^{i}|}\right)^{L}>0$.
The bound $p_{\min}$ is obtained through the following loose lower bounding
argument: beginning at $\bm{\pi}_{k}=\bm{\pi}$, the joint policy process
$\bm{\pi}_{k},\dots,\bm{\pi}_{k+L}$ follows the
$(\mathcal{V}^{*},\mathcal{W}^{*})$-subjective $\epsilon$-satificing path of
length $L$ described above with probability no less than the event that—at
each step—the “correct” unsatisfied player updates to the “correct” policy
uniformly at random, which occurs with probability no less than the
probability given by the ratio in the product.

Fix $u^{*}\in(0,1)$ such that

$\frac{u^{*}p_{\min}}{1-u^{*}+u^{*}p_{\min}}>1-\frac{\xi}{2}.$

Combining Lemma 23 with a union bound, we have that there exists
$\tilde{T}\in\mathbb{N}$ such that if $T_{l}>\tilde{T}$ for all $l\geq 0$,
then ${\rm Pr}\left(E_{k:k+L}|\bm{\pi}_{k}=\bm{\pi}\right)\geq u^{*}$ for all
$k\geq 0$ and any $\bm{\pi}\in\bm{\Pi}$. Thus, we have ${\rm
Pr}\left(E_{k:k+L}\middle|G_{k}\right)\geq u^{*}$ and ${\rm
Pr}\left(E_{k:k+L}\middle|G^{c}_{k}\right)\geq u^{*}$ for each $k\geq 0$.

We now lower bound ${\rm Pr}(G_{k+L})$ by conditioning on $G_{k}$ and
$G_{k}^{c}$ as follows.

$\displaystyle{\rm Pr}\left(G_{k+L}\right)={\rm Pr}(G_{k+L}|G_{k}){\rm
Pr}(G_{k})+{\rm Pr}(G_{k+L}|G_{k}^{c})(1-{\rm Pr}(G_{k})).$

We then lower bound each of the terms above by conditioning on $E_{k:k+L}$ and
invoking inequalities (10) and (11):

${\rm Pr}(G_{k+L})\geq 1\cdot{\rm Pr}(E_{k:k+L}|G_{k})\cdot{\rm
Pr}(G_{k})+p_{\min}{\rm Pr}(E_{k:k+L}|G_{k}^{c})(1-{\rm Pr}(G_{k})).$

Assume now that $T_{l}>\tilde{T}$ for all $l\geq 0$. We have

${\rm Pr}(G_{k+L})\geq u^{*}{\rm Pr}(G_{k})+u^{*}p_{\min}(1-{\rm
Pr}(G_{k})),\quad\forall k\geq 0.$

For each $k\in\\{0,1,\dots,L-1\\}$, define $y^{(k)}_{0}:={\rm Pr}(G_{k})$, and
for $m\geq 0$ define
$y^{(k)}_{m+1}:=u^{*}y^{(k)}_{m}+u^{*}p_{\min}\left(1-y^{(k)}_{m}\right).$ By
induction, one can show that

${\rm Pr}(G_{k+mL})\geq y^{(k)}_{m},\quad\forall m\geq 0.$ (12)

Observe that $y^{(k)}_{m+1}$ can be written as

$y^{(k)}_{m+1}=\left(u^{*}-u^{*}p_{\min}\right)^{m+1}y^{(k)}_{0}+u^{*}p_{\min}\sum_{j=0}^{m}\left(u^{*}-u^{*}p_{\min}\right)^{j}.$

Since $0<u^{*}-u^{*}p_{\min}<1$, we have that
$\lim_{m\to\infty}y^{(k)}_{m}=\frac{u^{*}p_{\min}}{1-u^{*}+u^{*}p_{\min}}>1-\frac{\xi}{2}$.
Then, by (12), we have that ${\rm Pr}(G_{k+mL})\geq 1-\xi/2$ holds for all
sufficiently large $m$.

In the case of global state observability, softness of policies in $\bm{\Pi}$
means that $G_{k}$ is equivalent to
$\bm{\pi}_{k}\in\bm{\Pi}\cap\bm{\Gamma}^{\epsilon{\rm\text{-}eq}}$, which
completes the proof.

∎

###  -E Proof of Theorem 7

This proof parallels that of Theorem 6, with only minor modifications: we
assume that for each $i\in\mathcal{N}$ we have $d^{i}\in(0,\bar{d}_{{\rm
MF}})$, where $\bar{d}_{{\rm MF}}$ is defined as $\bar{d}_{{\rm
MF}}:=\min\mathcal{O}_{{\rm MF}}$, where $\mathcal{O}_{{\rm MF}}:=S_{{\rm
MF}}\setminus\\{0\\}$, and $S_{{\rm MF}}$ is given by

$S_{{\rm
MF}}:=\left\\{\left|\epsilon-\left(V^{*i}_{\bm{\pi}}(y)-\min_{a^{i}\in\mathbb{U}}W^{*i}_{\bm{\pi}}(y,a^{i})\right)\right|:i\in\mathcal{N},\bm{\pi}\in\bm{\Pi},y\in\mathbb{Y}\right\\}.$

Then, the upper bound on learning accuracy $\Xi$ is here defined to be
$\Xi:=\frac{1}{2}\min_{i\in\mathcal{N}}\\{d^{i},\bar{d}_{{\rm MF}}-d^{i}\\}$.
The rest of the proof goes through essentially unchanged.

###  -F Proof of Theorem 8

As above, this proof also parallels that of Theorem 6, with only minor
modifications. We define $\bar{d}_{{\rm comp}}$ as $\bar{d}_{{\rm
comp}}:=\min\mathcal{O}_{{\rm comp}}$, where $\mathcal{O}_{{\rm
comp}}:=S_{{\rm comp}}\setminus\\{0\\}$ and

$S_{{\rm
comp}}:=\left\\{\left|\epsilon-\left(V^{*i}_{\bm{\pi}}(y)-\min_{a^{i}\in\mathbb{U}}W^{*i}_{\bm{\pi}}(y,a^{i})\right)\right|:i\in\mathcal{N},\bm{\pi}\in\bm{\Pi},y\in\mathbb{Y}\right\\}.$

The upper bound on learning accuracy $\Xi$ is here defined to be
$\Xi:=\frac{1}{2}\min_{i\in\mathcal{N}}\\{d^{i},\bar{d}_{{\rm
comp}}-d^{i}\\}$. The rest of the proof goes through essentially unchanged.

## References

  * [1] M. Huang, R. P. Malhamé, and P. E. Caines, “Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle,” Communications in Information & Systems, vol. 6, no. 3, pp. 221–252, 2006.
  * [2] M. Huang, P. E. Caines, and R. P. Malhamé, “Large-population cost-coupled LQG problems with nonuniform agents: individual-mass behavior and decentralized $\epsilon$-Nash equilibria,” IEEE Transactions on Automatic Control, vol. 52, no. 9, pp. 1560–1571, 2007.
  * [3] J.-M. Lasry and P.-L. Lions, “Mean field games,” Japanese journal of mathematics, vol. 2, no. 1, pp. 229–260, 2007.
  * [4] R. Carmona, F. Delarue, and D. Lacker, “Mean field games with common noise,” The Annals of Probability, vol. 44, no. 6, pp. 3740–3803, 2016.
  * [5] D. Lacker, “A general characterization of the mean field limit for stochastic differential games,” Probability Theory and Related Fields, vol. 165, no. 3, pp. 581–648, 2016.
  * [6] M. Fischer, “On the connection between symmetric $n$-player games and mean field games,” The Annals of Applied Probability, vol. 27, no. 2, pp. 757–810, 2017.
  * [7] N. Saldi, T. Basar, and M. Raginsky, “Markov–nash equilibria in mean-field games with discounted cost,” SIAM Journal on Control and Optimization, vol. 56, no. 6, pp. 4256–4287, 2018.
  * [8] S. Sanjari, N. Saldi, and S. Yüksel, “Optimality of independently randomized symmetric policies for exchangeable stochastic teams with infinitely many decision makers,” Mathematics of Operations Research, 2022\.
  * [9] O. Guéant, J.-M. Lasry, and P.-L. Lions, “Mean field games and applications,” in Paris-Princeton lectures on mathematical finance 2010, pp. 205–266, Springer, 2011.
  * [10] D. Bauso, H. Tembine, and T. Başar, “Robust mean field games with application to production of an exhaustible resource,” IFAC Proceedings Volumes, vol. 45, no. 13, pp. 454–459, 2012.
  * [11] D. Gomes, R. M. Velho, and M.-T. Wolfram, “Socio-economic applications of finite state mean field games,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 372, no. 2028, 2014.
  * [12] C. Yang, J. Li, P. Semasinghe, E. Hossain, S. M. Perlaza, and Z. Han, “Distributed interference and energy-aware power control for ultra-dense d2d networks: A mean field game,” IEEE Transactions on Wireless Communications, vol. 16, no. 2, pp. 1205–1217, 2016.
  * [13] G. Chevalier, J. Le Ny, and R. Malhamé, “A micro-macro traffic model based on mean-field games,” in 2015 American Control Conference (ACC), pp. 1983–1988, IEEE, 2015.
  * [14] M. L. Littman, “Markov games as a framework for multi-agent reinforcement learning,” in Machine Learning Proceedings 1994, pp. 157–163, Elsevier, 1994.
  * [15] K. Zhang, Z. Yang, and T. Başar, “Multi-agent reinforcement learning: A selective overview of theories and algorithms,” Handbook of Reinforcement Learning and Control, pp. 321–384, 2021.
  * [16] C. Berner, G. Brockman, B. Chan, V. Cheung, et al., “Dota 2 with large scale deep reinforcement learning,” arXiv preprint arXiv:1912.06680, 2019\.
  * [17] C. Daskalakis, D. J. Foster, and N. Golowich, “Independent policy gradient methods for competitive reinforcement learning,” arXiv preprint arXiv:2101.04233, 2021.
  * [18] M. O. Sayin, K. Zhang, D. S. Leslie, T. Basar, and A. Ozdaglar, “Decentralized Q-learning in zero-sum Markov games,” arXiv preprint arXiv:2106.02748, 2021.
  * [19] G. Arslan and S. Yüksel, “Decentralized Q-learning for stochastic teams and games,” IEEE Transactions on Automatic Control, vol. 62, no. 4, pp. 1545–1558, 2017.
  * [20] B. Yongacoglu, G. Arslan, and S. Yüksel, “Decentralized learning for optimality in stochastic dynamic teams and games with local control and global state information,” IEEE Transactions on Automatic Control, to appear.
  * [21] M. L. Littman and C. Szepesvári, “A generalized reinforcement-learning model: Convergence and applications,” in ICML, vol. 96, pp. 310–318, Citeseer, 1996.
  * [22] J. Hu and M. P. Wellman, “Nash q-learning for general-sum stochastic games,” Journal of Machine Learning Research, vol. 4, no. Nov, pp. 1039–1069, 2003\.
  * [23] M. L. Littman, “Friend-or-foe Q-learning in general-sum games,” in ICML, vol. 1, pp. 322–328, 2001.
  * [24] L. Wang, Z. Yang, and Z. Wang, “Breaking the curse of many agents: Provable mean embedding q-iteration for mean-field reinforcement learning,” in International Conference on Machine Learning, pp. 10092–10103, PMLR, 2020.
  * [25] L. Matignon, G. J. Laurent, and N. Le Fort-Piat, “Coordination of independent learners in cooperative Markov games.,” HAL preprint hal-00370889, 2009\.
  * [26] L. Matignon, G. J. Laurent, and N. Le Fort-Piat, “Independent reinforcement learners in cooperative markov games: a survey regarding coordination problems.,” Knowledge Engineering Review, vol. 27, no. 1, pp. 1–31, 2012\.
  * [27] A. Condon, “On algorithms for simple stochastic games.,” Advances in computational complexity theory, vol. 13, pp. 51–72, 1990.
  * [28] M. Tan, “Multi-agent reinforcement learning: Independent vs. cooperative agents,” in Proceedings of the Tenth International Conference on Machine Learning, pp. 330–337, 1993.
  * [29] S. Sen, M. Sekaran, and J. Hale, “Learning to coordinate without sharing information,” in Proceedings of the 12th National Conference on Artificial Intelligence, pp. 426–431, 1994.
  * [30] C. Claus and C. Boutilier, “The dynamics of reinforcement learning in cooperative multiagent systems,” in Proceedings of the Tenth Innovative Applications of Artificial Intelligence Conference, Madison, Wisconsin, pp. 746–752, 1998.
  * [31] P. Peng, Y. Wen, Y. Yang, Q. Yuan, Z. Tang, H. Long, and J. Wang, “Multiagent bidirectionally-coordinated nets: Emergence of human-level coordination in learning to play starcraft combat games,” arXiv preprint arXiv:1703.10069, 2017.
  * [32] R. Lowe, Y. I. Wu, A. Tamar, J. Harb, O. Pieter Abbeel, and I. Mordatch, “Multi-agent actor-critic for mixed cooperative-competitive environments,” Advances in neural information processing systems, vol. 30, 2017.
  * [33] J. Foerster, G. Farquhar, T. Afouras, N. Nardelli, and S. Whiteson, “Counterfactual multi-agent policy gradients,” in Proceedings of the AAAI conference on artificial intelligence, vol. 32, 2018.
  * [34] A. C. Kizilkale and P. E. Caines, “Mean field stochastic adaptive control,” IEEE Transactions on Automatic Control, vol. 58, no. 4, pp. 905–920, 2012\.
  * [35] H. Yin, P. G. Mehta, S. P. Meyn, and U. V. Shanbhag, “Learning in mean-field games,” IEEE Transactions on Automatic Control, vol. 59, no. 3, pp. 629–644, 2013.
  * [36] Y. Yang, R. Luo, M. Li, M. Zhou, W. Zhang, and J. Wang, “Mean field multi-agent reinforcement learning,” in International Conference on Machine Learning, pp. 5571–5580, PMLR, 2018.
  * [37] X. Guo, A. Hu, R. Xu, and J. Zhang, “Learning mean-field games,” Advances in Neural Information Processing Systems, vol. 32, 2019.
  * [38] J. Subramanian and A. Mahajan, “Reinforcement learning in stationary mean-field games,” in Proceedings of the 18th International Conference on Autonomous Agents and MultiAgent Systems, pp. 251–259, 2019.
  * [39] R. Elie, J. Perolat, M. Laurière, M. Geist, and O. Pietquin, “On the convergence of model free learning in mean field games,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 34, pp. 7143–7150, 2020\.
  * [40] Q. Xie, Z. Yang, Z. Wang, and A. Minca, “Learning while playing in mean-field games: Convergence and optimality,” in International Conference on Machine Learning, pp. 11436–11447, PMLR, 2021.
  * [41] M. A. uz Zaman, K. Zhang, E. Miehling, and T. Bașar, “Reinforcement learning in non-stationary discrete-time linear-quadratic mean-field games,” in 2020 59th IEEE Conference on Decision and Control (CDC), pp. 2278–2284, IEEE, 2020.
  * [42] B. Anahtarci, C. D. Kariksiz, and N. Saldi, “Learning in discrete-time average-cost mean-field games,” in 2021 60th IEEE Conference on Decision and Control (CDC), pp. 3048–3053, IEEE, 2021.
  * [43] O. Hernández-Lerma, Adaptive Markov control processes, vol. 79.  Springer Science & Business Media, 2012.
  * [44] C. Watkins, Learning from Delayed Rewards.  PhD thesis, Cambridge University, 1989.
  * [45] E. Wei and S. Luke, “Lenient learning in independent-learner stochastic cooperative games,” The Journal of Machine Learning Research, vol. 17, no. 1, pp. 2914–2955, 2016.
  * [46] B. Yongacoglu, G. Arslan, and S. Yüksel, “Satisficing paths and independent multi-agent reinforcement learning in stochastic games,” arXiv preprint arXiv:2110.04638, 2022.
  * [47] A. M. Fink, “Equilibrium in a stochastic $n$-person game,” Journal of Science of the Hiroshima University, Series AI (Mathematics), vol. 28, no. 1, pp. 89–93, 1964.
  * [48] P. Hernandez-Leal, M. Kaisers, T. Baarslag, and E. M. de Cote, “A survey of learning in multiagent environments: Dealing with non-stationarity,” arXiv preprint arXiv:1707.09183, 2017.
  * [49] A. D. Kara and S. Yuksel, “Convergence of finite memory q-learning for pomdps and near optimality of learned policies under filter stability,” arXiv preprint arXiv:2103.12158, 2021.
  * [50] J. N. Tsitsiklis, “Asynchronous stochastic approximation and Q-learning,” Machine Learning, vol. 16, no. 3, pp. 185–202, 1994.
  * [51] G. Arslan and S. Yüksel, “Subjective equilibria under beliefs of exogenous uncertainty in stochastic dynamic games,” arXiv preprint arXiv:2005.01640, 2020.
  * [52] S. Adlakha, R. Johari, and G. Y. Weintraub, “Equilibria of dynamic games with many players: Existence, approximation, and market structure,” Journal of Economic Theory, vol. 156, pp. 269–316, 2015.
  * [53] G. Y. Weintraub, L. Benkard, and B. Van Roy, “Oblivious equilibrium: A mean field approximation for large-scale dynamic games,” Advances in neural information processing systems, vol. 18, 2005.
  * [54] G. Y. Weintraub, C. L. Benkard, and B. Van Roy, “Markov perfect industry dynamics with many firms,” Econometrica, vol. 76, no. 6, pp. 1375–1411, 2008.
  * [55] S. Chien and A. Sinclair, “Convergence to approximate nash equilibria in congestion games,” Games and Economic Behavior, vol. 71, no. 2, pp. 315–327, 2011.
  * [56] G. C. Chasparis, A. Arapostathis, and J. S. Shamma, “Aspiration learning in coordination games,” SIAM J. Control and Optimization, vol. 51, no. 1, pp. 465–490, 2013.
  * [57] O. Candogan, A. Ozdaglar, and P. A. Parrilo, “Near-potential games: Geometry and dynamics,” ACM Transactions on Economics and Computation (TEAC), vol. 1, no. 2, pp. 1–32, 2013.
  * [58] P. Mertikopoulos, C. Papadimitriou, and G. Piliouras, “Cycles in adversarial regularized learning,” in Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 2703–2717, SIAM, 2018.
  * [59] E. Mazumdar, L. J. Ratliff, M. I. Jordan, and S. S. Sastry, “Policy-gradient algorithms have no guarantees of convergence in linear quadratic games,” arXiv preprint arXiv:1907.03712, 2019.
  * [60] E. Mazumdar, L. J. Ratliff, and S. S. Sastry, “On gradient-based learning in continuous games,” SIAM Journal on Mathematics of Data Science, vol. 2, no. 1, pp. 103–131, 2020.