#  Anomalous diffusion of the heavy quarks through the fractional Langevin
equation

Jai Prakash<EMAIL_ADDRESS>School of Physical Sciences, Indian
Institute of Technology Goa, Ponda-403401, Goa, India

###### Abstract

The dynamics of heavy quarks within the hot QCD medium have been revisited,
considering the influence of anomalous diffusion. This study has been
conducted using the framework of the fractional Langevin equation involving
the Caputo fractional derivative. We introduce a numerical scheme for the
fractional Langevin equation and demonstrate that the mean square displacement
of the particle exhibits anomalous diffusion, deviating from a linear
relationship with time. Our analysis calculates various entities, such as mean
squared momentum, momentum spread, and the nuclear suppression factor,
$R_{AA}$. Notably, our findings indicate that superdiffusion strongly
suppresses the $R_{AA}$ compared to normal diffusion in the hot QCD medium.
The possible impacts on other parameters are also discussed.

Heavy quark, quark-gluon plasma, fractional Langevin equation, anomalous
diffusion, Caputo fractional derivative.

##  I Introduction

The ultra-relativistic heavy-ion collisions (HICs) at the Relativistic Heavy-
Ion Collider (RHIC) and the Large Hadron Collider (LHC) have predicted the
existence of Quark-Gluon Plasma (QGP); a state of matter where quarks and
gluons are free to move beyond the nucleonic volume Adams  _et al._ (2006,
2005); Adcox  _et al._ (2005); Aamodt  _et al._ (2010); Arsene  _et al._
(2005). The QGP is short-lived, with an expected lifespan of a few fm/c
(approximately 4-5 fm/c at RHIC and 10-12 fm/c at LHC van Hees and Rapp
(2005); Rapp and van Hees (2010)). Conversely, investigating the
characteristics of the QGP by studying the dynamics of heavy quarks (charm and
beauty) is a subject of significant interest. In this context, The heavy
quarks (HQs) serve as prominent probes of QGP Song  _et al._ (2015); Andronic
_et al._ (2016); Dong and Greco (2019); Cao  _et al._ (2019); Beraudo  _et
al._ (2018); Prino and Rapp (2016); Aarts  _et al._ (2017); Uphoff  _et al._
(2011); Golam Mustafa  _et al._ (1998); Plumari  _et al._ (2018); Gossiaux and
Aichelin (2008); Prakash  _et al._ (2021, 2023); Prakash and Jamal (2023a);
Jamal  _et al._ (2023); Singh  _et al._ (2023); Kurian  _et al._ (2020); Cao
_et al._ (2016); Mazumder  _et al._ (2011); Zhang  _et al._ (2023); Jamal  _et
al._ (2021); Jamal and Mohanty (2021a, b); Sun  _et al._ (2023); Plumari  _et
al._ (2020); Prakash and Jamal (2023b); Debnath  _et al._ (2023); Du and Qian
(2023); Shaikh  _et al._ (2021); Kumar  _et al._ (2022); Das  _et al._ (2022).
Due to their large masses ($M_{c}/M_{b}$ $\sim$ 1.3/4.5 GeV), HQs are
generated at the early stages of HICs. Moreover, the thermalization of the HQs
is delayed compared to the light partons of the bulk medium by a factor on the
order of $\sim$ $M/T$. This delay renders the thermalization time of the HQs
comparable to the lifetime of the QGP fireball. Since the HQs are not expected
to achieve complete thermalization, they retain the memory of their
interaction history. Hence, they serve as a novel probe for observing the
complete evolution of the QGP medium and act as non-equilibrium entities
within the equilibrated QGP. The standard Langevin equation (LE) is used to
study Brownian motion under the assumption of normal diffusion, which means a
linear increase in the mean-square displacement of the particle with time for
a sufficiently large time Kubo  _et al._ (2012). By using the Langevin
equation Young  _et al._ (2012); Lang  _et al._ (2016); Cao and Bass (2011);
van Hees  _et al._ (2008); Das  _et al._ (2014) in the context of the HQs
momentum and position evolution in the hot QCD medium and numerous studies
available for the experimental observables related to the HQs, such as the
nuclear suppression factor, $R_{AA}$ and elliptic flow, $v_{2}$ He  _et al._
(2023, 2012); Das  _et al._ (2015); Zhang  _et al._ (2023). Several processes
exhibit anomalous diffusion, a phenomenon where the mean squared displacement
of the particle does not vary linearly with time as predicted by Einstein’s
law. Specifically, anomalous diffusion is characterized by a non-standard time
dependence of the mean squared displacement ($\langle x(t)^{2}\rangle$)
written as follows Balescu (1995); Bouchaud and Georges (1990); Metzler and
Klafter (2000); Caspi  _et al._ (2000); Weber  _et al._ (2010); Vitali  _et
al._ (2018); Oliveira  _et al._ (2019); Sokolov and Klafter (2005); Lutz
(2001); Lim and Teo (2009); Wang and Chen (2023),

$\displaystyle\langle x(t)^{2}\rangle\propto t^{\mu},\ {\mu}\neq 1,$ (1)

where $t$ is the evolution time of the particle. The process describes normal
diffusion, corresponding to the case where ${\mu}=1$. Subdiffusion occurs when
${\mu}<1$, whereas superdiffusion occurs when ${\mu}>1$ Metzler and Klafter
(2000). Despite the crucial role played by the LE across several fields, it
fails to accurately describe certain behaviours, such as anomalous diffusion
(i.e., superdiffusion and subdiffusion). Therefore, fractional Langevin
equations (FLE) have been proposed to study anomalous diffusion Li  _et al._
(2012); Kobelev and Romanov (2000); Coffey  _et al._ (2004). Mainardi et al.,
Gorenflo and Mainardi (1997); Mainardi and Pironi (1996); Mainardi (1997)
introduced an FLE in their groundbreaking research. The Langevin method has
primarily been applied to study anomalous diffusion using Caputo fractional
derivative Guo  _et al._ (2013); Li and Zeng (2013); Miller and Ross (1993).
On the other hand, there are several distinct formulations for fractional
derivatives, such as Riemann-Liouville derivative Carpinteri and Mainardi
(2014), Riesz derivative Agrawal (2007), Feller derivative Feller (1971), and
others. In the subsequent, we will only use the Caputo fractional derivative
in our analysis to study the anomalous diffusion of the particle in the
medium.

Another valuable approach to studying anomalous diffusion involves the
investigation of the fractional diffusion equation, fractional Fokker-Planck
equation Metzler  _et al._ (1999), and generalized Chapman-Kolmogorov equation
Metzler (2000). Anomalous diffusion has been used in many fields, including
molecular chemistry Zumofen and Klafter (1994), biology Lomholt  _et al._
(2005), and anomalous diffusion, polymer transport theory Shlesinger  _et al._
(1993). Also, when accounting for interactions between the Brownian particle
and the constituent particles of the medium, the fluctuations are influenced
by their prior states, a phenomenon termed as memory. The current movement is
affected by previous movements via a memory kernel in the generalized Langevin
equation. Such memory effects can result in anomalous diffusion for a
particular memory time Wang and Tokuyama (1999); Wang (1992a); Porrà  _et al._
(1996); Wang (1992b); Fa (2006); Siegle  _et al._ (2010); Chen  _et al._
(2023, 2023) and also plays a role in the QGP Hammelmann  _et al._ (2019);
Kapusta  _et al._ (2012); Ruggieri  _et al._ (2022); Pooja  _et al._ (2023).

Recently, the transverse momentum broadening of a fast parton via
superdiffusion in the QCD matter has been studied in Ref. Caucal and Mehtar-
Tani (2022). With this motivation for the applications of superdiffusion on
the HQs dynamics in hot QCD matter, we use FLE, a generalized form of the LE.
Unlike the LE, the FLE replaces integer-order derivatives with fractional-
order derivatives, specifically of Caputo type Caputo (1967). This may be the
first attempt where we present an anomalous diffusion for the HQ dynamics in
the QGP medium. In this article, we specifically study the effect of
superdiffusion on the HQ dynamics in the QGP medium. The subdiffusion is not
considered within the scope of the current analysis.

We anticipate a key outcome in our study, specifically that in the presence of
superdiffusion, there will be an increase in the energy loss experienced by
the HQs within the QGP medium. This conclusion is confirmed by studying
various key parameters, including the mean squared displacement, mean squared
momentum, momentum spread, $dN/dp_{T}$, and the $R_{AA}$ of the HQs. Notably,
the latter holds particular significance in the phenomenology of the HQs
within the QGP. The strong suppression observed in $R_{AA}$ is anticipated to
contribute significantly to developing a substantial elliptic flow ($v_{2}$)
for the HQs.

The article is structured as follows: Section II introduces the formalism; we
analytically solve the FLE of non-relativistic heavy particles using the
Laplace technique, along with presenting a numerical scheme for solving the
relativistic FLE. Section III is dedicated to the presentation of our results.
Lastly, Section IV provides a comprehensive summary of our conclusions.

##  II The fraction Langevin equation for Brownian particles

For illustrative purposes, we discuss the scenario of a heavy particle with
mass $M$ undergoing one-dimensional (1-D) motion within the nonrelativistic
limit. We use the FLE to describe the dynamics of this particle. In this
framework, the evolution of the particle’s position and momentum is
characterized by fractional derivatives of orders $\beta$ and $\alpha$ Kobelev
and Romanov (2000); Li  _et al._ (2012),

${}^{C}D^{\beta}_{0+}x(t)=\frac{p(t)}{M},$ (2)
${}^{C}D^{\alpha}_{0+}p(t)=-\gamma p(t)+\xi(t),$ (3)

where the momentum of a particle at time $t$ is denoted as $p(t)$, its
position as $x(t)$. ${}^{C}D_{0+}^{\alpha}$ denotes Caputo fractional
derivative, $\alpha$ and $\beta$ are the fractional parameter, with
$n-1<\alpha\leq n$, and ${}^{C}D_{0+}^{\beta}$ with $n-1<\beta\leq n$,
$n\in\mathbb{N}$ ($n$ is a natural number). Brownian particle encounters two
distinct forces: the dissipative force, characterized by the drag coefficient,
denoted by $\gamma$, and the stochastic force, denoted as $\xi(t)$. The latter
governs the random noise, commonly referred to as white Gaussian noise. White
noise gives rise to a fluctuating field without memory, characterized by
instantaneous decay in correlations of the white noise, often referred to as a
$\delta$ correlation. The random force satisfies certain properties, such as:

$\displaystyle\langle\xi(t)\xi(t^{\prime})\rangle=2\mathcal{D}\delta(t-t^{\prime}),$
(4) $\displaystyle\langle\xi(t)\rangle=0,$ (5)

$\mathcal{D}$ is the diffusion coefficient of the heavy particle in a medium
of temperature, $T$. The drag coefficient ($\gamma$) is related to the
diffusion coefficient through the Fluctuation-Dissipation Theorem (FDT) as
follows:

$\displaystyle\gamma=\frac{\mathcal{D}}{MT}.$ (6)

The Caputo fractional derivative Caputo (1967) is defined as,

$\displaystyle^{C}D^{\nu}_{0+}u(t)=\frac{1}{\Gamma({n-\nu})}\int_{0}^{t}\frac{u^{(n)}(s)}{(t-s)^{1+\nu-n}}ds,$
(7)

where $u^{(n)}$ denotes the $n^{th}$ derivative of $u$. The fractional
derivative which corresponds to the superdiffusion, i.e., when $1<\nu\leq 2$,
is given by

$\displaystyle^{C}D^{\nu}_{0+}u(t)=\frac{1}{\Gamma({2-\nu})}\int_{0}^{t}\frac{u^{(2)}(s)}{(t-s)^{\nu-1}}ds,$
(8)

where $\Gamma(\cdot)$ denotes the gamma function. In the following subsection,
we solve analytically the FLE of non-relativistic Brownian particles.

###  II.1 Analytical solutions

Analytically, we solve the FLE governing the motion of non-relativistic
Brownian particles. To achieve this, we use the Laplace technique, providing
analytical expressions for $\langle p^{2}(t)\rangle$ and $\langle
x^{2}(t)\rangle$. Performing the Laplace transformation of Caputo derivative
for $n-1<\alpha\leq n$ and $n-1<\beta\leq n$ Podlubny ,

$\displaystyle\mathcal{L}\left[{}^{C}D^{\alpha}_{0+}p(t)\right](z)=z^{\alpha}\widehat{p}(z)-\sum_{k=0}^{n-1}z^{\alpha-k-1}\left[p^{k}(0)\right].$
(9)

In order to calculate the solution of $p(t)$ and $x(t)$ of the heavy particle
for a super-diffusive process, i.e., when $1<\alpha\leq 2$ and $1<\beta\leq
2$, one can take the Laplace transform Eqs. (2), (3) and using Eq. (9), it can
be obtained easily as

$\displaystyle
z^{\beta}\widehat{x}(z)-z^{\beta-1}x(t_{0})-z^{\beta-2}x^{\prime}(t_{0})=\frac{\widehat{p}(z)}{M},$
(10)
$\displaystyle\widehat{p}(z)=\frac{z^{\alpha-1}}{z^{\alpha}+\gamma}p(t_{0})+\frac{z^{\alpha-2}}{z^{\alpha}+\gamma}p^{\prime}(t_{0})+\frac{\widehat{\xi}(z)}{z^{\alpha}+\gamma}.$
(11)

Now, substituting Eq. (11) into Eq. (10), for simplicity, we take $M=1$ Sandev
_et al._ (2012), we obtain

$\displaystyle\widehat{x}(z)=$
$\displaystyle\;\frac{x(t_{0})}{z}+\frac{x^{\prime}(t_{0})}{z^{2}}+\frac{z^{\alpha-\beta-1}}{z^{\alpha}+\gamma}p(t_{0})+\frac{z^{\alpha-\beta-2}}{z^{\alpha}+\gamma}p^{\prime}(t_{0})$
$\displaystyle+\frac{z^{-\beta}}{z^{\alpha}+\gamma}\widehat{\xi}(z).$ (12)

Taking the inverse Laplace transform of Eqs. (11) and (II.1), we obtain

$\displaystyle p(t)=$ $\displaystyle\;E_{\alpha,1}(-\gamma
t^{\alpha})p(t_{0})+tE_{\alpha,2}(-\gamma t^{\alpha})p^{\prime}({t_{0}})$
$\displaystyle+\int_{0}^{t}(t-s)^{\alpha-1}E_{\alpha,\alpha}(-\gamma(t-s)^{\alpha})\xi(s)ds,$
(13)

and

$\displaystyle x(t)=$
$\displaystyle\;x_{0}+tx^{\prime}_{0}+t^{\beta}E_{\alpha,\beta+1}(-\gamma
t^{\alpha})p(t_{0})$ $\displaystyle+t^{\beta+1}E_{\alpha,\beta+2}(-\gamma
t^{\alpha})p^{\prime}({t_{0}})$
$\displaystyle+\int_{0}^{t}(t-s)^{\beta+\alpha-1}E_{\alpha,\beta+\alpha}(-\gamma(t-s)^{\alpha})\xi(s)ds,$
(14)

Figure 1: $\langle p^{2}(t)\rangle$ versus time (left panel) and $\langle
x^{2}(t)\rangle$ versus time (right panel) for a 1-D, purely diffusive motion
for different values of $\alpha$ and $\beta$. Analytic solutions are
represented by solid lines, and numerical (Numer.) solutions are represented
by dashed lines.

where

$E_{\alpha,\beta}(z)=\sum^{\infty}_{j=0}\frac{z^{j}}{\Gamma(j\alpha+\beta)}$

is the two-parameter Mittag-Leffler function for $\alpha>0$, $\beta>0$, and
$z$ is a complex number Erdélyi  _et al._ (1953) and the Laplace transform of
the two-parameter Mittag-Leffler function is Hilfer (2000),

$\displaystyle\mathcal{L}\left[x^{\beta-1}E_{\alpha,\beta}(-\lambda
x^{\alpha})\right](z)=\frac{z^{\alpha-\beta}}{z^{\alpha}+\lambda}.$ (15)

Next we calculate $\langle p^{2}(t)\rangle$, and $\langle x^{2}(t)\rangle$
analytically for the Brownian particle.

####  II.1.1 Purely diffusive motion

We are calculating the $\langle p^{2}(t)\rangle$, and $\langle
x^{2}(t)\rangle$ for a purely diffusive (when $\gamma=0$ in the FLE), 1-D
motion. Then Eq. (3) simplified to

${}^{C}D^{\alpha}_{0+}p(t)=\xi(t),$ (16)

with initial momentum, $p(t_{0})=0$. Then to calculate $\langle
p^{2}(t)\rangle$ we follows Eq. (13) and Eq. (16), then we get

$\displaystyle\langle
p^{2}(t)\rangle=\frac{1}{\Gamma(\alpha)^{2}}\int_{0}^{t}\int_{0}^{t}(t-s_{1})^{\alpha-1}(t-s_{2})^{\alpha-1}$
$\displaystyle\langle\xi(s_{1})\xi(s_{2})\rangle$ $\displaystyle
ds_{1}ds_{2}.$ (17)

Solving further Eq. (II.1.1), using the correlation properties defined in Eq.
(4), we have,

$\displaystyle\langle p^{2}(t)\rangle$
$\displaystyle=\frac{2\mathcal{D}}{\Gamma(\alpha)^{2}}\int_{0}^{t}(t-s)^{2\alpha-2}ds,$
(18) $\displaystyle\langle p^{2}(t)\rangle$
$\displaystyle=\frac{2\mathcal{D}}{\Gamma(\alpha)^{2}}\frac{t^{2\alpha-1}}{2\alpha-1}.$
(19)

Similarly, one can calculate analytical solutions of $\langle x^{2}(t)\rangle$
for $\gamma=0$, and with initial conditions, $x(t_{0})=0$, using Eqs. (II.1),
(16) as

$\displaystyle\langle x^{2}(t)\rangle\approx$
$\displaystyle\;\frac{2\mathcal{D}}{\Gamma(\alpha+\beta)^{2}}\int_{0}^{t}(t-s)^{2\alpha+2\beta-2}ds$
$\displaystyle=\frac{2\mathcal{D}}{\Gamma(\alpha+\beta)^{2}}\frac{t^{2\alpha+2\beta-1}}{2\alpha+2\beta-1}.$
(20)

We have plotted Fig. 1 for the case of purely diffusive motion, we depict the
variation of $\langle p^{2}(t)\rangle$ over time (left panel) and $\langle
x^{2}(t)\rangle$ (right panel), each for different values of $\alpha$ and
$\beta$, while keeping $\mathcal{D}$ constant at $0.1$ GeV2/fm. For
simplicity, a non-relativistic case has been considered with $M=1$ GeV. In
Fig. 1 (left panel), results obtained from Eq. (19) show that $\langle
p^{2}(t)\rangle$ exhibits nearly linear evolution with time when
$\alpha\rightarrow 1$ (as shown by black solid lines), which is giving
$\langle p^{2}(t)\rangle=2\mathcal{D}t$. In this scenario, the Caputo
fractional derivative reduces to an ordinary first-order derivative; in this
limit of $\alpha$, the FLE is reduced back to the LE. On the other hand, when
the value of $\alpha>1$, the linear increase in $\langle p^{2}(t)\rangle$ with
time converts into nonlinear growth over time. It becomes evident that the
diffusion process evolves more gradually during the initial time, particularly
up to $1.5$ fm/c. In contrast, around $2$ fm/c, the pattern is reversed, now a
higher value of $\alpha$, giving a faster diffusion. Similarly, with the same
input parameter, Fig. 1 (right panel) results obtained from Eq. (20) depicts
the mean square displacement, $\langle x^{2}(t)\rangle$ varies with time. In
the limit, $\alpha,\beta\rightarrow 1$(as shown by black solid lines),
$\langle x^{2}(t)\rangle$ varies with $t^{3}$ (as shown by black solid lines),
showing the normal diffusion. However, with value higher values of $\alpha$
and $\beta$, $\langle x^{2}(t)\rangle$ varies with the power of time greater
than cubic power.

###  II.2  A numerical scheme for the fractional Langevin equation

In classical stochastic differential equations driven by Brownian motion, the
Itô or Stratonovich stochastic calculus is typically employed for solutions.
However, these methods are not applicable to the FLE driven by fractional
Brownian motion, as it is not a semimartingale (see Ref. Rogers (1997) for
detailed proof). Although the Monte Carlo method is a reliable solution for
stochastic differential equations, it is not appropriate in this case because
it relies on independent sequences. Still, the sequences in fractional noise
are dependent. Therefore, we have employed L2 numerical schemes, as detailed
in Oldham and Spanier (1974). These schemes are among the most effective for
discretizing the Caputo fractional derivative and provide a key contribution
to this article. We numerically solve the FLE in Eqs. (2), (3). To validate
our numerical computations, we compared the numerical results with the
analytical ones for the 1-D FLE, considering $\gamma$ = 0. This step is
crucial to ensure the accuracy of our numerical method, which will be employed
in solving the relativistic Langevin equation that will be discussed in the
following subsection. To solve the FLE, numerical methods are commonly
categorized as indirect or direct. Since time-fractional differential
equations can generally be reformulated into integro-differential equations,
solving such equations corresponds to indirect methods. On the other hand,
direct methods focus on approximating the time-fractional derivative itself.
This aspect constitutes the core of our work, where our goal is to discretize
the fractional derivative directly without transforming the associated
differential equation into its integral form. Our numerical algorithm for
solving the FLE is based on the three-step scheme using the central difference
method shown in Appendix VI.

####  II.2.1 The fractional Langevin equation of relativistic Brownian
particles: QGP

We extended the solution of FLE defined in Eqs. (2),(3), to describe the HQs
dynamics in the QGP medium in the relativistic limit as follows,

${}^{C}D^{\beta}_{0+}x(t)=\frac{p(t)}{E(t)},$ (21)
${}^{C}D^{\alpha}_{0+}p(t)=-\gamma p(t)+\xi(t),$ (22)

where, $E=\sqrt{p^{2}+M^{2}}$ represents the energy, and $p$ denotes the
momentum of the HQs. The $\gamma$ can be related to the diffusion coefficient
through an FDT as Walton and Rafelski (2000); Moore and Teaney (2005);
Mazumder  _et al._ (2014),

$\displaystyle\gamma=\frac{\mathcal{D}}{ET}.$ (23)

Since both Eq. (21) and Eq. (22) are non-linear differential integral
equations, preventing analytical solutions through methods such as Laplace
transformations, as done for the non-relativistic case for Eqs. (2),(3). The
solution of relativistic FLE defined in Eq. (22) and Eq. (21) is possible
solely through numerical methods. We utilize the numerical approach to compute
$\langle x^{2}(t)\rangle$, $\langle p^{2}(t)\rangle$, the momentum
distribution, and $R_{AA}$ of the HQs. For superdiffusion, the relativistic
FLE can be written in the discrete form using Eq. (30) ( given in appendix VI)
as follows,

$\displaystyle\begin{cases}x(t_{1})=x(t_{0})+\left[\displaystyle\frac{p(t_{0})}{E(t_{0})}\right]k\;&:\;n=1,\\\
x(t_{2})=2x(t_{1})-x(t_{0})+\left[\displaystyle\frac{p(t_{1})}{E(t_{1})}\right]k\;&:\;n=2,\\\
x(t_{n})=2x(t_{n-1})-x(t_{n-2})+\Bigg{[}\displaystyle\frac{p(t_{n-1})}{E(t_{n-1})}&\\\
-\displaystyle\sum_{j=1}^{n-2}b_{j}\left(x(t_{n-j})-2x(t_{n-j-1})+x(t_{n-j-2})\right)+\displaystyle\frac{n^{2-\alpha}-(n-1)^{2-\alpha}}{k}(x(t_{1})-x(t_{0}))\Bigg{]}k\;&:\;n\geq
3,\end{cases}$ (24)

and

$\displaystyle\begin{cases}p(t_{1})=p(t_{0})+\left[-\gamma
p(t_{0})+\sqrt{2\mathcal{D}/\Delta{t}}\eta\right]k\;&:\;n=1,\\\
p(t_{2})=2p(t_{1})-p(t_{0})+\left[-\gamma
p(t_{1})+\sqrt{2\mathcal{D}/\Delta{t}}\eta\right]k\;&:\;n=2,\\\
p(t_{n})=2p(t_{n-1})-p(t_{n-2})+\Bigg{[}-\gamma
p(t_{n-1})+\sqrt{2\mathcal{D}/\Delta{t}}\eta&\\\
-\displaystyle\sum_{j=1}^{n-2}b_{j}\left(p(t_{n-j})-2p(t_{n-j-1})+p(t_{n-j-2})\right)+\displaystyle\frac{n^{2-\alpha}-(n-1)^{2-\alpha}}{k}(p(t_{1})-p(t_{0}))\Bigg{]}k\;&:\;n\geq
3,\end{cases}$ (25)

where $k=\Delta t^{\alpha}\Gamma(3-\alpha)$, with $t_{1}=t_{0}+\Delta{t}$,
$t_{n}=t_{n-1}+\Delta t$. Here $t_{0}$ is the initial time, and $\Delta{t}$ is
the time step.

Figure 2: $\langle p^{2}(t)\rangle$ versus time, for the various values of the
$\alpha$ and $\mathcal{D}=0.1$ GeV2/fm at $T=250$ MeV.

##  III Results

###  III.1  A check for non-relativistic

To analyze the numerical scheme of FLE, we compute $\langle x^{2}(t)\rangle$
from Eq. (24) and $\langle p^{2}(t)\rangle$ from Eq. (25), subsequently
comparing the results with the analytical solutions derived from Eq. (20) and
Eq. (19), respectively. The FLE has been solved for the 1-D motion of the
heavy particle within the nonrelativistic limit, as explained in Sec. II. For
illustrative purposes, we have considered constant diffusion coefficient,
$\mathcal{D}$ = 0.1 GeV2/fm and mass of the heavy particle, $M$ = 1 GeV and
$\gamma=0$ (i.e., pure diffusion case) in FLE. In the left panel of Fig. 1,
the numerically computed $\langle p^{2}(t)\rangle$ for $\alpha=1.001$ (brown
dashed line), $\alpha=1.2$ (grey dashed line), $\alpha=1.4$ (red dashed line),
and $\alpha=1.6$ (orange dashed line) is presented. Notably, these numerical
results align with the analytic solutions from Eq. (19) and match the other
values of the $\alpha$. For an additional test of our numerical approach, we
have calculated $\langle x^{2}(t)\rangle$ of the heavy particle as shown in
Fig. 1 (right panel). In the same figure, the numerically calculated $\langle
x^{2}(t)\rangle$ for $\alpha=\beta=1.001$ (brown dashed line) and
$\alpha=\beta=1.2$ (grey dashed line), $\alpha=\beta=1.4$ (red dashed line),
and $\alpha=\beta=1.6$ (orange dashed line) is displayed. Again, the numerical
results match the analytic results from Eq. (20). One can notice that the
numerical simulation agrees with the analytical result, indicating that our
numerical scheme works correctly.

###  III.2 The evolution of $\langle p^{2}(t)\rangle$ and $\langle
x^{2}(t)\rangle$ of the HQs.

For our analysis, the definition we used for $\langle p^{2}(t)\rangle$ is
given by

$\displaystyle\langle p^{2}(t)\rangle=\langle
p_{x}^{2}(t)+p_{y}^{2}(t)\rangle.$ (26)

We have computed $\langle p^{2}(t)\rangle$ over time for four different values
of $\alpha$, namely $\alpha=1.001$ (black line), $\alpha=1.2$ (red line),
$\alpha=1.4$ (green line), and $\alpha=1.6$ (blue line) at $T=250$ MeV, $M$ =
1.3 GeV and $\mathcal{D}=0.1$ GeV2/fm correspond to the charm quark. The
corresponding results are shown in Fig. 2. In this scenario, the initial
momentum is set to be $p_{x}(t_{0})=p_{y}(t_{0})=0$.

From the left panel of Fig. 2, it is evident that as the magnitude of $\alpha$
increases, the behaviour of the process indicates superdiffusion for the charm
quark within the hot QCD medium. The impact of superdiffusion becomes more
noticeable at higher values of $\alpha$. It is noticeable that as
$\alpha\rightarrow 1$, the superdiffusion process converges back to normal
diffusion. Additionally, in the later stages, the mean squared momentum
$\langle p^{2}(t)\rangle$ tends to approach $3MT$ Moore and Teaney (2005), and
also the FLE defined in Eq. (25) simplifies to the standard LE described in
Ref. Moore and Teaney (2005); Das  _et al._ (2014); He  _et al._ (2012); Young
_et al._ (2012); Lang  _et al._ (2016). This transition emphasizes the
connection between superdiffusive behaviour governed by the FLE and the
conventional diffusion process characterized by the standard LE. In the right
panel of Fig. 2, we present a subset of the results depicted in the left
panel, focusing on the early-time evolution of $\langle p^{2}(t)\rangle$ for
four different values of $\alpha$. It is worth noting that in the initial
time, the diffusion process gradually evolves for higher values of $\alpha$.
However, for times beyond $3$ fm/c, the trend reverses, and larger $\alpha$
corresponds to a faster diffusion, as anticipated in Fig. 2 (left panel).

In Fig. 3, we have calculated the evolution of $\langle x^{2}(t)\rangle$ over
time for three values of $\alpha=\beta=1.001,1.2,1.4$, maintaining other
parameters consistent with those illustrated in Fig. 2. The initial conditions
are set to be $x(t_{0})=y(t_{0})=0$. It can be noticed that both
$\alpha=\beta=1.2,1.4$, a distinct shift towards superdiffusion is observed
(shown in Fig. 3) as we discussed initially in Eq. (1). The larger values of
$\alpha$ and $\beta$ contribute to this notable change in behaviour,
underlining the complex dynamics of the system. As $\alpha\rightarrow 1$ and
$\beta\rightarrow 1$, the behaviour reflects normal diffusion. In this limit,
at a later time, $\langle x^{2}(t)\rangle$ exhibits a proportional
relationship with $t$, as described in Moore and Teaney (2005); Svetitsky
(1988). This puzzling behaviour of $\langle p^{2}(t)\rangle$ and $\langle
x^{2}(t)\rangle$ due to the anomalous diffusion, which was not explained
before in the context of the HQ dynamics in a hot QCD medium.

Figure 3: $\langle x^{2}(t)\rangle$ versus time, for the various values of the
$\alpha$, $\beta$ and $\mathcal{D}=0.1$ GeV2/fm at $T=250$ MeV.

In the following sections, we have calculated $R_{AA}$ and the charm quark
momentum distribution, $dN/dp_{T}$, to see the effect of superdiffusion
processes. To determine the interaction of the HQ with the thermalized bath
consisting of massless quarks and gluons, we employ perturbative Quantum
Chromodynamics (pQCD) transport coefficients for elastic processes with the
well-established diffusion coefficients Svetitsky (1988). Here, the Debye
mass, $m_{D}=g_{T}T$ screens the infrared divergence associated with the
$t$-channel diagrams.

Figure 4:  The $R_{AA}$ as a function of $p_{T}$, with $t$ = 6 fm/c, and at
two different temperatures ($T$ = 250 MeV (left panel) and $T$ = 350 MeV
(right panel) with four different values of $\alpha$.

###  III.3 Nuclear modification factor

To analyze the impact of the superdiffusion on the experimental observable, we
have calculated the nuclear suppression factor, $R_{AA}(p_{T})$ for the HQs,
which is defined as follows Moore and Teaney (2005),

$\displaystyle R_{AA}(p_{T})=\frac{f_{\tau_{f}}(p_{T})}{f_{\tau_{0}}(p_{T})}.$
(27)

The momentum spectrum, $f_{\tau_{f}}(p)$, of charm quarks calculated for the
time evolution, $\tau_{f}$ = $6$ fm/c in our computational results and
$f_{\tau_{i}}(p)$ is initial momentum distribution of the charm quark. The
initial momentum spectra, $f_{\tau_{i}}(p)$, is taken according to the fixed
order + next-to-leading log (FONLL) calculations, which has been shown to be
capable of reproducing the spectra of D-mesons produced in proton-proton
collisions through fragmentation Cacciari  _et al._ (2005, 2012), the initial
momentum spectrum is written as,

$\frac{dN}{d^{2}p{{}_{T}}}=\frac{x_{0}}{(x_{1}+p_{T})^{x_{2}}},$ (28)

where the parameters are estimated as follow; $x_{0}=6.365480\times{10}^{8}$,
$x_{1}=9.0$ and $x_{2}=10.27890$. $R_{AA}(p_{T})$ $\neq$ 1 implies that charm
quarks undergo interactions with the medium. These interactions lead to
modifications in the spectrum of charm quarks.

In Fig. 4, the behaviour of $R_{AA}$ is depicted as a function of $p_{T}$,
which are calculated using the FLE as defined in Eq. (25) for different values
of $\alpha$. The calculations are performed at two distinct temperatures, $T$=
250 MeV (left panel) and $T$= 350 MeV (right panel). At $T$= 250 MeV, the
dominating influence is the drag force across the entire range of $p_{T}$. As
$p_{T}$ increases, the significance of energy loss becomes more pronounced,
which can be noticed in Fig. 4 (left panel). When $\alpha>$ 1, the normal
diffusion converts into superdiffusion, it is evident for $\alpha=1.2$ (red
line), $\alpha=1.4$ (blue line), and $\alpha=1.6$ (green line). With an
increasing value of $\alpha$, there is a notable decrease in magnitude
$R_{AA}$ (more suppression) at high $p_{T}$. For $\alpha\rightarrow 1$
(depicted by the black line), the behaviour of the $R_{AA}$ aligns with normal
diffusion, consistent with findings available in the literature for the same
input parameters Lang  _et al._ (2016); Das  _et al._ (2014, 2015). Fig. 4
(right panel), corresponds to $T=350$ MeV, the observed behaviour is
attributed to the diffusion-dominated propagation of the HQs within the hot
QCD medium. This dominance of diffusion mechanisms effectively leads to the
diffusion of low-momentum charm quarks to higher momentum states. The high
temperature, $T=350$ MeV, enhances the significance of diffusion processes,
resulting in a distinct pattern in the $R_{AA}$ as compared to $R_{AA}$ at
$T$= 250 MeV. For $\alpha>$ 1, a significant reduction in the magnitude of
$R_{AA}$ is observed, indicating more suppression at high $p_{T}$. Notably,
for the highest considered value of $\alpha=1.6$, a larger proportion of
particles tends to remain at low $p_{T}$, a consequence of the superdiffusion
process. This behaviour underscores the complex dynamics associated with
superdiffusion and its impact on the $R_{AA}$ of the HQs in the QGP medium.

###  III.4 Momentum spread of HQs

Figure 5:  The evolution of charm quark momentum distribution at $\tau_{f}$ =
6 fm/c as a function of $p_{T}$, for various $\alpha$ and at two different
temperatures ($T$ = 250 MeV (above panel) and $T$ = 350 MeV (below panel).
Assuming an initial delta distribution centered at $p(t_{0})$ = 10 GeV.

We show the evolution of charm momentum distribution, $dN/dp_{T}$, at static
temperatures, $T$ = 250 MeV (top panel) and $T$ = 350 MeV (bottom panel) in
Fig. 5. The $dN/dp_{T}$ evolution is performed for various values of $\alpha$
at a final evolution time of $\tau_{f}$ = 6 fm/c. To understand the impact of
superdiffusion on charm quarks within the QGP, we take initial conditions
where all charm quarks are concentrated within an extremely narrow $p_{T}$
bin, creating a delta-like distribution at $p(t_{0})$ = 10 GeV (magenta line).
It is observed that the interaction of the HQs with the QGP medium results in
the spreading of $dN/dp_{T}$. Subsequently, the evolution of this distribution
is analyzed using the FLE as defined in Eq. (25). We have observed the
evolution of $dN/dp_{T}$ at higher values of $\alpha$ for the case of
superdiffusion, where, $\alpha=1.2$ (red line), $\alpha=1.4$ (blue line), and
$\alpha=1.6$ (green line). As depicted in Fig. 5, when $\alpha>1$, the
distribution $dN/dp_{T}$ undergoes a notable spread and average momentum
shifts towards lower values of $p_{T}$ under the influence of diffusion and
drag coefficients, respectively. Specifically, for $\alpha=1.6$ (represented
by the blue line), the extent of spreading is more pronounced, and the average
momentum shifts towards the lower $p_{T}$ compared to other $\alpha$ values
such as 1.001, 1.2, and 1.4. However, for $\alpha$ = 1.001 (black line),
corresponding to the normal diffusion coefficient as explained in Ref. Cao and
Bass (2011); Das  _et al._ (2014). At the same time, the total area under the
curve remains constant for all values of $\alpha$ and both temperatures.

##  IV Conclusion and outlooks

In this paper, we have discussed anomalous diffusion through the FLE with the
Caputo fractional derivative, specifically focusing on superdiffusion in the
context of the HQs in the QGP medium. Notably, the mean squared displacement
of the particle exhibits a power-law dependence on time (as shown in Fig. 3).
Our analysis discussed the scenarios where the values of $\alpha$ and $\beta$
were taken as 1.001, 1.2, 1.4, and 1.6, showcasing the effects of
superdiffusion. We have demonstrated that as $\alpha$ approaches 1, the
superdiffusion reverts to normal diffusion, verified by numerical and
analytical calculations for $\langle x^{2}(t)\rangle$ and $\langle
p^{2}(t)\rangle$ in the non-relativistic 1-D case (see Fig. 1). Several key
quantities characterizing the dynamics of the HQs under superdiffusion have
been computed, including the $\langle p^{2}(t)\rangle$, $\langle
x^{2}(t)\rangle$. Extending our analysis, we incorporated physical observables
$R_{AA}$ of the HQs in the QGP medium. We then shifted our focus to the
momentum spread, $dN/dp_{T}$, utilizing an initial momentum distribution at
$p(t_{0})=10$ GeV. The FLE, with the HQs moving under dissipative and random
forces, was solved with transport coefficients serving as input parameters.
Our findings indicate that superdiffusion results in more suppression in
$R_{AA}$. To select specific values for $\alpha$ and $\beta$ within the ranges
of $1<\alpha\leq 2$ and $1<\beta\leq 2$, it is essential to be able to
simultaneously describe the experimental observables, such as $R_{AA}$ and
$v_{2}$, for the entire measured range of $p_{T}$.

In the future, we plan to explore the impact of superdiffusion on the HQ
dynamics in the QGP medium, especially considering the time correlation of
thermal noise. This study devotes the groundwork for more realistic
conditions, which should incorporate an exact initial geometry and an
expanding medium in the near future. It is interesting to study the $v_{2}$ of
the HQ in the context of superdiffusion within an expanding QGP medium.
Superdiffusion might impact various observables, such as HQs directed flow,
particle correlations, etc. Given the simplifying assumptions made in our
current study, it is challenging to anticipate the specific modifications. A
more comprehensive and quantitative analysis will be performed in near future
investigations to understand these phenomena further.

##  V Acknowledgements

I gratefully acknowledge Dr. Santosh Kumar Das for his invaluable advice,
Aditi Tomar for insightful discussions and inspirational encouragement, and
Mohammad Yousuf Jamal for his numerous informative contributions, collectively
enriching the content of this paper. Additionally, I acknowledge the support
provided by the BRNS Project No.: 2021/EMR/SKD/026.

##  VI Appendix

Consider a partition $\\{t_{n}=\frac{nT}{N},;0\leq n\leq N\\}$ of time
interval $[0,T]$. The three-step numerical scheme called L2 approximation is
used for the case of superdiffusion. The second derivative $u^{(2)}$ is
approximated using a central difference formula, and the resulting numerical
scheme involves the values of $u$ at the previous three-time points $u^{n-1}$,
$u^{n-2}$, and $u^{n-3}$.

For $1<\nu\leq 2$, superdiffusion

${}^{C}D^{\nu}_{0+}u(t_{n})$
$\displaystyle=\frac{1}{\Gamma({2-\nu})}\int_{0}^{t_{n}}\frac{u^{(2)}(s)}{(t_{n}-s)^{\nu-1}}ds$
$\displaystyle=\frac{1}{\Gamma({2-\nu})}\sum_{j=1}^{n}\int_{t_{j-1}}^{t_{j}}\frac{u^{(2)}(s)}{(t_{n}-s)^{\nu-1}}ds$
$\displaystyle\approx\frac{1}{\Gamma({2-\nu})}\sum_{j=1}^{n}\frac{u(t_{j})-2u(t_{j-1})+u(t_{j-2})}{\Delta
t^{2}}\int_{t_{j-1}}^{t_{j}}\frac{1}{(t_{n}-s)^{\nu-1}}ds$
$\displaystyle=\frac{1}{\Gamma({3-\nu})}\sum_{j=1}^{n}\frac{((t_{n}-t_{j-1})^{2-\nu}-(t_{n}-t_{j})^{2-\nu})}{\Delta
t^{2}}(u(t_{j})-2u(t_{j-1})+u(t_{j-2}))$
$\displaystyle=\sum_{j=0}^{n-2}b_{j}(u(t_{n-j})-2u(t_{n-j-1})+u(t_{n-j-2}))$
$\displaystyle+\frac{n^{2-\nu}-(n-1)^{2-\nu}}{\Gamma(3-\nu)\Delta
t^{\nu}}(u(t_{1})-u(t_{0})).$ (29)

$\displaystyle^{C}D^{\nu}_{n}u(t_{n})=\begin{cases}\displaystyle\frac{u^{1}-u^{0}}{\Delta
t^{\nu}\Gamma(3-\nu)}\;&:\;n=1\\\
\displaystyle\frac{u^{2}-2u^{1}+u^{0}}{\Delta
t^{\nu}\Gamma(3-\nu)}\;&:\;n=2\\\
\displaystyle\frac{u^{n}-2u^{n-1}+u^{n-2}}{\Delta
t^{\nu}\Gamma(3-\nu)}+\sum_{j=1}^{n-2}b_{j}(u^{n-j}-2u^{n-j-1}+u^{n-j-2})&\\\
\displaystyle+\;\frac{n^{2-\nu}-(n-1)^{2-\nu}}{\Gamma(3-\nu)\Delta
t^{\nu}}(u^{1}-u^{0})\;&:\;n\geq 3\end{cases}$ (30)

where the coefficients
$b_{j}=\frac{(j+1)^{2-\nu}-j^{2-\nu}}{\Gamma(3-\nu)\Delta t^{\nu}}$ in the
scheme is determined by the difference formula for the second derivative and
are used to account for the fractional order.

## References

  * Adams  _et al._ (2006) J. Adams _et al._ (STAR), Phys. Rev. Lett. 97, 162301 (2006), arXiv:nucl-ex/0604018 .
  * Adams  _et al._ (2005) J. Adams _et al._ (STAR), Nucl. Phys. A 757, 102 (2005), arXiv:nucl-ex/0501009 .
  * Adcox  _et al._ (2005) K. Adcox _et al._ (PHENIX), Nucl. Phys. A 757, 184 (2005), arXiv:nucl-ex/0410003 .
  * Aamodt  _et al._ (2010) K. Aamodt _et al._ (ALICE), Phys. Rev. Lett. 105, 252301 (2010), arXiv:1011.3916 [nucl-ex] .
  * Arsene  _et al._ (2005) I. Arsene _et al._ (BRAHMS), Phys. Rev. C 72, 014908 (2005), arXiv:nucl-ex/0503010 .
  * van Hees and Rapp (2005) H. van Hees and R. Rapp, Phys. Rev. C 71, 034907 (2005), arXiv:nucl-th/0412015 .
  * Rapp and van Hees (2010) R. Rapp and H. van Hees (2010) pp. 111–206, arXiv:0903.1096 [hep-ph] .
  * Song  _et al._ (2015) T. Song, H. Berrehrah, D. Cabrera, J. M. Torres-Rincon, L. Tolos, W. Cassing, and E. Bratkovskaya, Phys. Rev. C 92, 014910 (2015), arXiv:1503.03039 [nucl-th] .
  * Andronic  _et al._ (2016) A. Andronic _et al._ , Eur. Phys. J. C 76, 107 (2016), arXiv:1506.03981 [nucl-ex] .
  * Dong and Greco (2019) X. Dong and V. Greco, Prog. Part. Nucl. Phys. 104, 97 (2019).
  * Cao  _et al._ (2019) S. Cao _et al._ , Phys. Rev. C 99, 054907 (2019), arXiv:1809.07894 [nucl-th] .
  * Beraudo  _et al._ (2018) A. Beraudo _et al._ , Nucl. Phys. A 979, 21 (2018), arXiv:1803.03824 [nucl-th] .
  * Prino and Rapp (2016) F. Prino and R. Rapp, J. Phys. G 43, 093002 (2016), arXiv:1603.00529 [nucl-ex] .
  * Aarts  _et al._ (2017) G. Aarts _et al._ , Eur. Phys. J. A 53, 93 (2017), arXiv:1612.08032 [nucl-th] .
  * Uphoff  _et al._ (2011) J. Uphoff, O. Fochler, Z. Xu, and C. Greiner, Phys. Rev. C 84, 024908 (2011), arXiv:1104.2295 [hep-ph] .
  * Golam Mustafa  _et al._ (1998) M. Golam Mustafa, D. Pal, and D. Kumar Srivastava, Phys. Rev. C 57, 889 (1998), [Erratum: Phys.Rev.C 57, 3499–3499 (1998)], arXiv:nucl-th/9706001 .
  * Plumari  _et al._ (2018) S. Plumari, V. Minissale, S. K. Das, G. Coci, and V. Greco, Eur. Phys. J. C 78, 348 (2018), arXiv:1712.00730 [hep-ph] .
  * Gossiaux and Aichelin (2008) P. B. Gossiaux and J. Aichelin, Phys. Rev. C 78, 014904 (2008), arXiv:0802.2525 [hep-ph] .
  * Prakash  _et al._ (2021) J. Prakash, M. Kurian, S. K. Das, and V. Chandra, Phys. Rev. D 103, 094009 (2021), arXiv:2102.07082 [hep-ph] .
  * Prakash  _et al._ (2023) J. Prakash, V. Chandra, and S. K. Das, Phys. Rev. D 108, 096016 (2023), arXiv:2306.07966 [hep-ph] .
  * Prakash and Jamal (2023a) J. Prakash and M. Y. Jamal, (2023a), arXiv:2304.04003 [nucl-th] .
  * Jamal  _et al._ (2023) M. Y. Jamal, J. Prakash, I. Nilima, and A. Bandyopadhyay, (2023), arXiv:2304.09851 [hep-ph] .
  * Singh  _et al._ (2023) M. Singh, M. Kurian, S. Jeon, and C. Gale, Phys. Rev. C 108, 054901 (2023), arXiv:2306.09514 [nucl-th] .
  * Kurian  _et al._ (2020) M. Kurian, M. Singh, V. Chandra, S. Jeon, and C. Gale, Phys. Rev. C 102, 044907 (2020), arXiv:2007.07705 [hep-ph] .
  * Cao  _et al._ (2016) S. Cao, T. Luo, G.-Y. Qin, and X.-N. Wang, Phys. Rev. C 94, 014909 (2016), arXiv:1605.06447 [nucl-th] .
  * Mazumder  _et al._ (2011) S. Mazumder, T. Bhattacharyya, J.-e. Alam, and S. K. Das, Phys. Rev. C 84, 044901 (2011), arXiv:1106.2615 [nucl-th] .
  * Zhang  _et al._ (2023) C. Zhang, L. Zheng, S. Shi, and Z.-W. Lin, Phys. Lett. B 846, 138219 (2023), arXiv:2210.07767 [nucl-th] .
  * Jamal  _et al._ (2021) M. Y. Jamal, S. K. Das, and M. Ruggieri, Phys. Rev. D 103, 054030 (2021).
  * Jamal and Mohanty (2021a) M. Y. Jamal and B. Mohanty, Eur. Phys. J. C 81, 616 (2021a), arXiv:2101.00164 [nucl-th] .
  * Jamal and Mohanty (2021b) M. Y. Jamal and B. Mohanty, Eur. Phys. J. Plus 136, 130 (2021b), arXiv:2002.09230 [nucl-th] .
  * Sun  _et al._ (2023) Y. Sun, S. Plumari, and S. K. Das, Phys. Lett. B 843, 138043 (2023), arXiv:2304.12792 [nucl-th] .
  * Plumari  _et al._ (2020) S. Plumari, G. Coci, V. Minissale, S. K. Das, Y. Sun, and V. Greco, Phys. Lett. B 805, 135460 (2020), arXiv:1912.09350 [hep-ph] .
  * Prakash and Jamal (2023b) J. Prakash and M. Y. Jamal, Journal of Physics G: Nuclear and Particle Physics 51, 025101 (2023b).
  * Debnath  _et al._ (2023) M. Debnath, R. Ghosh, M. Y. Jamal, M. Kurian, and J. Prakash, (2023), arXiv:2311.16005 [hep-ph] .
  * Du and Qian (2023) X. Du and W. Qian, arXiv preprint arXiv:2312.16294 (2023).
  * Shaikh  _et al._ (2021) A. Shaikh, M. Kurian, S. K. Das, V. Chandra, S. Dash, and B. K. Nandi, Phys. Rev. D 104, 034017 (2021), arXiv:2105.14296 [hep-ph] .
  * Kumar  _et al._ (2022) A. Kumar, M. Kurian, S. K. Das, and V. Chandra, Phys. Rev. C 105, 054903 (2022), arXiv:2111.07563 [hep-ph] .
  * Das  _et al._ (2022) S. K. Das _et al._ , International Journal of Modern Physics E , 2250097 (2022).
  * Kubo  _et al._ (2012) R. Kubo, M. Toda, and N. Hashitsume, 31 (2012).
  * Young  _et al._ (2012) C. Young, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 86, 034905 (2012).
  * Lang  _et al._ (2016) T. Lang, H. van Hees, G. Inghirami, J. Steinheimer, and M. Bleicher, Phys. Rev. C 93, 014901 (2016).
  * Cao and Bass (2011) S. Cao and S. A. Bass, Phys. Rev. C 84, 064902 (2011).
  * van Hees  _et al._ (2008) H. van Hees, M. Mannarelli, V. Greco, and R. Rapp, Phys. Rev. Lett. 100, 192301 (2008), arXiv:0709.2884 [hep-ph] .
  * Das  _et al._ (2014) S. K. Das, F. Scardina, S. Plumari, and V. Greco, Phys. Rev. C 90, 044901 (2014), arXiv:1312.6857 [nucl-th] .
  * He  _et al._ (2023) M. He, H. van Hees, and R. Rapp, Prog. Part. Nucl. Phys. 130, 104020 (2023), arXiv:2204.09299 [hep-ph] .
  * He  _et al._ (2012) M. He, R. J. Fries, and R. Rapp, Phys. Rev. C 86, 014903 (2012).
  * Das  _et al._ (2015) S. K. Das, F. Scardina, S. Plumari, and V. Greco, Phys. Lett. B 747, 260 (2015), arXiv:1502.03757 [nucl-th] .
  * Balescu (1995) R. Balescu, Physical Review E 51, 4807 (1995).
  * Bouchaud and Georges (1990) J.-P. Bouchaud and A. Georges, Physics reports 195, 127 (1990).
  * Metzler and Klafter (2000) R. Metzler and J. Klafter, Physics reports 339, 1 (2000).
  * Caspi  _et al._ (2000) A. Caspi, R. Granek, and M. Elbaum, Phys. Rev. Lett. 85, 5655 (2000).
  * Weber  _et al._ (2010) S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Phys. Rev. Lett. 104, 238102 (2010).
  * Vitali  _et al._ (2018) S. Vitali, V. Sposini, O. Sliusarenko, P. Paradisi, G. Castellani, and G. Pagnini, Journal of The Royal Society Interface 15, 20180282 (2018).
  * Oliveira  _et al._ (2019) F. A. Oliveira, R. M. Ferreira, L. C. Lapas, and M. H. Vainstein, Frontiers in Physics 7, 18 (2019).
  * Sokolov and Klafter (2005) I. M. Sokolov and J. Klafter, Chaos: An Interdisciplinary Journal of Nonlinear Science 15 (2005).
  * Lutz (2001) E. Lutz, Phys. Rev. E 64, 051106 (2001).
  * Lim and Teo (2009) S. Lim and L. Teo, Journal of Statistical Mechanics: Theory and Experiment 2009, P08015 (2009).
  * Wang and Chen (2023) X. Wang and Y. Chen, Phys. Rev. E 107, 024105 (2023).
  * Li  _et al._ (2012) C. Li, F. Zeng, and F. Liu, Fractional Calculus and Applied Analysis 15, 383 (2012).
  * Kobelev and Romanov (2000) V. Kobelev and E. Romanov, Progress of Theoretical Physics Supplement 139, 470 (2000).
  * Coffey  _et al._ (2004) W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron, World Scientific Publishing Co., Inc., River Edge, NJ 105 (2004), 10.1142/5343.
  * Gorenflo and Mainardi (1997) R. Gorenflo and F. Mainardi, CISM Courses and Lect., 378, 223 (1997).
  * Mainardi and Pironi (1996) F. Mainardi and P. Pironi, Extracta Math. 11, 140 (1996).
  * Mainardi (1997) F. Mainardi,  CISM Courses and Lect., 378, 291 (1997).
  * Guo  _et al._ (2013) P. Guo, C. Zeng, C. Li, and Y. Chen, Fractional Calculus and Applied Analysis 16, 123 (2013).
  * Li and Zeng (2013) C. Li and F. Zeng, Numerical Functional Analysis and Optimization 34, 149 (2013).
  * Miller and Ross (1993) K. S. Miller and B. Ross, John Wiley & Sons, Inc., New York , xvi+366 (1993).
  * Carpinteri and Mainardi (2014) A. Carpinteri and F. Mainardi, Fractals and Fractional Calculus in Continuum Mechanics 378 (2014).
  * Agrawal (2007) O. P. Agrawal, Journal of Physics A: Mathematical and Theoretical 40, 6287 (2007).
  * Feller (1971) W. Feller, John Wiley and Sons (1971).
  * Metzler  _et al._ (1999) R. Metzler, E. Barkai, and J. Klafter, Phys. Rev. Lett. 82, 3563 (1999).
  * Metzler (2000) R. Metzler, Phys. Rev. E 62, 6233 (2000).
  * Zumofen and Klafter (1994) G. Zumofen and J. Klafter, Chemical physics letters 219, 303 (1994).
  * Lomholt  _et al._ (2005) M. A. Lomholt, T. Ambjörnsson, and R. Metzler, Phys. Rev. Lett. 95, 260603 (2005).
  * Shlesinger  _et al._ (1993) M. F. Shlesinger, G. M. Zaslavsky, and J. Klafter, Nature 363, 31 (1993).
  * Wang and Tokuyama (1999) K. Wang and M. Tokuyama, Physica A: Statistical Mechanics and its Applications 265, 341 (1999).
  * Wang (1992a) K. Wang, Physical Review A 45, 833 (1992a).
  * Porrà  _et al._ (1996) J. M. Porrà, K.-G. Wang, and J. Masoliver, Phys. Rev. E 53, 5872 (1996).
  * Wang (1992b) K. G. Wang, Phys. Rev. A 45, 833 (1992b).
  * Fa (2006) K. S. Fa, Phys. Rev. E 73, 061104 (2006).
  * Siegle  _et al._ (2010) P. Siegle, I. Goychuk, and P. Hänggi, Phys. Rev. Lett. 105, 100602 (2010).
  * Chen  _et al._ (2023) W. Chen, C. Greiner, and Z. Xu, Phys. Rev. E 107, 064131 (2023).
  * Hammelmann  _et al._ (2019) J. Hammelmann, J. M. Torres-Rincon, J.-B. Rose, M. Greif, and H. Elfner, Phys. Rev. D 99, 076015 (2019).
  * Kapusta  _et al._ (2012) J. I. Kapusta, B. Müller, and M. Stephanov, Phys. Rev. C 85, 054906 (2012).
  * Ruggieri  _et al._ (2022) M. Ruggieri, Pooja, J. Prakash, and S. K. Das, Phys. Rev. D 106, 034032 (2022), arXiv:2203.06712 [hep-ph] .
  * Pooja  _et al._ (2023) Pooja, S. K. Das, V. Greco, and M. Ruggieri, Phys. Rev. D 108, 054026 (2023), arXiv:2306.13749 [hep-ph] .
  * Caucal and Mehtar-Tani (2022) P. Caucal and Y. Mehtar-Tani, Phys. Rev. D 106, L051501 (2022), arXiv:2109.12041 [hep-ph] .
  * Caputo (1967) M. Caputo, Geophysical Journal International 13, 529 (1967).
  * (89) I. Podlubny, Fractional Differential Equations, Aca- demic Press, San Diego, 1999 .
  * Sandev  _et al._ (2012) T. Sandev, R. Metzler, and . Tomovski, Fractional Calculus and Applied Analysis 15, 426 (2012).
  * Erdélyi  _et al._ (1953) A. Erdélyi, W. Magnus, F. Oberhettinger, and F. Tricomi, Higher transcendental functions 2, 133 (1953).
  * Hilfer (2000) R. Hilfer, World Scientific Publishing Co., Inc., River Edge, NJ , viii+463 (2000).
  * Rogers (1997) L. C. G. Rogers, Mathematical finance 7, 95 (1997).
  * Oldham and Spanier (1974) K. B. Oldham and J. Spanier, Mathematics in Science and Engineering 111 (1974), 2105.14296 .
  * Walton and Rafelski (2000) D. B. Walton and J. Rafelski, Phys. Rev. Lett. 84, 31 (2000).
  * Moore and Teaney (2005) G. D. Moore and D. Teaney, Phys. Rev. C 71, 064904 (2005), arXiv:hep-ph/0412346 .
  * Mazumder  _et al._ (2014) S. Mazumder, T. Bhattacharyya, and J.-e. Alam, Phys. Rev. D 89, 014002 (2014), arXiv:1305.6445 [nucl-th] .
  * Svetitsky (1988) B. Svetitsky, Phys. Rev. D 37, 2484 (1988).
  * Cacciari  _et al._ (2005) M. Cacciari, P. Nason, and R. Vogt, Phys. Rev. Lett. 95, 122001 (2005), arXiv:hep-ph/0502203 .
  * Cacciari  _et al._ (2012) M. Cacciari, S. Frixione, N. Houdeau, M. L. Mangano, P. Nason, and G. Ridolfi, JHEP 10, 137 (2012), arXiv:1205.6344 [hep-ph] .


Further author information: (Send correspondence to Jack Lashner)
Jack Lashner: E-mail<EMAIL_ADDRESS>
#  The Simons Observatory: A Minimum-Cost Matching Algorithm for Pairing
Measured Resonances with Designed Detectors

Jack Lashner  Wright Laboratory, Department of Physics, Yale University, New
Haven, Connecticut 06511  0000-0002-0001-2151 Kaiwen Zheng  Joseph Henry
Laboratories of Physics, Jadwin Hall, Princeton University, Princeton, NJ
08544  0000-0003-4645-7084 Kevin T. Crowley  Department of Physics, University
of California, San Diego, La Jolla, CA  0000-0001-5068-1295 Nicholas Galitzki
Department of Physics, University of Texas at Austin, Austin, TX, 78712, USA
Weinberg Institute for Theoretical Physics, Texas Center for Cosmology and
Astroparticle Physics, Austin, TX 78712, USA  0000-0001-7225-6679 Kathleen
Harrington  Argonne National Laboratory, High Energy Physics Division. 9700 S
Cass Ave, Lemont, IL 60439  University of Chicago, Department of Astronomy and
Astrophysics. 5801 S Ellis Ave, Chicago, IL 60637  0000-0003-1248-9563
Hironobu Nakata  Department of Physics, Faculty of Science, Kyoto University,
Kyoto  0000-0002-6300-1495 Max Silva-Feaver  Wright Laboratory, Department of
Physics, Yale University, New Haven, Connecticut 06511  0000-0001-7480-4341

###### Abstract

The Simons Observatory (SO) is a ground-based cosmic microwave background
experiment currently being deployed to Cerro Toco in the Atacama Desert of
Chile. The initial deployment of SO, consisting of three 0.46m-diameter small-
aperture telescopes and one 6m-primary large-aperture telescope, will field
over 60,000 transition-edge sensors that will observe at frequencies between
30 GHz and 280 GHz. SO will read out its detectors using Superconducting
Quantum Interference Device (SQUID) microwave-frequency multiplexing
($\mu$mux), a form of frequency division multiplexing where an RF-SQUID
couples each TES bolometer to a superconducting resonator tuned to a unique
frequency. Resonator frequencies are spaced roughly every 2 MHz between 4 and
6 GHz, allowing for multiplexing factors on the order of 1000. One challenge
of $\mu$mux is matching each tracked resonator with its corresponding physical
detector. Variations in resonator fabrication, and frequency shifts between
cooldowns caused by trapped flux can cause the measured resonance frequencies
to deviate significantly from their designed values. In this study, we
introduce a method for pairing measured and designed resonators by
constructing a bipartite graph based on the two resonator sets, and assigning
edge weights based on measured resonator and detector properties such as
resonance frequency, detector pointing, and assigned bias lines. Finding the
minimum-cost matching for a given set of edge weights is a well-studied
problem that can be solved very quickly, and this matching tells us the best
assignment of measured resonators to designed detectors for our input
parameters. We will present results based on the first on-sky measurements
from SAT1, the first SO MF small-aperture telescope.

##  1 Introduction

Millimeter wavelength radiation from the cosmic microwave background (CMB)
contains a wealth of information about the early universe, and is key to
understanding its structure and evolution. The Simons Observatory (SO) is a
collection of millimeter-wave telescopes, sited on Cerro Toco in the Atacama
Desert of Chile, which will observe the temperature and polarization of the
CMB in order to place constraints on cosmological parameters, and to map out
the large-scale structure of the universe [1, 2]. Initially, SO will consist
of three small-aperture telescopes[3] and one large-aperture telescope[4, 5],
which will be used to study the CMB from degree to arc-minute scales. SO will
use on order of 60,000 transition-edge sensing detectors (TES) to measure the
incoming radiation in six different optical bands between 30 and 280 GHz.

To read out the 60,000 detectors, SO uses SQUID microwave multiplexing
($\mu$mux), a form of frequency-division multiplexing which achieves
multiplexing factors on order of 1000 by inductively coupling each TES to a
superconducting resonator with resonance frequency between 4 and 6 GHz via an
RF-SQUID [6, 7]. Variations in the intensity of incoming radiation cause
shifts in the TES resistance that are transduced into changes in resonance
frequency. We use the SLAC Microresonator RF (SMURF) warm electronics[8] to
monitor and interpret these changes in resonance frequency, to convert them
into measurements of the TES current.

The focal planes for SO instruments consist of tiled, compact assemblies
called universal focal-plane modules (UFMs) that contain the TES detectors,
the optical coupling mechanism, and the cold readout electronics [9, 10]. Each
UFM contains two transmission lines, where each transmission line is coupled
to 924 superconducting resonators tuned to frequencies between 4-6 GHz. These
transmission lines can be used to read out up to 1756 detectors, 1720 with
optical coupling and 36 dark detectors. The remaining uncoupled resonators are
used for diagnostic purposes, such as to measure readout and two-level system
noise. The UFM also routes a number of DC bias lines, including two flux-ramp
lines and 12 TES bias lines.

Figure 1:  The complex transmission for a subsection of the 2 GHz band of the
transmission line. Resonators show up as dips in the transmission with width
around 30 kHz, and a designed spacing of 2 MHz.

During operation, we need to tune the readout electronics to make sure it is
outputting and tracking tones at the resonant frequencies of the
superconducting resonators. A tune, or a measurement of the resonance
frequencies, is created by measuring the complex transmission of a probe-tone
as it sweeps across the 2 GHz band, and using a peak-finding algorithm to
determine the location of each resonator. The complex transmission for a small
segment of the bandwidth can be seen in figure 1.

Calibration operations such as IVs are performed to determine resonator and
detector characteristics, such as which bias line is connected to the TES, the
normal resistance of the TES, and the operating bias voltages that should be
used for each bias line.

A critical step in the analysis pipeline is determining what physical detector
is connected to each tracked resonator. This mapping is required to compare
detector measurements across different tunes, and to combine data from
multiple detectors into a single map of the sky. The determination of this
mapping is non-trivial, since even though we know the intended resonator
frequencies of each detector, variations in the fabrication process and per-
cooldown shifts in resonator frequency due to flux trapped in the wires and
SQUIDs make it difficult to directly assign a measured resonance to a detector
without additional information. We must use a combination of resonator and
detector properties, such as the resonance frequency, the coupled bias line,
and the detector pointing information to determine the correct mapping from
measured resonators to designed detectors.

In the following sections we will describe the algorithm used to generate the
mapping between measured resonators and their physical detectors, and will
present results based on initial observations showing its efficacy.

##  2 Methodology

In order to create maps, we need to be able to accurately determine which
physical detector corresponds to each of the resonators that is being read
out. This determination is made based on a number of measured properties which
can be compared to the corresponding designed properties for each of the
detectors, including:

  * •

Resonance frequency of the superconducting resonator

  * •

Detector pointing angle relative to boresight center

  * •

TES bias line.

Previously, detector matches have been generated by hand based on lab
measurements of resonance frequency and bias line assignments, which is
challenging for a number of reasons. Because lab measurements do not provide
detector pointing information, these assignments are generated with
considerably less information than is available during deployment. Frequency
matching without additional information can be time-consuming because the
measured resonator frequencies scatter around the designed frequencies,
resulting in clusters and crosses. In some cases, due to resonator crossings
and missing resonators, it may be impossible to determine the correct match
without additional information. This method of manual mapping is not feasible
for deployment as it is time-consuming to generate for many UFMs, and we will
need to regularly create new matches as more calibration data are obtained.
Furthermore, the inclusion of pointing information in deployment greatly
reduces mapping uncertainty, enabling us to create an automated and robust
software solution for detector mapping.

To automate detector matching, we model the problem as a complete bipartite
graph where each resonator is a node, and the disjoint sets represent the set
of measured resonators and the set of designed ones. Edge weights between the
two sub-graphs are determined by a cost function that depends on resonator
properties such as those listed above. The best match is represented by the
minimum-cost matching of the graph, the finding of which is a well-studied
problem that can be solved very quickly for a given set of edge weights. Due
to the possibility that designed resonators are not present, or that spurious
peaks are mistaken for real resonators in the measured data, we need a way to
leave resonators on both sides unassigned. To do this, we add a fixed number
of unassigned nodes to each side with edge weights that penalize leaving a
resonator unassigned, where the penalty cost depends on the properties of the
unassigned resonator. A depiction of the bipartite graph and matching can be
seen in figure 2. Our implementation can be found in sotodlib, the SO library
for time-ordered data
analysis111https://github.com/simonsobs/sotodlib/blob/master/sotodlib/coords/det_match.py.

We generate our match based on edge weights using scipy’s
linear_sum_assignment
function222https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.linear_sum_assignment.html,
which implements a modified Jonker-Volgenant algorithm[11, 12]. For a full UFM
of around 1800 resonators per side, and an additional 1000 unassigned nodes
per side, this algorithm takes roughly 0.5 seconds to find the optimal match.

Figure 2:  The complete graph representation of the matching problem. The two
groups of nodes represent two sets of resonators to be matched, typically one
being the set of measured resonators and the other being the set of designed
detectors. Each colored node represents a different resonator and each grey
square represents an unassigned node, used to represent when a resonator is
not matched. Black and light grey lines show the resonator-to-node, and
resonator-to-unassigned-node edges in the graph respectively. The red lines
show an example matching between the two resonator sets.

The cost function that we use depends primarily on the three properties listed
above: the difference in the resonance frequencies, the difference in the
detector pointing angles, and whether the assigned bias lines based on
calibration measurements matches the designed bias line. Resonator frequencies
are determined by sweeping a probe-tone from 4-6 GHz on each transmission
line, measuring the complex transmission of the tone. A peak-finding algorithm
to determine the approximate location of each resonator, and then a finer
sweep about each peak is used to hone in on the resonance frequency, and cut
out resonators that don’t have the expected shape.

Pointing measurements are made by scanning a celestial source, typically the
Moon or Jupiter for the small-aperture telescopes, and finding the angle for
which each detector sees the source relative to the boresight of the
telescope. The designed pointing angles are derived using an optical model of
the instrument. Detector pointing is parametrized by the angles $\xi$ and
$\eta$: the angles of the detector beam center relative to the boresight in
the elevation and azimuth directions. They can be approximated as:

$\xi\approx\Delta\text{elevation}\qquad\eta\approx-\Delta\text{azimuth}/\cos(\text{elevation})$
(1)

Bias lines are assigned to each measured detector by stepping the bias voltage
up and down for each line one at a time while detectors are superconducting,
and measuring which bias line has the largest correlation with TES current for
each detector.

Below, we show the matching cost function and unassigned cost function that we
use to determine the edge weights of the graph:

⬇

1

2def matching_cost(r1, r2):

3 if r1.coax_line != r2.coax_line:

4 return np.inf

5

6 cost = np.exp(np.abs(r1.freq - r2.freq) / delta_freq)

7

8 theta_diff = np.sqrt((r1.xi - r2.xi)**2 + (r1.eta - r2.eta)**2)

9 cost += np.exp(theta_diff / delta_theta)

10

11 if r1.bias_line != r2.bias_line:

12 if r1.bias_line is None or r2.bias_line is None:

13 cost += bl_unassigned_mismatch_penalty

14 else:

15 cost += bl_assigned_mismatch_penalty

16 return cost

17

18def unassigned_cost(r):

19 if r.bias_line is None:

20 return bl_unassigned_unmatched_penalty

21 else:

22 return bl_assigned_unmatched_penalty

Here the matching costs depend on the resonator properties coax_line, the
transmission line that the resonator is connected to; freq, the resonance
frequency; bias_line, the bias line providing DC voltage bias to the TES; and
xi and eta, the pointing angles of the detector. The parameters delta_freq,
delta_theta, and the various bias line mismatch penalties are tunable
parameters that can be adjusted to change the behavior and priorities of the
matching algorithm. We force a matching cost of $\infty$ for resonators that
we know are on different transmission lines. Different penalties are used for
resonators where we were able to determine an assigned bias line, and
resonators where we were not able to determine a bias-line, since it is more
likely that measured resonators without proper bias-line assignments are
spurious peaks rather than real resonators.

The matching algorithm is typically run in one of two modes of operation. The
first is to find the best possible solution between measured resonators of a
UFM and designed detectors, called a solution match. This takes into account
any calibration data available, including pointing measurements generated from
scanning a celestial sources, bias line assignments. The solution match can be
thought of as a fingerprint that is unique to each UFM: though most UFMs share
the same designed resonator frequencies the fabrication variation is different
for every UFM, so the solution match is unique. For the solution match it is
important to run preprocessing on the measured data to account for any global
shifts and systematic offsets in both frequency and pointing space, as any
systematic offset in the datasets can easily lead to incorrect matches. The
methods that we use for data preprocessing and the matching parameters used to
generate the solution match are discussed in Appendix A.

The second mode of operation is to match a set of measured resonance
frequencies to an existing solution match. While the cooldown-to-cooldown
variation in resonance frequency is much smaller, on the order of tens of kHz,
this variation can still effect the SMuRF channel assignment in such a way
that it makes detector comparisons impossible across different tunes without a
matching. To account for this, we generate a new match to the measurement set
used to generate the solution match for each new tune, allowing us to
determine the designed detector for each measured resonator in the new tune.
Unlike the solution match, this mode of operation matches two measured
resonator sets with each other, and so there is only cooldown-to-cooldown
variation and no fabrication variation. This means we are able to rely just on
resonance frequencies for this match, and don’t need to incorporate other data
sets.

New tunes are created once at the beginning of each cooldown, and occasionally
during the cooldown if we observe higher noise measurements and suspect we can
benefit from a new measurement of the resonance frequencies. So far during the
deployment, we have not found it necessary to regularly create new tunes
during a cooldown, and tunes have remained stable for months at a time.

##  3 Results

Figure 3:  Matching results for the central UFM in SAT-MF1. (Left) Pointing
angles of the designed detectors (blue) and measured detectors (orange), and
the best matches between the two (black). Measured pointing angles are based
on beam maps made by scanning the Moon. Dark design detectors are depicted
with “x” markers. Due to the brightness of the Moon, the source can sometimes
be seen in the dark-detector data due to its thermal effect. Note that these
detector pointing angles have been preprocessed and aligned using the method
described in Appendix A. (Right) Resonator frequencies of the designed
detectors (top) and measured detectors (bottom), and the best matches between
them (light-grey lines in the middle) for a small region of frequency space.
Different colors of the resonator lines represent different assigned bias
lines, where grey means that the resonator is not coupled to a bias line,
either intentionally or due to wire-bonding issues. This region of frequency
space contains only two separate bias lines in addition to resonators with no
bias line assignments.

Here we show results from a single, central UFM in the SAT-MF1 instrument, as
this UFM has good detector yield and nearly complete pointing coverage. To
assess the performance of the matching algorithm and the accuracy of the
solution matches, we look at various metrics, including the percentage of
detectors which were matched, the percentage of resonators with matching bias
line assignments, and the differences in the pointing angles between the
measured and designed datasets. A depiction of differences between the two
detector sets in both pointing and frequency space can be seen in figure 3.
Statistics for these results can be seen in table 1. The distribution of the
differences in pointing angles between the measured and designed detectors can
be seen in figure 4.

After performing the data preprocessing described in Appendix A, detector
pointing angles align very well between the measured and designed datasets, as
can be seen in the left panel of figure 3. This allows us to match using a
small tolerance in pointing angle difference, resulting in fewer mismatches.
Detectors that are left unassigned are those for which we were not able to
make good pointing measurements, either due to detector data quality issues or
because the detectors were out of range of our source scans. We expect
pointing angles for these detectors to be filled in as we take and analyze
more calibration data, leading to better solutions.

Within a single pixel, we can be confident that we don’t have cross-
polarization or cross-band mis-matches, because resonators within a single
pixel are designed to have very different resonance frequencies. Additionally,
any given bias line only controls detectors observing at a single frequency
band, helping further avoid cross-band mismatches. If we believe that we still
see cross-band or cross-polarization mismatches, we may look into adding
additional datasets to the matching algorithm, such as wiregrid polarization
angle measurements, or detector time constants to break the degeneracy.

Figure 4:  Distribution of angular separation between measured and designed detector pointing angles for the best match. Before matching, deformations in the measured pointing angles are corrected using the procedure described in Appendix A. Angle differences of all detector pairs fall well below 0.5 degrees, which is the separation between neighboring pixels on the focal-plane. This signifies that there are few cross-pixel mismatches.  Designed Resonators | 1783
---|---
Measured Resonators | 1495
Matched Resonators | 1427
Table 1:  Table comparing the number of resonators in each set, the number of
matched resonators, and the number of resonators which have pointing
estimations that are left unmatched. Few detectors with good pointing
information are unable to find matches. For most unmatched detectors with
pointing information, even though we have pointing fits, these fits have high
uncertainties and are not suitable for a match with narrow pointing
tolerances. As we obtain better detector pointing data, we can expect the
percentage of matched detectors to increase. The large difference between the
number of measured and designed resonators is partially caused by two missing
mux-chips on this particular UFM, accounting for 132 missing resonators.

We can also validate detector matching results against datasets that were not
involved in the match creation process. One such dataset is a measurement of
our detector polarization angles measured using a rotating sparse wire grid
that reflects linearly polarized light along the wire direction[13]. We
compare the detector polarization angles measured using the wire grid
calibration with the angles of the designed detectors on the focal-plane and
see good agreement after subtracting a global angular offset from the measured
data, as can be seen in figure 5.

Figure 5:  The distribution in the difference between the detector
polarization angle measured with a wire grid. Here, we are taking only
detectors with good pointing information, as those matches are the most
reliable. We subtract out a global offset for each frequency band, to account
for a difference in definition of the polarization angle, and to account for
rotation that is induced by the detector time constant. We must subtract the
medians for the frequency bands separately because the detector time constants
are fairly different between the two frequency bands. Due to the arrangement
of the focal-plane, a mismatch between detectors that have different
polarization angles would appear at angles of $n\pi/12$ radians, which is
shown by the vertical dashed line for $n=1$. We can see a small cluster of
mismatched 90 GHz resonators, however 95% of detectors lie below the mismatch
threshold.

##  4 Conclusion

SO will be analyzing data from 60,000 detectors to generate detailed maps of
the CMB. SO will be using SQUID microwave multiplexing to readout on order of
1000 detectors per transmission line, by coupling each TES to a
superconducting resonator with a unique resonance frequency between 4 and 6
GHz. A critical component of the analysis is determining which physical
detector is mapped to each measured resonance frequency along the transmission
line. This is non-trivial due to differences in the design and true resonances
caused by fabrication variation as well as cooldown-to-cooldown shifts in
resonance frequency.

In this paper, we have presented our automated method for finding this mapping
using properties such as resonance frequency, TES bias line, and detector
pointing angles. We show matching statistics for the central UFM of the small-
aperture telescope. We show that, given good detector pointing data, we are
able to find detector matches for over 95% of measured resonators. We assess
the quality of the match by comparing the angular difference in measured and
designed detector pointing angles and detector polarization angles for matched
pairs. Based on these datasets, we find that mismatches are rare for detectors
with good pointing information. We expect that as we obtain more calibration
data, we will be able to generate correct match solutions for all UFMs.

##  5 Acknowledgements

This work was supported in part by a grant from the Simons Foundation (Award
#457687, B.K.).

## References

  * [1] Galitzki, N., Ali, A., Arnold, K. S., Ashton, P. C., Austermann, J. E., Baccigalupi, C., Baildon, T., Barron, D., Beall, J. A., Beckman, S., Bruno, S. M. M., Bryan, S., Calisse, P. G., Chesmore, G. E., Chinone, Y., Choi, S. K., Coppi, G., Crowley, K. D., Crowley, K. T., Cukierman, A., Devlin, M. J., Dicker, S., Dober, B., Duff, S. M., Dunkley, J., Fabbian, G., Gallardo, P. A., Gerbino, M., Goeckner-Wald, N., Golec, J. E., Gudmundsson, J. E., Healy, E. E., Henderson, S., Hill, C. A., Hilton, G. C., Ho, S.-P. P., Howe, L. A., Hubmayr, J., Jeong, O., Keating, B., Koopman, B. J., Kuichi, K., Kusaka, A., Lashner, J., Lee, A. T., Li, Y., Limon, M., Lungu, M., Matsuda, F., Mauskopf, P. D., May, A. J., McCallum, N., McMahon, J., Nati, F., Niemack, M. D., Orlowski-Scherer, J. L., Parshley, S. C., Piccirillo, L., Rao, M. S., Raum, C., Salatino, M., Seibert, J. S., Sierra, C., Silva-Feaver, M., Simon, S. M., Staggs, S. T., Stevens, J. R., Suzuki, A., Teply, G., Thornton, R., Tsai, C., Ullom, J. N., Vavagiakis, E. M., Vissers, M. R., Westbrook, B., Wollack, E. J., Xu, Z., and Zhu, N., “The Simons Observatory: Instrument Overview,” Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX , 3 (July 2018).
  * [2] Ade, P., Aguirre, J., Ahmed, Z., Aiola, S., Ali, A., Alonso, D., Alvarez, M. A., Arnold, K., Ashton, P., Austermann, J., Awan, H., Baccigalupi, C., Baildon, T., Barron, D., Battaglia, N., Battye, R., Baxter, E., Bazarko, A., Beall, J. A., Bean, R., Beck, D., Beckman, S., Beringue, B., Bianchini, F., Boada, S., Boettger, D., Bond, J. R., Borrill, J., Brown, M. L., Bruno, S. M., Bryan, S., Calabrese, E., Calafut, V., Calisse, P., Carron, J., Challinor, A., Chesmore, G., Chinone, Y., Chluba, J., Cho, H.-M. S., Choi, S., Coppi, G., Cothard, N. F., Coughlin, K., Crichton, D., Crowley, K. D., Crowley, K. T., Cukierman, A., D’Ewart, J. M., Dünner, R., de Haan, T., Devlin, M., Dicker, S., Didier, J., Dobbs, M., Dober, B., Duell, C. J., Duff, S., Duivenvoorden, A., Dunkley, J., Dusatko, J., Errard, J., Fabbian, G., Feeney, S., Ferraro, S., Fluxà, P., Freese, K., Frisch, J. C., Frolov, A., Fuller, G., Fuzia, B., Galitzki, N., Gallardo, P. A., Ghersi, J. T. G., Gao, J., Gawiser, E., Gerbino, M., Gluscevic, V., Goeckner-Wald, N., Golec, J., Gordon, S., Gralla, M., Green, D., Grigorian, A., Groh, J., Groppi, C., Guan, Y., Gudmundsson, J. E., Han, D., Hargrave, P., Hasegawa, M., Hasselfield, M., Hattori, M., Haynes, V., Hazumi, M., He, Y., Healy, E., Henderson, S. W., Hervias-Caimapo, C., Hill, C. A., Hill, J. C., Hilton, G., Hilton, M., Hincks, A. D., Hinshaw, G., Hložek, R., Ho, S., Ho, S.-P. P., Howe, L., Huang, Z., Hubmayr, J., Huffenberger, K., Hughes, J. P., Ijjas, A., Ikape, M., Irwin, K., Jaffe, A. H., Jain, B., Jeong, O., Kaneko, D., Karpel, E. D., Katayama, N., Keating, B., Kernasovskiy, S. S., Keskitalo, R., Kisner, T., Kiuchi, K., Klein, J., Knowles, K., Koopman, B., Kosowsky, A., Krachmalnicoff, N., Kuenstner, S. E., Kuo, C.-L., Kusaka, A., Lashner, J., Lee, A., Lee, E., Leon, D., Leung, J. S.-Y., Lewis, A., Li, Y., Li, Z., Limon, M., Linder, E., Lopez-Caraballo, C., Louis, T., Lowry, L., Lungu, M., Madhavacheril, M., Mak, D., Maldonado, F., Mani, H., Mates, B., Matsuda, F., Maurin, L., Mauskopf, P., May, A., McCallum, N., McKenney, C., McMahon, J., Meerburg, P. D., Meyers, J., Miller, A., Mirmelstein, M., Moodley, K., Munchmeyer, M., Munson, C., Naess, S., Nati, F., Navaroli, M., Newburgh, L., Nguyen, H. N., Niemack, M., Nishino, H., Orlowski-Scherer, J., Page, L., Partridge, B., Peloton, J., Perrotta, F., Piccirillo, L., Pisano, G., Poletti, D., Puddu, R., Puglisi, G., Raum, C., Reichardt, C. L., Remazeilles, M., Rephaeli, Y., Riechers, D., Rojas, F., Roy, A., Sadeh, S., Sakurai, Y., Salatino, M., Rao, M. S., Schaan, E., Schmittfull, M., Sehgal, N., Seibert, J., Seljak, U., Sherwin, B., Shimon, M., Sierra, C., Sievers, J., Sikhosana, P., Silva-Feaver, M., Simon, S. M., Sinclair, A., Siritanasak, P., Smith, K., Smith, S. R., Spergel, D., Staggs, S. T., Stein, G., Stevens, J. R., Stompor, R., Suzuki, A., Tajima, O., Takakura, S., Teply, G., Thomas, D. B., Thorne, B., Thornton, R., Trac, H., Tsai, C., Tucker, C., Ullom, J., Vagnozzi, S., van Engelen, A., Lanen, J. V., Winkle, D. D. V., Vavagiakis, E. M., Vergès, C., Vissers, M., Wagoner, K., Walker, S., Ward, J., Westbrook, B., Whitehorn, N., Williams, J., Williams, J., Wollack, E. J., Xu, Z., Yu, B., Yu, C., Zago, F., Zhang, H., and and, N. Z., “The Simons Observatory: Science goals and forecasts,” Journal of Cosmology and Astroparticle Physics 2019, 056–056 (Feb. 2019).
  * [3] Galitzki, N., Tsan, T., Spisak, J., Randall, M., Silva-Feaver, M., Seibert, J., Lashner, J., Adachi, S., Adkins, S. M., Alford, T., Arnold, K., Ashton, P. C., Austermann, J. E., Baccigalupi, C., Bazarko, A., Beall, J. A., Bhimani, S., Bixler, B., Coppi, G., Corbett, L., Crowley, K. D., Crowley, K. T., Day-Weiss, S., Dicker, S., Dow, P. N., Duell, C. J., Duff, S. M., Gerras, R. G., Groh, J. C., Gudmundsson, J. E., Harrington, K., Hasegawa, M., Healy, E., Henderson, S. W., Hubmayr, J., Iuliano, J., Johnson, B. R., Keating, B., Keller, B., Kiuchi, K., Kofman, A. M., Koopman, B. J., Kusaka, A., Lee, A. T., Lew, R. A., Lin, L. T., Link, M. J., Lucas, T. J., Lungu, M., Mangu, A., McMahon, J. J., Miller, A. D., Moore, J. E., Morshed, M., Nakata, H., Nati, F., Newburgh, L. B., Nguyen, D. V., Niemack, M. D., Page, L. A., Sakaguri, K., Sakurai, Y., Rao, M. S., Saunders, L. J., Shroyer, J. E., Sugiyama, J., Tajima, O., Takeuchi, A., Bua, R. T., Teply, G., Terasaki, T., Ullom, J. N., Van Lanen, J. L., Vavagiakis, E. M., Vissers, M. R., Walters, L., Wang, Y., Xu, Z., Yamada, K., and Zheng, K., “The Simons Observatory: Design, integration, and testing of the small aperture telescopes,” (May 2024).
  * [4] Zhu, N., Bhandarkar, T., Coppi, G., Kofman, A. M., Orlowski-Scherer, J. L., Xu, Z., Adachi, S., Ade, P., Aiola, S., Austermann, J., Bazarko, A. O., Beall, J. A., Bhimani, S., Bond, J. R., Chesmore, G. E., Choi, S. K., Connors, J., Cothard, N. F., Devlin, M., Dicker, S., Dober, B., Duell, C. J., Duff, S. M., Dünner, R., Fabbian, G., Galitzki, N., Gallardo, P. A., Golec, J. E., Haridas, S. K., Harrington, K., Healy, E., Ho, S.-P. P., Huber, Z. B., Hubmayr, J., Iuliano, J., Johnson, B. R., Keating, B., Kiuchi, K., Koopman, B. J., Lashner, J., Lee, A. T., Li, Y., Limon, M., Link, M., Lucas, T. J., McCarrick, H., Moore, J., Nati, F., Newburgh, L. B., Niemack, M. D., Pierpaoli, E., Randall, M. J., Sarmiento, K. P., Saunders, L. J., Seibert, J., Sierra, C., Sonka, R., Spisak, J., Sutariya, S., Tajima, O., Teply, G. P., Thornton, R. J., Tsan, T., Tucker, C., Ullom, J., Vavagiakis, E. M., Vissers, M. R., Walker, S., Westbrook, B., Wollack, E. J., and Zannoni, M., “The Simons Observatory Large Aperture Telescope Receiver,” The Astrophysical Journal Supplement Series 256, 23 (Sept. 2021).
  * [5] Gudmundsson, J. E., Gallardo, P. A., Puddu, R., Dicker, S. R., Adler, A. E., Ali, A. M., Bazarko, A., Chesmore, G. E., Coppi, G., Cothard, N. F., Dachlythra, N., Devlin, M., Dünner, R., Fabbian, G., Galitzki, N., Golec, J. E., Ho, S.-P. P., Hargrave, P. C., Kofman, A. M., Lee, A. T., Lee, A. T., Limon, M., Matsuda, F. T., Mauskopf, P. D., Moodley, K., Moodley, K., Nati, F., Niemack, M. D., Niemack, M. D., Orlowski-Scherer, J., Page, L. A., Partridge, B., Puglisi, G., Reichardt, C. L., Sierra, C. E., Simon, S. M., Teply, G. P., Tucker, C., Wollack, E. J., Xu, Z., and Zhu, N., “The Simons Observatory: Modeling optical systematics in the Large Aperture Telescope,” Applied Optics 60, 823–837 (Feb. 2021).
  * [6] Dober, B., Becker, D. T., Bennett, D. A., Bryan, S. A., Duff, S. M., Gard, J. D., Hays-Wehle, J. P., Hilton, G. C., Hubmayr, J., Mates, J. A. B., Reintsema, C. D., Vale, L. R., and Ullom, J. N., “Microwave SQUID Multiplexer demonstration for Cosmic Microwave Background Imagers,” Applied Physics Letters 111, 243510 (Dec. 2017).
  * [7] Henderson, S. W., Ahmed, Z., Austermann, J., Becker, D., Bennett, D. A., Brown, D., Chaudhuri, S., Cho, H.-M. S., D’Ewart, J. M., Dober, B., Duff, S. M., Dusatko, J. E., Fatigoni, S., Frisch, J. C., Gard, J. D., Halpern, M., Hilton, G. C., Hubmayr, J., Irwin, K. D., Karpel, E. D., Kernasovskiy, S. S., Kuenstner, S. E., Kuo, C.-L., Li, D., Mates, J. A. B., Reintsema, C. D., Smith, S. R., Ullom, J., Vale, L. R., Van Winkle, D. D., Vissers, M., and Yu, C., “Highly-multiplexed microwave SQUID readout using the SLAC Microresonator Radio Frequency (SMuRF) Electronics for Future CMB and Sub-millimeter Surveys,” Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX , 43 (July 2018).
  * [8] Yu, C., Ahmed, Z., Frisch, J. C., Henderson, S. W., Silva-Feaver, M., Brown, D., Connors, J., Cukierman, A. J., D’Ewart, J. M., Dober, B. J., Haller, G., Herbst, R., Hilton, G. C., Hubmayr, J., Irwin, K. D., Mates, J. A. B., Ruckman, L., Ullom, J., Vale, L., Winkle, D. D. V., Vasquez, J., and Young, E., “SLAC Microresonator RF (SMuRF) Electronics: A tone-tracking readout system for superconducting microwave resonator arrays,” 29.
  * [9] McCarrick, H., Healy, E., Ahmed, Z., Arnold, K., Atkins, Z., Austermann, J. E., Bhandarkar, T., Beall, J. A., Bruno, S. M., Choi, S. K., Connors, J., Cothard, N. F., Crowley, K. D., Dicker, S., Dober, B., Duell, C. J., Duff, S. M., Dutcher, D., Frisch, J. C., Galitzki, N., Gralla, M. B., Gudmundsson, J. E., Henderson, S. W., Hilton, G. C., Ho, S.-P. P., Huber, Z. B., Hubmayr, J., Iuliano, J., Johnson, B. R., Kofman, A. M., Kusaka, A., Lashner, J., Lee, A. T., Li, Y., Link, M. J., Lucas, T. J., Lungu, M., Mates, J. A. B., McMahon, J. J., Niemack, M. D., Orlowski-Scherer, J., Seibert, J., Silva-Feaver, M., Simon, S. M., Staggs, S., Suzuki, A., Terasaki, T., Thornton, R., Ullom, J. N., Vavagiakis, E. M., Vale, L. R., Lanen, J. V., Vissers, M. R., Wang, Y., Wollack, E. J., Xu, Z., Young, E., Yu, C., Zheng, K., and Zhu, N., “The Simons Observatory Microwave SQUID Multiplexing Detector Module Design,” The Astrophysical Journal 922, 38 (Nov. 2021).
  * [10] Wang, Y., Zheng, K., Atkins, Z., Austermann, J., Bhandarkar, T., Choi, S. K., Duff, S. M., Dutcher, D., Galitzki, N., Healy, E., Huber, Z. B., Hubmayr, J., Johnson, B. R., Lashner, J., Li, Y., McCarrick, H., Niemack, M. D., Seibert, J., Staggs, S. T., Vavagiakis, E., and Xu, Z., “Simons Observatory Focal-Plane Module: In-lab Testing and Characterization Program,” arXiv:2111.11301 [astro-ph, physics:physics] (Nov. 2021).
  * [11] Ramshaw, L. and Tarjan, R. E., “On Minimum-Cost Assignments in Unbalanced Bipartite Graphs,” (June 2012).
  * [12] Crouse, D. F., “On implementing 2D rectangular assignment algorithms,” IEEE Transactions on Aerospace and Electronic Systems 52, 1679–1696 (Aug. 2016).
  * [13] Murata, M., Nakata, H., Iijima, K., Adachi, S., Seino, Y., Kiuchi, K., Matsuda, F., Randall, M. J., Arnold, K., Galitzki, N., Johnson, B. R., Keating, B., Kusaka, A., Lloyd, J. B., Seibert, J., Silva-Feaver, M., Tajima, O., Terasaki, T., and Yamada, K., “The Simons Observatory: A fully remote controlled calibration system with a sparse wire grid for cosmic microwave background telescopes,” Review of Scientific Instruments 94, 124502 (Dec. 2023).

##  Appendix A Data Preprocessing

When generating a solution match, it is important to preprocess the data
before running the matching algorithm, as systematic shifts in the resonance
frequency or detector pointing datasets can easily lead to incorrect matches.
In addition to random variation in the resonance frequencies, fabrication
error can cause entire mux-chips, containing hundreds of resonators clumped
together in frequency space, to be offset by a constant frequency. Resonance
frequencies can also shift between cooldowns due to changes in the amount of
flux trapped in the transmission line and SQUID loops. For detector pointing,
the measured pointing angles and designed pointing angles may be off by an
affine transformation. To correct for these transformations, we apply a
sequence of iterative matches, freeing up different axes of parameter space in
each iteration by adjusting the tunable parameters, and using the results to
correct the measured data. Once resonator frequency and pointing angles are
correctly aligned, we can narrow tolerances in the cost function to produce a
final solution match.

First, we apply a global shift to the measured detector pointing angles, to
account for global systematic offsets. In most cases, the shift can be found
by finding the difference between the centroid of the measured and designed
pointing angle clusters. However, in cases where significant portions of the
focal-plane are missing measured pointing data, we may need to set this global
offset by eye.

Second, in order to correct for global and per-mux-chip frequency offsets, we
greatly relax the resonance-frequency penalty and apply a rough match based on
pointing angles and bias line assignments. Here we tune $\Delta\theta$ in the
listing 2 to be several degrees, much larger than the difference in angle
between neighboring detectors, which is roughly 0.5 degrees. We subtract out a
rolling median of the frequency error between resonators and their match over
large chunks of frequency space (we use a box size of 100 MHz) to correct for
global offsets, as can be seen in figure 6.

Figure 6:  Resonator frequency offset between measured and designed resonators
for a region of frequency space before any corrections are applied. The orange
line is the rolling median of the offset using a box-size of 100 MHz. This is
what is subtracted by the first correction step to roughly align measured and
design frequencies.

Next, we use a similar method to correct for large-scale pointing
deformations. We perform a match based primarily on resonance frequency and
bias line assignment, relaxing the pointing restriction such that it is much
larger than the difference in angle between neighboring detectors. We then
compute the median pointing error between a resonator and its match over large
swaths of pointing space (we use a radius of 2 degrees), and subtract it out
to correct for local pointing deformations. The local deformations that are
applied to the measured detector pointing before the final match can be seen
in figure 7. This process corrects for residual pointing offsets that are not
accounted for by the global shift, including any small scaling and rotational
transformations. Corrected resonator frequencies and detector angles can from
this process can be seen in figure 3, where we see that the datasets are now
well-aligned.

Figure 7:  Detector pointing corrections between measured and designed
detectors. The detector pointing angles shown here are those after an initial
global correction that aligns the centers of the two distributions. The black
arrows show the local interpolated pointing offsets for the pre-corrected
match, computed by taking the median in the pointing angle difference between
design and measured detectors about a circle centered at each point with a
radius of 2 degrees. These are the final pointing corrections that are applied
before the final match, shown in figure 3.