Heavy quarkonia in a baryon asymmetric strongly magnetized hot quark matter

Salman Ahamad<EMAIL_ADDRESS>Mujeeb
<EMAIL_ADDRESS>and Binoy Krishna Patra†<EMAIL_ADDRESS>† Department of Physics, Indian Institute of Technology Roorkee, Roorkee 247
667, India

###### Abstract

Recently there is a resurrection in the study of heavy quark bound states in a
hot and baryonless matter with an ambient magnetic field but the matter
produced at heavy-ion collider experiments is not perfectly baryonless, so we
wish to explore the effect of small baryon asymmetry on the in-medium
properties of the heavy quark bound states submerged in a strongly magnetized
hot quark matter. Therefore, we have first given a revisit to the general
covariant tensor structure of gluon self-energy in the above environment to
compute the resummed propagator for gluons. This resummed propagator embodies
the properties of medium, which gets translated into the (complex) potential
between $Q$ and $\bar{Q}$ placed in the medium. We observe that the baryon
asymmetry makes the real-part of potential slightly more attractive and
weakens the imaginary-part. This opposing effects thus lead to the enhancement
of binding energies and the reduction of thermal widths of $Q\bar{Q}$ ground
states, respectively. Finally, the properties of quarkonia thus deciphered
facilitate to compute the dissociation points of $J/\psi$ and $\Upsilon$,
which are found to have slightly larger values in the presence of baryon
asymmetry. For example, $J/\psi$ gets dissociated at $1.64~{}T_{c}$,
$1.69~{}T_{c}$, and $1.75~{}T_{c}$, whereas $\Upsilon$ is dissociated at
$1.95~{}T_{c}$, $1.97~{}T_{c}$ and $2.00~{}T_{c}$, for $\mu=0,60$ and $100$
MeV, respectively. This observation prevents early dissociation of quarkonia
in the matter produced at ultrarelativistic heavy ion collisions with a small
net baryon number, compared to the ideal baryonless matter.

Keywords: Strong magnetic field; Perturbative QCD; Quark chemical potential;
Resummed propagator; Heavy quark potential;

##  1 Introduction

In the presence of high temperatures and/or densities, colourless hadrons get
melted into a plasma of asymptotically free quarks and gluons known as Quark-
Gluon Plasma (QGP). This type of extreme environment is created in the
ultrarelativistic heavy-ion collisions (URHICs) experiments at RHIC [1, 2] and
LHC [3, 4] and is planned to be created in the Compressed Baryonic Matter
(CBM) experiment at Facility for Antiproton and Ion Research (FAIR) [5]. In
noncentral events in URHICs, the relative motion of the spectator quarks
generates a strong magnetic field (SMF) at initial phases of the collisions.
The magnitude of the produced magnetic field ($B$) is around $m_{\pi}^{2}$
($\sim 10^{18}$ Gauss) at RHIC to 10 $m_{\pi}^{2}$ at LHC [6, 7]. This strong
magnetic field is a short pulse. Earlier, people did not pay much attention to
the effects of magnetic field on thermal QCD medium because they thought that
the life-span of the strong magnetic field was too short to have any
observable effects in the phenomenology of the heavy-ion collisions. They
observed that the produced strong magnetic field decays very fast in vacuum.
But some recent theoretical calculations predict that a thermal medium could
have produced as early as the production of magnetic field. As a consequence,
the life span of the magnetic field is elongated due to the finite electrical
conductivity of the medium [8, 9, 10]. The decay of the magnetic field in
medium becomes slower and it remains strong for much larger time. However, its
magnitude becomes weaker as the time is elapsed. Many theoretical studies have
been conducted in recent years to decode the effects of this intense $B$ on
thermal and magnetic properties of QGP [11, 12, 13, 14] and chiral and axial
magnetic effects [15, 16, 17, 18].

The timescale for the generation of the SMF and the heavy quark pairs in
URHICs is almost similar. Apart from this, it was found in few studies that a
large quark chemical potential ($\mu$) (upto 100 MeV) is also produced near
the deconfining temperature [19, 20, 21] and in the SMF its strength reaches
upto 200 MeV [22]. The experimental set up at FAIR will provide facilities to
investigate the rich physics of the deconfined phase of matter at high baryon
densities. The study of the effect of baryon asymmetry on the properties of
the QCD medium is very interesting and have applications to few branches of
cosmology, the physics of the core of neutron stars, and URHICs. The
$Q\bar{Q}$ bound state of the heavy quarks is one of the very promising
signatures to visualize the properties of QGP, therefore the study of the
effects of baryon asymmetry on $Q\bar{Q}$ potential becomes a necessary task.
Two of us recently investigated the $Q\bar{Q}$ potential and dissociation
temperatures for the bound states in magnetized QCD medium at zero chemical
potential [23]. In present work, we will examine the effects of the finite
quark chemical potential on the various properties (binding energy, decay
width and dissociation temperature) of $Q\bar{Q}$ bound states.

In past few years, heavy quarkonium physics has seen many developments. To
describe the heavy quark bound system, many effective field theories (EFTs)
have been derived using the hierarchy of the various scales of $Q\bar{Q}$
system. One of such theory is non-relativistic QCD (NRQCD) [24] which has been
formulated by integrating out the mass ($m_{Q}$) . Similarly, potential NRQCD
[25], is derived by integrating out the momentum exchange ($m_{Q}v$). Since
the hierarchy of scales in these effective theories is not very noticeable,
lattice QCD approach [26] is formulated to overcome the difficulties in EFTs.
In parallel, phenomenological potential models have been also derived which
provide an substitute to probe the various properties of the quarkonia. It has
been observed in [27] that the $Q\bar{Q}$ potential possess both real as well
as imaginary part. The real part gets screened because of the color charges
[28] whereas the resonance acquires a thermal width as a result of the
imaginary part [29]. Earlier, it was believed that the dissociation of the
quarkonia was because of the color screening phenomena but now a days it is
well understood that the $Q\bar{Q}$ bound states are dissolved mainly because
of the widening of the thermal width either due to Landau damping [27] or
color singlet-octet transition [30]. The bound state gets dissociated at lower
temperatures (in comparison to the binding energy) even if there is very weak
screening present in the medium as a consequence of the imaginary part. One of
us has studied the medium modification to imaginary-part by taking both the
perturbative and non-perturbative components and has observed that the
imaginary part becomes smaller [31, 32], compared to the perturbative term
alone [33]. The imaginary part has been also calculated in the framework of
Gauge-gravity duality in [34, 35] and in lattice studies [36]. The (complex)
$Q\bar{Q}$ potential has been explored using the generalized Gauss law in [37]
and the divergence in the imaginary-part was tamed by choosing the Debye mass
as the regulation scale. In a recent work, the string part of the potential
has been calculated through a phenomenological term which includes the effects
due to the low frequency modes subsumed in the dimension two gluon condensates
[38]. The above mentioned works have been carried out in the absence of the
$B$. The effect of the SMF on the heavy quarkonia has been investigated
recently in many groups. The influence of constant external magnetic field on
the QCD bound states have been studied for the case of a harmonic interaction
and for Cornell potential plus a spin spin interaction term in [39, 40]. The
$Q\bar{Q}$ bound states have been extensively investigated in SMF in [41, 42].
Two of us have explored the effects of SMF on the in-medium properties of the
$Q\bar{Q}$ bound states in hot medium and studied the dissociation through
both color screening and Landau damping phenomenon [43, 23]. Authors in [44]
have also investigated the complex $Q\bar{Q}$ potential using the generalized
Gauss law. An attempt has been made to study the anisotropic nature of the
inter-quark potential in [45] and the dissociation of the heavy quark bound
states in the weak magnetic field in [46]. Apart from theses works on the
$Q\bar{Q}$ bound states, the effects of $B$ on the propagation of the heavy
quarks in the thermal QCD medium have been investigated in some other works
like wakes in [47]. The collisional energy loss in [48] and the anisotropic
nature of the diffusion and the drag coefficients using the Langevindynamics
in [49]. In addition to this the general structure of the gauge boson two
point functions have been investigated in the magnetized hot material medium
in [50].

In this article, we will study the implications of the baryon asymmetry on the
in-medium behaviour of $Q\bar{Q}$ bound states submerged in the strongly
magnetized hot QCD medium. We have reassessed the structure of the gluon two
point functions in magnetized thermal medium. We, then evaluate the gluon
self-energy and relevant form factors in the frame work of the imaginary-time
formalism. Using these form factors, we have calculated the real and imaginary
parts of the resummed (full) gluon propagator which will later facilitate the
calculation of the (complex) $Q\bar{Q}$ potential. The real-part is plugged
into the Schrödinger equation to get the binding energy whereas the imaginary-
part gives the thermal width. At last, we study the effect of quark chemical
potential on the quasi-free dissociation of the heavy quarkonia by calculating
the dissociation temperatures for the $J/\psi$ and $\Upsilon$ states.

The paper has been organized as follows: In section 2, we have revisited the
general covariant tensor structure of the gluon self-energy in the extreme
environment of temperature, density and SMF. In subsection 2.1, we calculate
the real and imaginary parts of the form factor $b(p_{0},p)$ in imaginary-
time formalism which give the complex resummed gluon propagator in subsection
3.1. Next, we calculate the medium modification to the complex $Q\bar{Q}$
potential in subsection 3.2. We will discuss the binding energy in subsection
4.1, decay width in 4.2 and will explore the quasi-free dissociation process
of $J/\psi$ and $\Upsilon$ in subsection 4.3. Finally, we present the
conclusion of this study in section 5.

##  2  Gluon self-energy tensor and resummed gluon propagator in hot and dense
strongly magnetized medium

In this section, we will review the covariant tensor structure of gluon self-
energy in magnetized hot and dense QCD medium. The covariant structure is
given by [50]

$\displaystyle\Pi^{\mu\nu}(P)=b(P)B^{\mu\nu}(P)+c(P)R^{\mu\nu}(P)+d(P)M^{\mu\nu}(P)+a(P)N^{\mu\nu}(P),$
(1)

the projection tensors used in the construction of self-energy tensor in Eq.
(1) are defined as

$\displaystyle B^{\mu\nu}(P)$ $\displaystyle=$
$\displaystyle\frac{{\bar{u}}^{\mu}{\bar{u}}^{\nu}}{{\bar{u}}^{2}},$ (2)
$\displaystyle R^{\mu\nu}(P)$ $\displaystyle=$ $\displaystyle
g_{\perp}^{\mu\nu}-\frac{P_{\perp}^{\mu}P_{\perp}^{\nu}}{P_{\perp}^{2}},$ (3)
$\displaystyle M^{\mu\nu}(P)$ $\displaystyle=$
$\displaystyle\frac{{\bar{n}}^{\mu}{\bar{n}}^{\nu}}{{\bar{n}}^{2}},$ (4)
$\displaystyle N^{\mu\nu}(P)$ $\displaystyle=$
$\displaystyle\frac{{\bar{u}}^{\mu}{\bar{n}}^{\nu}+{\bar{u}}^{\nu}{\bar{n}}^{\mu}}{\sqrt{{\bar{u}}^{2}}\sqrt{{\bar{n}}^{2}}},$
(5)

where $u^{\mu}=(1,0,0,0)$ is the four velocity of the heat bath and
$n_{\mu}=(0,0,0,1)$ represents the direction of $B$. ${\bar{u}}^{\mu}$ and
${\bar{n}}^{\mu}$ are constructed as

$\displaystyle\bar{u}^{\mu}$ $\displaystyle=$
$\displaystyle\left(g^{\mu\nu}-\frac{P^{\mu}P^{\nu}}{P^{2}}\right)u_{\nu},$
(6) $\displaystyle\bar{n}^{\mu}$ $\displaystyle=$
$\displaystyle\left(\tilde{g}^{\mu\nu}-\frac{\tilde{P}^{\mu}\tilde{P}^{\nu}}{\tilde{P}^{2}}\right)n_{\nu},$
(7)

where ${\tilde{g}}^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ and
$\tilde{P}^{\mu}=P^{\mu}-(P.u)~{}u^{\mu}$. The form factors defined in (1) can
be evaluated using the contraction properties as

$\displaystyle b(P)$ $\displaystyle=$ $\displaystyle
B^{\mu\nu}(P)~{}\Pi_{\mu\nu}(P),$ (8) $\displaystyle c(P)$ $\displaystyle=$
$\displaystyle R^{\mu\nu}(P)~{}\Pi_{\mu\nu}(P),$ (9) $\displaystyle d(P)$
$\displaystyle=$ $\displaystyle M^{\mu\nu}(P)~{}\Pi_{\mu\nu}(P),$ (10)
$\displaystyle a(P)$ $\displaystyle=$
$\displaystyle\frac{1}{2}N^{\mu\nu}(P)~{}\Pi_{\mu\nu}(P).$ (11)

The general covariant structure of the full gauge boson (gluon) propagator in
magnetized hot and dense medium can be written in Landau gauge as [50]

$\displaystyle
D^{\mu\nu}(P)=\frac{(P^{2}-d)B^{\mu\nu}}{(P^{2}-b)(P^{2}-d)-a^{2}}+\frac{R^{\mu\nu}}{P^{2}-c}+\frac{(P^{2}-b)M^{\mu\nu}}{(P^{2}-b)(P^{2}-d)-a^{2}}+\frac{aN^{\mu\nu}}{(P^{2}-b)(P^{2}-d)-a^{2}}.$
(12)

Since, we are interested in the static inter-quark potential, we need only the
“00”-component of the full gluon propagator. The “00”-component is given by

$\displaystyle
D^{00}(P)=\frac{(P^{2}-d)\bar{u}^{2}}{(P^{2}-b)(P^{2}-d)-a^{2}},$ (13)

because $R^{00}=M^{00}=N^{00}=0$. The $Q\bar{Q}$ potential is given by the
static limit of the full gluon propagator. Since form factor $a(p_{0},p)$
vanishes in the static limit (see the appendix A). The resummed propagator
(13) becomes

$\displaystyle D^{00}(p_{0}=0,p)=-\frac{1}{(p^{2}+b(p_{0}=0,p))}.$ (14)

Thus, we are left only with the form factor $b(p_{0},p)$ which is needed to be
calculated in the static limit. So the next subsection is devoted to the
calculation of $b(p_{0},p)$.

###  2.1 Real and imaginary parts of the form factor $b(p_{0},p)$

In order to find the complex full gluon propagator, we will work out the real
and imaginary parts of the form factor $b(p_{0},p)$. Using the contraction
property of Eq. (8), $b(p_{0},p)$ is given by

$\displaystyle b(P)$ $\displaystyle=$ $\displaystyle
B^{\mu\nu}(P)\Pi_{\mu\nu}(P),$ $\displaystyle b(P)$ $\displaystyle=$
$\displaystyle\frac{{\bar{u}}^{\mu}{\bar{u}}^{\nu}}{{\bar{u}}^{2}}\Pi_{\mu\nu}(P),$
(15) $\displaystyle=$
$\displaystyle\left[\frac{u^{\mu}u^{\nu}}{{\bar{u}}^{2}}-\frac{(P\cdot
u)u^{\nu}P^{\mu}}{\bar{u}^{2}P^{2}}-\frac{(P\cdot
u)u^{\mu}P^{\nu}}{\bar{u}^{2}P^{2}}+\frac{(P\cdot
u)^{2}P^{\nu}P^{\mu}}{\bar{u}^{2}P^{4}}\right]\Pi_{\mu\nu}(P),$
$\displaystyle=$
$\displaystyle\frac{u^{\mu}u^{\nu}}{{\bar{u}}^{2}}\Pi_{\mu\nu}(P),$

where we have exploited transversality condition
$P^{\mu}\Pi_{\mu\nu}(P)=P^{\nu}\Pi_{\mu\nu}(P)=0$.

We will now evaluate the gluon self-energy in a strong magnetic field at
finite temperature and density. As we know only the quark loop gets influenced
in SMF while gluon loop will give only thermal contribution. The quark loop
contribution to the gluon self-energy is given as

$\displaystyle\Pi^{ab}_{\mu\nu}(P)$ $\displaystyle=$ $\displaystyle
i\int\frac{d^{4}K}{(2\pi)^{4}}{\rm
Tr}\left[gt_{b}\gamma_{\mu}S(K)gt_{a}\gamma_{\nu}S(Q)\right],$ (16)
$\displaystyle=$
$\displaystyle\sum_{f}\frac{ig^{2}\delta_{ab}}{2}\int\frac{d^{4}K}{(2\pi)^{4}}{\rm
Tr}\left[\gamma_{\mu}S(K)\gamma_{\nu}S(Q)\right],$

where $Q=(K-P)$ and ${\rm Tr}(t_{a}t_{b})=\frac{\delta_{ab}}{2}$ . The quark
propagator $S(K)$ in SMF limit reads as [50, 51]

$\displaystyle
iS(K)=ie^{-\frac{K^{2}_{\perp}}{|q_{f}B|}}\frac{(\not{K}_{\parallel}+m_{f})}{(K^{2}_{\parallel}-{m_{f}}^{2})}(1-i\gamma_{1}\gamma_{2}),$
(17)

where $q_{f}$ and $m_{f}$ refers to the charge and mass of the $f^{th}$ quark
flavor, respectively. We choose the metric tensor as

$\displaystyle g^{\mu\nu}_{\parallel}$ $\displaystyle=$ $\displaystyle{\rm
diag}(1,0,0-1),$ $\displaystyle~{}g^{\mu\nu}_{\perp}$ $\displaystyle=$
$\displaystyle{\rm diag}(0,-1,-1,0),$

and the four-momentum can be decomposed as

$\displaystyle K^{\mu}_{\parallel}$ $\displaystyle=$
$\displaystyle(k_{0},0,0,k_{z}),$ (18) $\displaystyle K^{\mu}_{\perp}$
$\displaystyle=$ $\displaystyle(0,k_{x},k_{y},0),$ (19) $\displaystyle
K_{\parallel}^{2}$ $\displaystyle=$ $\displaystyle k_{0}^{2}-k_{z}^{2},$ (20)
$\displaystyle K_{\perp}^{2}$ $\displaystyle=$ $\displaystyle
k_{x}^{2}+k_{y}^{2},$ (21)

In the SMF, the momentum integration can be decomposed into longitudinal
($\parallel$ ) and transverse ($\perp$) components with respect to the
magnetic field, so the gluon self-energy (16) can be factorized into
($\parallel$ ) and ($\perp$) components as

$\displaystyle\Pi_{\mu\nu}(P)=X(K_{\perp})i\int\frac{d^{2}K_{\parallel}}{(2\pi)^{2}}{\rm
Tr}\left[\gamma_{\mu}\frac{(\not{K}_{\parallel}+m_{f})}{(K^{2}_{\parallel}-{m_{f}}^{2})}(1-i\gamma_{1}\gamma_{2})\gamma_{\nu}\frac{(\not{Q}_{\parallel}+m_{f})}{(Q^{2}_{\parallel}-{m_{f}}^{2})}(1-i\gamma_{1}\gamma_{2})\right],$
(22)

where the transverse part is given by

$\displaystyle X(K_{\perp})$ $\displaystyle=$
$\displaystyle\sum_{f}\frac{g^{2}}{2}\int\frac{d^{2}K_{\perp}}{(2\pi)^{2}}~{}e^{\frac{-K^{2}_{\perp}-Q^{2}_{\perp}}{|q_{f}B|}},$
(23) $\displaystyle=$
$\displaystyle\sum_{f}e^{-\frac{P^{2}_{\perp}}{2|q_{f}B|}}\frac{g^{2}|q_{f}B|}{2\pi}.$

In the strong magnetic field $(|q_{f}B|>>P^{2}_{\perp})$, we can approximate
$e^{-\frac{P^{2}_{\perp}}{2|q_{f}B|}}\approx 1$. After substituting transverse
part from Eq. (23), the Eq. (22) becomes,

$\displaystyle\Pi_{\mu\nu}(P)$ $\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi}T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\frac{L_{\mu\nu}}{(K^{2}_{\parallel}-{m_{f}}^{2})(Q^{2}_{\parallel}-{m_{f}}^{2})},$
(24)

where

$\displaystyle L_{\mu\nu}$ $\displaystyle=$ $\displaystyle
u_{\mu}u_{\nu}(k_{0}q_{0}+k_{3}q_{3}+m_{f}^{2})+n_{\mu}n_{\nu}(k_{0}q_{0}+k_{3}q_{3}-m_{f}^{2})$
(25) $\displaystyle+(u_{\mu}n_{\nu}+n_{\mu}u_{\nu})(k_{0}q_{3}+k_{3}q_{0}).$

Here the strong coupling $g$ is the function of temperature, chemical
potential and magnetic field. It is given by [52]

$\displaystyle\alpha_{s}(\Lambda^{2},eB)$ $\displaystyle=$
$\displaystyle\frac{g^{2}}{4\pi}$ (26) $\displaystyle=$
$\displaystyle\frac{\alpha_{s}(\Lambda^{2})}{1+b_{1}\alpha_{s}(\Lambda^{2})\ln\left(\frac{\Lambda^{2}}{\Lambda^{2}+eB}\right)},$

with

$\displaystyle\alpha_{s}(\Lambda^{2})=\frac{1}{b_{1}\ln\left(\frac{\Lambda^{2}}{\Lambda_{\overline{MS}}^{2}}\right)},$
(27)

where $\Lambda$ is set at $2\pi\sqrt{T^{2}+\frac{\mu^{2}}{\pi^{2}}}$ for
quarks and $2\pi T$ for gluons, $b_{1}=\frac{11N_{c}-2N_{f}}{12\pi}$ and
$\Lambda_{\overline{MS}}=0.176GeV$.

Now substituting $\Pi_{\mu\nu}(P)$ from Eq. (24) in Eq. (15) we get

$\displaystyle b(P)$ $\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi\bar{u}^{2}}T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\frac{(k_{0}q_{0}+k_{3}q_{3}+m_{f}^{2})}{(K^{2}_{\parallel}-{m_{f}}^{2})(Q^{2}_{\parallel}-{m_{f}}^{2})}.$
(28)

In the static limit, the real and imaginary parts of $b(p_{0},p)$ from the
quark loop are obtained as (see the appendix B)

$\displaystyle{\rm Re}~{}b(p_{0}=0,p)$ $\displaystyle=$
$\displaystyle\sum_{f}\frac{g^{2}|q_{f}B|}{4\pi^{2}T}\int_{0}^{\infty}{dk_{3}}~{}\bigg{\\{}n^{+}(E_{1})(1-n^{+}(E_{1}))$
(29) $\displaystyle\quad+n^{-}(E_{1})(1-n^{-}(E_{1}))\bigg{\\}}$
$\displaystyle\left[\frac{{\rm Im}~{}b(p_{0},p)}{p_{0}}\right]_{p_{0}=0}$
$\displaystyle=$ $\displaystyle\sum_{f}g^{2}\frac{|q_{f}B|m_{f}^{2}}{16\pi
T(\frac{p_{3}^{2}}{4}+m_{f}^{2})}\bigg{\\{}n^{+}(\Omega)(1-n^{+}(\Omega))$
(30) $\displaystyle\quad+n^{-}(\Omega)(1-n^{-}(\Omega))\bigg{\\}},$

where $E_{1}=\sqrt{k_{3}^{2}+m_{f}^{2}}$ and
$\Omega=\sqrt{\frac{p_{3}^{2}}{4}+m_{f}^{2}}$. The distribution functions
$n^{+}(E_{1})$ and $n^{-}(E_{1})$ for quarks and anti-quarks, respectively are
given as

$\displaystyle n^{\pm}(E_{1})=\frac{1}{e^{\beta(E_{1}\mp\mu)}+1}.$ (31)

Eq. (30) can be further simplified using the identity

$\displaystyle
n^{\pm}(E)(1-n^{\pm}(E))=\frac{1}{2\left(1+\cosh(\beta(E\mp\mu)\right)},$ (32)

as (neglecting $O(\frac{\mu^{2}}{T^{2}})$ and higher order terms since we are
working in the limit $\mu<T$ )

$\displaystyle\left[\frac{{\rm Im}~{}b(p_{0},p)}{p_{0}}\right]_{p_{0}=0}$
$\displaystyle=$ $\displaystyle\sum_{f}g^{2}\frac{|q_{f}B|m_{f}^{2}}{8\pi
T}~{}\frac{1}{p_{3}^{2}}.$ (33)

The “00” component of gluon self-energy tensor as a consequence of the gluon-
loop is given by [53, 54]

$\displaystyle\Pi_{00}(p_{0},p)=-g^{2}T^{2}\frac{N_{c}}{3}\left(\frac{p_{0}}{2p}\ln\frac{p_{0}+p+i\epsilon}{p_{0}-p+i\epsilon}-1\right),$
(34)

we extract the real and imaginary parts of Eq. (34) which are found to be

$\displaystyle{\rm Re}~{}b_{0}(p_{0}=0)$ $\displaystyle=$ $\displaystyle
g^{2}T^{2}\left(\frac{N_{c}}{3}\right),$ (35) $\displaystyle\left[\frac{{\rm
Im}~{}b_{0}(p_{0},p)}{p_{0}}\right]_{p_{0}=0}$ $\displaystyle=$ $\displaystyle
g^{2}T^{2}\left(\frac{N_{c}}{3}\right)\frac{\pi}{2p}.$ (36)

The square of the Debye mass in the strong magnetic field at finite
temperature and chemical potential is given by

$\displaystyle m_{D}^{2}$ $\displaystyle=$
$\displaystyle(b+b_{0})\bar{u}^{2}|_{p_{0}=0},$ (37) $\displaystyle=$
$\displaystyle m_{\rm q,D}^{2}+m_{\rm g,D}^{2},$ (38)

where the quark-loop contribution ($m_{\rm q,D}^{2}$) to the Debye mass is

$\displaystyle m_{\rm
q,D}^{2}(T,\mu;B)=\sum_{f}g^{2}\frac{|q_{f}B|}{4\pi^{2}T}\int_{0}^{\infty}{dk_{3}}~{}\\{n^{+}(E_{1})(1-n^{+}(E_{1}))+n^{-}(E_{1})(1-n^{-}(E_{1}))\\}.$
(39)

It is worthwhile to mention here that the dependence of chemical potential
($\mu$) in $m_{\rm q,D}$ is only manifested in the finite (physical) quark
masses otherwise it simply reduces to the known result [43]

$\displaystyle m_{\rm q,D}^{2}(B)=\sum_{f}\frac{g^{2}~{}|q_{f}B|}{4\pi^{2}}.$
(40)

On the other hand, the gluon contribution is as usual given by

$\displaystyle m_{\rm g,D}^{2}(T)=\frac{N_{C}}{3}g^{2}T^{2}$ (41)

To visualize the effect of finite baryon asymmetry on the collective modes of
a strongly magnetized hot QCD medium, we have plotted the Debye mass ($m_{D}$)
as a function of temperature (in units of $T_{c}$) with the increasing quark
chemical potentials at a fixed magnetic field strength ($eB=15~{}m_{\pi}^{2}$)
(left panel of Fig. 1). We have seen that $m_{D}$ increases with $T$, as
expected but on the contrary it decreases with $\mu$, which is more
significant in the low temperature region [55]. This finding can be better
understood if we plot the same with respect to $\mu$ at a fixed temperature,
$T=200$ MeV and magnetic field strength, $eB=15~{}m_{\pi}^{2}$, wherein Debye
mass decreases with $\mu$ and this trend is pronounced at large $\mu$ under
consideration. The Debye mass gets reduced in the presence of the strong $B$
in comparison to the $B=0$ case.

|
---|---
(a) | (b)
Figure 1: Variation of the Debye mass with the a) temperature, b) quark
chemical potential

##  3 Medium modified $Q\bar{Q}$ potential

In this section, we will explore the effect of surplus of baryons over anti-
baryons to the medium modification of $Q\bar{Q}$ potential immersed in a
deconfined medium with an ambient strong magnetic field. The inverse Fourier
transform of the resummed gluon propagator in the static limit serves the
desired medium-modification in the coordinate space as [56]

$\displaystyle
V(r;T,B,\mu)=C_{F}~{}g^{2}\int\frac{d^{3}p}{(2\pi)^{3}}(e^{ip.r}-1)~{}D^{00}(p_{0}=0,p),$
(42)

where $C_{F}(=4/3)$ is the Casimir factor and $D^{00}(p_{0}=0,p)$ is the
static limit of the complex full gluon propagator, whose real and imaginary
parts are needed to obtain the complex inter-quark potential. The
$r$-independent term (which is the perturbative free energy of quarkonium at
infinite separation) has been subtracted to renormalize the heavy quark free
energy.

###  3.1 The real and imaginary parts of the resummed gluon propagator

The static limit of the real-part of “00”-component of full gluon propagator
is given by using Eq. (14) and Eq. (37) as

$\displaystyle{\rm Re}~{}D^{00}(p_{0}=0)=-\frac{1}{p^{2}+m_{D}^{2}}.$ (43)

Similarly, the imaginary-part reads [57]

$\displaystyle{\rm Im}~{}D^{00}(p_{0},p)=\frac{2T}{p_{0}}\frac{{\rm
Im}~{}b(p_{0},p)}{(P^{2}-{\rm Re}~{}b(p_{0},p))^{2}+({\rm
Im}~{}b(p_{0},p))^{2}},$ (44)

which can be further simplified as

$\displaystyle{\rm Im}~{}D^{00}(p_{0},p)=2T\frac{\left[\frac{{\rm
Im}~{}b(p_{0},p)}{p_{0}}\right]}{(P^{2}-{\rm
Re}~{}b(p_{0},p))^{2}+\left(p_{0}\left[\frac{{\rm
Im}~{}b(p_{0},p)}{p_{0}}\right]\right)^{2}},$ (45)

the above Eq. (45) in static limit ($p_{0}=0$) reduces to

$\displaystyle{\rm Im}~{}D^{00}(p_{0}=0)=2T\frac{\left[\frac{{\rm
Im}~{}b(p_{0},p)}{p_{0}}\right]_{p_{0}=0}}{(p^{2}+m_{D}^{2})^{2}},$ (46)

where we have exploited ${\rm Re}~{}b(p_{0}=0,p)=m_{D}^{2}$. Using Eq. (33),
the imaginary part of $D^{00}(p_{0}=0,p)$ is given as

$\displaystyle{\rm
Im}~{}D^{00}(p_{0}=0,p)=\sum_{f}\frac{g^{2}|q_{f}B|m_{f}^{2}}{4\pi}\frac{1}{p_{3}^{2}(p^{2}+m_{D}^{2})^{2}}.$
(47)

Now, we will revisit the procedure to handle the large distance behaviour of
the $Q\bar{Q}$ potential. A phenomenological model has been proposed to study
the string part of the potential in [38] where the authors have added a
phenomenological nonperturbative term to the HTL full gluon propagator in
order to include the effects due to the low frequency modes incorporated in
the dimension two gluon condensates. The real and imaginary parts of the
phenomenological nonperturbative (NP) term are given by

$\displaystyle{\rm Re}~{}D^{00}_{\rm
NP}(p_{0}=0,p)=-\frac{m_{G}^{2}}{(p^{2}+m_{D}^{2})^{2}},$ (48)
$\displaystyle{\rm Im}~{}D^{00}_{\rm NP}(p_{0}=0,p)=\frac{2\pi
Tm^{2}_{g}m_{G}^{2}}{p(p^{2}+m_{D}^{2})^{3}},$ (49)

where the dimension two constant $m_{G}^{2}$ can be expressed in terms of the
string tension as $\sigma=\alpha m_{G}^{2}/2$. Eqs. (48) and (49) will induce
the string contribution in the complex $Q\bar{Q}$ potential. Finally, the real
and imaginary parts of the “00”-component of the full gluon propagator can be
written as

$\displaystyle{\rm Re}~{}D^{00}(p_{0}=0,p)$ $\displaystyle=$
$\displaystyle-\frac{1}{p^{2}+m_{D}^{2}}-\frac{m_{G}^{2}}{(p^{2}+m_{D}^{2})^{2}},$
(50) $\displaystyle{\rm Im}~{}D^{00}(p_{0}=0,p)$ $\displaystyle=$
$\displaystyle\sum_{f}\frac{g^{2}|q_{f}B|m_{f}^{2}}{4\pi}\frac{1}{p_{3}^{2}(p^{2}+m_{D}^{2})^{2}}+\frac{\pi
Tm_{g}^{2}}{p(p^{2}+m_{D}^{2})^{2}}+\frac{2\pi
Tm_{g}^{2}m_{G}^{2}}{p(p^{2}+m_{D}^{2})^{3}}.$ (51)

We will use Eqs. (50) and (51) to derive the (complex) $Q\bar{Q}$ potential in
the next subsection.

###  3.2 Real and Imaginary parts of the $Q\bar{Q}$ potential

In this subsection, we will compute the real- and imaginary-parts of inter
quark potential between $Q\bar{Q}$ in a strongly magnetized hot quark matter
with finite chemical potential by substituting the real- and imaginary-parts
of the full gluon propagator, respectively, into the definition (42). Thus,
the real-part of $Q\bar{Q}$ potential (with $\hat{r}=rm_{D}$) is obtained as

$\displaystyle{\rm{Re}}~{}V(r;T,B,\mu)$ $\displaystyle=$
$\displaystyle-\frac{4}{3}\alpha_{s}\left(\frac{e^{-\hat{r}}}{r}+m_{D}(T,\mu,B)\right)+\frac{4}{3}\frac{\sigma}{m_{D}(T,\mu,B)}\left(1-e^{-\hat{r}}\right),$
(52)

where the dependences of temperature, chemical potential and magnetic field in
Debye mass get translated into the medium modified inter-quark potential . The
string term in Eq. (52) comes from the nonperturbative part in Eq. (50).
Mainly we wish to visualize the modification due to the sole effect of baryon
asymmetry on the real-part of $Q\bar{Q}$ potential as a function of inter-
quark distance ($r$) for increasing $\mu$’s (in Fig. 2). While plotting the
real-part, we have abandoned the $r$-independent terms, which are needed in
the potential (52) to obtain its form in $T\rightarrow 0$ limit. In the Fig. 2
(a), we have displayed ${\rm Re}~{}V(r;T,B,\mu)$ for $\mu=0$, $60$ and $100$
MeV at fixed temperature $T=200$ MeV and strong magnetic field
$eB=15~{}m_{\pi}^{2}$. We have observed that the real-part becomes more
attractive at finite chemical potential in comparison to $\mu=0$. This strong
nature of the $Q\bar{Q}$ potential can be attributed to the less screening in
the presence of baryon asymmetry ($\mu\neq 0$) in the strongly magnetized QCD
medium. We have displayed the ${\rm Re}~{}V(r;T,B,\mu)$ considering the same
values of $\mu$ at $T=250$ MeV (in the Fig. 2 (b)) and have found the same
behavior. It is evident from the Fig 2 (b) that as the temperature rises the
effect of $\mu$ diminishes.

We will now evaluate the imaginary-part of the inter-quark potential using the
imaginary-part of full gluon propagator from Eq. (51) into Eq. (42), which is
separable into perturbative and nonperturbative (NP) parts as

$\displaystyle{\rm Im}~{}V(r;T,B,\mu)$ $\displaystyle=$ $\displaystyle{\rm
Im}~{}V_{\rm perturbative}(r;T,B,\mu)+{\rm Im}~{}V_{\rm NP}(r;T,B,\mu).$ (53)

The perturabtive part is separated into quark-loop (q) and gluon-loop (g)
contributions as

$\displaystyle{\rm Im}~{}V_{\rm perturbative}(r;T,B,\mu)={\rm
Im}~{}V_{q}(r;T,B,\mu)+{\rm Im}~{}V_{g}(r;T,B,\mu),$ (54)

where the quark-loop contribution has been calculated as

$\displaystyle{\rm Im}~{}V_{q}(r,T,B,\mu)$ $\displaystyle=$
$\displaystyle\sum_{f}\alpha_{s}g^{2}m_{f}\frac{|q_{f}B|}{3\pi^{2}}\bigg{[}\frac{\pi}{2m_{D}^{3}}-\frac{\pi
e^{-\hat{r}}}{2m_{D}^{3}}-\frac{\pi\hat{r}e^{-\hat{r}}}{2m_{D}^{3}}$ (55)
$\displaystyle-\frac{2\hat{r}}{m_{D}}\int_{0}^{\infty}\frac{pdp}{(p^{2}+m_{D}^{2})^{2}}~{}\int_{0}^{pr}\frac{\sin{t}}{t}dt\bigg{]}.$

and the gluon-loop contribution is

$\displaystyle{\rm Im}~{}V_{g}(r;T,B,\mu)$ $\displaystyle=$
$\displaystyle-\frac{4}{3}\frac{\alpha_{s}Tm_{g}^{2}}{m_{D}^{2}}\psi_{1}(\hat{r}),$
(56)

where the function $\psi_{1}(\hat{r})$ is given by [38, 46]

$\displaystyle\psi_{1}(\hat{r})$ $\displaystyle=$ $\displaystyle
2\int_{0}^{\infty}\frac{zdz}{(z^{2}+1)^{2}}\left[1-\frac{\sin
z\hat{r}}{z\hat{r}}\right],$ (57)

which can be further simplified in the small $\hat{r}$ limit as

$\displaystyle\psi_{1}(\hat{r})$ $\displaystyle\approx$
$\displaystyle-\frac{1}{9}{\hat{r}}^{2}\left(3\ln\hat{r}-4+3\gamma_{E}\right).$
(58)

Similarly, we calculate the imaginary-part of the string part of the
$Q\bar{Q}$ potential using the nonperturbative term in the full gluon
propagator from Eq. (51) in Eq. (42), we get

$\displaystyle{\rm Im}~{}V_{\rm NP}(r;T,B,\mu)$ $\displaystyle=$
$\displaystyle-\frac{16\sigma Tm_{g}^{2}}{3m_{D}^{4}}\psi_{2}(\hat{r}),$ (59)

where the function $\psi_{2}(\hat{r})$ is given in [38, 46]

$\displaystyle\psi_{2}(\hat{r})$ $\displaystyle=$ $\displaystyle
2\int_{0}^{\infty}\frac{zdz}{(z^{2}+1)^{3}}\left[1-\frac{\sin
z\hat{r}}{z\hat{r}}\right],$ (60)

which further takes the form in limit $(\hat{r}\ll 1)$

$\displaystyle\psi_{2}(\hat{r})$ $\displaystyle\approx$
$\displaystyle\frac{{\hat{r}}^{2}}{12}+\frac{{\hat{r}}^{4}}{900}\left(15\ln\hat{r}-23+15\gamma_{E}\right).$
(61)

While the real-part will explore the effect of baryon asymmetry on the binding
energy, the effect on the dissociation will be understood through the
imaginary-part (in Fig. 3). The magnitude of imaginary-part is found to
decrease in baryon asymmetric matter, $\mu\neq 0$, compared to its counterpart
at $\mu=0$ and it decreases further as $\mu$ rises. We have conducted a
similar investigation at $T=250$ MeV in the Fig. 3 (b) and have found the same
behavior, however the effect of quark chemical potential on imaginary-part is
less pronounced at high temperatures because the effect of the $\mu$ on the
Debye mass is not much visible.

Now we will compare our results of the $Q\bar{Q}$ potential in the presence of
the strong magnetic field ($eB=15m_{\pi}^{2}$) with those in the absence of
magnetic field ($B=0$). The real and imaginary parts of the resummed gluon
propagator (in the static limit) are given in the absence of magnetic field as

$\displaystyle{\rm Re}~{}D^{00}(p_{0}=0,p)$ $\displaystyle=$
$\displaystyle-\frac{1}{p^{2}+{m^{\prime}}_{D}^{2}}-\frac{m_{G}^{2}}{(p^{2}+{m^{\prime}}_{D}^{2})^{2}},$
(62) $\displaystyle{\rm Im}~{}D^{00}(p_{0}=0,p)$ $\displaystyle=$
$\displaystyle\frac{\pi
T{m^{\prime}}_{D}^{2}}{p(p^{2}+{m^{\prime}}_{D}^{2})^{2}}+\frac{2\pi
T{m^{\prime}}_{D}^{2}m_{G}^{2}}{p(p^{2}+{m^{\prime}}_{D}^{2})^{3}}.$ (63)

respectively. Using Eq. (42), the real and imaginary parts of $Q\bar{Q}$
potential are found to be

$\displaystyle{\rm{Re}}~{}V(r;T,\mu)$ $\displaystyle=$
$\displaystyle-\frac{4}{3}\alpha^{\prime}_{s}\left(\frac{e^{-\hat{r}}}{r}+{m^{\prime}}_{D}\right)+\frac{4}{3}\frac{\sigma}{{m^{\prime}}_{D}}\left(1-e^{-\hat{r}}\right),$
(64) $\displaystyle{\rm{Im}}~{}V(r;T,\mu)$ $\displaystyle=$
$\displaystyle-\frac{4}{3}\alpha^{\prime}_{s}T\psi_{1}(\hat{r})-\frac{16\sigma
T}{3{m^{\prime}}_{D}^{2}}\psi_{2}(\hat{r}),$ (65)

respectively. Here $\alpha^{\prime}$ is the strong coupling constant, which is
given by

$\displaystyle\alpha^{\prime}_{s}(T)=\frac{6\pi}{(33-2N_{f})\ln\left(\frac{Q}{\Lambda_{QCD}}\right)},$
(66)

where $Q=2\pi\sqrt{T^{2}+\frac{\mu^{2}}{\pi^{2}}}$ and ${m^{\prime}}_{D}$ is
the Debye mass, which reads as

$\displaystyle{m^{\prime}}_{D}^{2}=g^{\prime
2}T^{2}\left\\{\frac{N_{c}}{3}+\frac{N_{f}}{6}\left(1+\frac{3\mu^{2}}{\pi^{2}T^{2}}\right)\right\\}.$
(67)

We notice from Fig 1 (a) that magnitude of the Debye mass in the absence of
$B$ is greater in comparison to that in strong $B$ environment. So the real
part of the potential gets screened at higher amount and becomes less
attractive in comparison to $B\neq 0$ case. [seen in Fig. 2 (a)]. Similar
observation we notice when the temperature of the medium is 250 MeV [seen in
Fig. 2 (b)]. On the other hand, the magnitude of the imaginary part gets
enhanced in comparison to $B\neq 0$ case [seen in Fig. 3 (a) and (b)]. The
impact of $\mu$ on the Debye mass (and hence on the $Q\bar{Q}$ potential) in
the absence of magnetic field is not much visible for $\mu=60$ and $100$ MeV.
The Debye mass gets slightly increased as the strength of $\mu$ is raised.
Consequently, the real part becomes less attractive and the magnitude of the
imaginary part increases.

|
---|---
(a) | (b)
Figure 2:  Variation of ${\rm Re}~{}V(r,T,B,\mu)$ with inter-quark separation
($r$) at different strengths of the quark chemical potential ($\mu$).

|
---|---
(a) | (b)
Figure 3:  Variation of ${\rm Im}~{}V(r,T,B,\mu)$ with inter-quark separation
($r$) at different strengths of the quark chemical potential ($\mu$).

##  4 Properties of quarkonia

We will now explore how the presence of the baryon asymmetry in the strongly
magnetized hot QCD medium affects the properties of heavy quarkonia. We will
compute the binding energy and decay width of the $Q\bar{Q}$ system with the
help of the real and imaginary parts of the $Q\bar{Q}$ potential,
respectively.

###  4.1 Binding energy (BE)

We have solved the radial part of the Schrödinger equation numerically
exploiting the real-part of the potential to obtain the energy eigenvalues
which are utilized to calculate the binding energy of quarkonia. We have
examined the effect of quark chemical potential on the binding energy of the
heavy quarkonium states in Fig. 4. For that purpose, we have computed the BE
of $J/\psi$ and $\Upsilon$ at $\mu=0,60$ and $100$ MeV while fixing
$eB=15m_{\pi}^{2}$. We have observed that BE decreases with $T$, which is
justified since the screening mass increases with $T$. The magnitude of the
binding energy is slightly higher for $\mu\neq 0$ in comparison to $\mu=0$
case. This behavior can be understood in terms of the softening of the the
screening mass in the presence of the baryon asymmetry in the medium which
leads to the stronger nature of the real-part of inter-quark potential hence
slightly enhanced values of the binding energy.

|
---|---
(a) | (b)
Figure 4:  Variation of binding energy of $J/\psi$ (a) and $\Upsilon$ (b) with
$T$ at different strengths of the quark chemical potential ($\mu$).

###  4.2 Thermal width

We will now explore the broadening of the thermal width of the $Q\bar{Q}$
bound states in a strongly magnetized hot QCD medium through the imaginary
part of the $Q\bar{Q}$ potential. In small distance limit, the imaginary part
of the potential can be treated as the perturbation to the vaccum potential
which gives the thermal width ($\Gamma$) for a particular resonance state as

$\displaystyle\Gamma({T,B,\mu})=-2\int_{0}^{\infty}{\rm
Im}~{}V(r;T,B,\mu)|\Psi(r)|^{2}d\tau,$ (68)

where we choose $\Psi(r)$ as the Coloumbic wave function which reads

$\displaystyle\Psi(r)=\frac{1}{\sqrt{\pi a_{0}^{3}}}e^{-r/a_{0}}.$ (69)

Here $a_{0}$ refers to the Bohr radius of the $Q\bar{Q}$ bound state.

|
---|---
(a) | (b)
Figure 5: Variation of Decay width of $J/\psi$ (a) and $\Upsilon$ (b) with $T$
(in units of $T_{c}$) at different strengths of $\mu$.

In order to decipher the effect of the finite chemical potential on the
thermal width ($\Gamma$) of the $Q\bar{Q}$ bound states, we have evaluated
$\Gamma$ numerically for $J/\psi$ and $\Upsilon$ with respect to $T$ (in Fig.
5) for $\mu=0,60$ and $100$ MeV. We have observed that $\Gamma$ increases with
the temperature while it gets decreased in the presence of baryon asymmetry
($\mu\neq 0$) in the medium. This behavior can be explained in terms of ${\rm
Im}~{}V(r;T,B,\mu)$ whose magnitude gets decreased in the presence of $\mu$.

###  4.3 Dissociation of quarkonia

In this section we will study the dissociation process of the $Q\bar{Q}$ bound
states in a baryon asymmetric strongly magnetized thermal QCD medium and will
see how the dissociation temperature ($T_{D}$) of quarkonia is affected by the
presence of finite amount of $\mu$. We have used the criterion ($\Gamma$):
$\Gamma\geq 2$ binding energy [58] to evaluate the values of the dissociation
points for $J/\psi$ and $\Upsilon$.

|
---|---
(a) | (b)
Figure 6:  Competition between the $\Gamma$ and 2× BE for $J/\psi$ (a) and $\Upsilon$ (b) with respect to $T$ at different strengths of the quark chemical potential ($\mu$)  |  $T_{D}$ (in terms of $T_{c}$), $eB=15~{}m_{\pi}^{2}$
---|---
State | $J/\psi$ | $\Upsilon$
$\mu$ =0 | 1.64 | 1.95
$\mu=60$ | 1.69 | 1.97
$\mu=100$ | 1.75 | 2.00
Table 1: $T_{D}$’s of $J/\psi$ and $\Upsilon$ at different strengths of
chemical potential ($\mu$).

We have plotted thermal width and twice of the binding energy with the
temperature in Fig. 6 and have found that $T_{D}$’s of $J/\psi$ and $\Upsilon$
increase slightly in baryon asymmetric QCD medium in comparison to baryonless
$(\mu=0)$ medium. The dissociation temperatures for $J/\psi$ are found to be
$1.64~{}T_{c}$ ,$1.69~{}T_{c}$,and $1.75~{}T_{c}$ at the $\mu=0,60$ and $100$
MeV respectively whereas $\Upsilon$ is dissociated at $1.95~{}T_{c}$,
$1.97~{}T_{c}$ and $2.00~{}T_{c}$ for $\mu=0,60$ and $100$ MeV respectively.

##  5 Conclusions

To conclude, we have examined the effects of quark chemical potential on the
properties of the quarkonia in baryon asymmetric strongly magnetized hot QCD
medium. First, we have given a revisit to the general covariant tensor
structure of gluon self-energy in above mentioned medium and computed the
relevant form factors. We use these form factors in the calculation of the
(complex) full gluon propagator which further get translated into the
(complex) $Q\bar{Q}$ potential. We have added a phenomenological non-
perturbative term induced by the dimension two gluon condensate to the usual
HTL resummed propagator to evaluate the medium modification to string part of
the $Q\bar{Q}$ potential. The real-part becomes more attractive while
magnitude of the imaginary-part gets decreased in the baryon asymmetric
medium. We evaluate the binding energy of the $Q\bar{Q}$ bound states solving
the Schrödinger equation numerically considering the real-part of the
potential whereas the imaginary-part gives the thermal width. The binding
energy of $J/\psi$ and $\Upsilon$ get enhanced while decay width gets
decreased at finite $\mu$ in comparison to baryonless medium $(\mu=0)$. This
increment in the binding energy is attributed to the stronger nature of the
$Q\bar{Q}$ potential in presence of $\mu$. we have finally explored the
dissociation process of heavy quark bound states in the above mentioned medium
and observed that $T_{D}$’s for $J/\psi$ and $\Upsilon$ now attains slightly
higher values in the baryon asymmetric medium. $J/\psi$ is dissociated at
$1.64~{}T_{c}$, $1.69~{}T_{c}$ and $1.75~{}T_{c}$ for $\mu=0,60,100$ MeV,
respectively whereas $\Upsilon$ is dissociated at $1.95~{}T_{c}$,
$1.97~{}T_{c}$ and $2.00~{}T_{c}$ for the same strengths of $\mu$. This study
leads to the conclusion that baryon asymmetry in strongly magnetized hot QCD
medium prevents slightly early dissociation of quarkonia in comparison to
baryonless medium.

## Acknowledgements

One of us BKP acknowledges the financial assistance from the CSIR (Grant No.03
(1407)/17/EMR-II), Government of India. We present explicit calculation of the
form factor $a(p_{0},p)$ and $b(p_{0},p)$ in the following appendices.

##  Appendix A Calculation of the form factor $a(P)$

We will calculate the form factor $a(p_{0},p)$ using the imaginary time
formalism of the finite temperature field theory. We can write from Eq. (11)

$\displaystyle a(p_{0},p)$ $\displaystyle=$ $\displaystyle
i\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi}\frac{p_{0}p_{3}}{\sqrt{\bar{u}^{2}}\sqrt{\bar{n}^{2}}\tilde{P}^{2}}\int\frac{d^{2}K_{\parallel}}{(2\pi)^{2}}\frac{\left[k_{0}^{2}+k_{3}^{2}+m_{f}^{2}\right]}{(K^{2}_{\parallel}-m^{2}_{f})(Q^{2}_{\parallel}-m_{f}^{2})}$
(A.70) $\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi}\frac{p_{0}p_{3}}{\sqrt{\bar{u}^{2}}\sqrt{\bar{n}^{2}}\tilde{P}^{2}}T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\frac{\left[k_{0}^{2}+k_{3}^{2}+m_{f}^{2}\right]}{(K^{2}_{\parallel}-m^{2}_{f})(Q^{2}_{\parallel}-m_{f}^{2})}$
$\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi}\frac{p_{0}p_{3}}{\sqrt{\bar{u}^{2}}\sqrt{\bar{n}^{2}}\tilde{P}^{2}}T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\left(\frac{1}{(K^{2}_{\parallel}-m^{2}_{f})}+\frac{2(k^{2}_{3}+m^{2}_{f})}{(K^{2}_{\parallel}-m^{2}_{f})(Q^{2}_{\parallel}-m_{f}^{2})}\right)$
$\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi}\frac{p_{0}p_{3}}{\sqrt{\bar{u}^{2}}\sqrt{\bar{n}^{2}}\tilde{P}^{2}}(I_{1}(P)+I_{2}(P))$

where

$\displaystyle I_{1}(P)$ $\displaystyle=$ $\displaystyle
T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\frac{1}{(K^{2}_{\parallel}-m^{2}_{f})},$
(A.71) $\displaystyle=$
$\displaystyle-\int\frac{dk_{3}}{2\pi}~{}\frac{(1-n^{+}(E_{1})-n^{-}(E_{1}))}{2E_{1}},$

the first term gives the nonleading contribution in $T$, retaining only
leading order term we get

$\displaystyle
I_{1}(P)=\int\frac{dk_{3}}{2\pi}~{}\frac{(n^{+}(E_{1})+n^{-}(E_{1}))}{2E_{1}},$
(A.72)

where $n^{+}(E_{1})$ and $n^{-}(E_{1})$ are the distribution function for the
quarks and anti quarks, respectively and are given as

$\displaystyle n^{+}(E_{1})=\frac{1}{e^{\beta(E_{1}-\mu)}+1}$ (A.73)
$\displaystyle n^{-}(E_{1})=\frac{1}{e^{\beta(E_{1}+\mu)}+1}.$ (A.74)

Now taking the second term in (A.70)

$\displaystyle I_{2}(P)$ $\displaystyle=$ $\displaystyle
T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\frac{2(k^{2}_{3}+m^{2}_{f})}{(K^{2}_{\parallel}-m^{2}_{f})(Q^{2}_{\parallel}-m_{f}^{2})},$
(A.75) $\displaystyle=$
$\displaystyle-\int\frac{dk_{3}}{2\pi}~{}\frac{2(k^{2}_{3}+m^{2}_{f})}{4E_{1}E_{2}}\left(\frac{(1-n^{+}(E_{1})-n^{-}(E_{2}))}{(i\omega-
E_{1}-E_{2})}+\frac{(n^{-}(E_{1})-n^{-}(E_{2}))}{(i\omega+E_{1}-E_{2})}\right.$
$\displaystyle\left.+\frac{(n^{+}(E_{1})-n^{+}(E_{2}))}{(i\omega-
E_{1}+E_{2})}-\frac{(1-n^{-}(E_{1})-n^{+}(E_{2}))}{(i\omega+E_{1}+E_{2})}\right),$

where $E_{1}=\sqrt{k_{3}^{2}+m_{f}^{2}}$ and
$E_{2}=\sqrt{(k_{3}-p_{3})^{2}+m_{f}^{2}}$. In the HTL approximation the Eq.
(A.75) reduces to

$\displaystyle
I_{2}(P)=-\int\frac{dk_{3}}{2\pi}\left(\frac{n^{+}(E_{1})+n^{-}(E_{1}))}{2E_{1}}-\frac{dn^{-}(E_{1})}{dE_{1}}\frac{p_{3}}{2(p_{0}+p_{3})}+\frac{dn^{+}(E_{1})}{dE_{1}}\frac{p_{3}}{2(p_{0}-p_{3})}\right).$
(A.76)

Adding (A.72) and (A.76)

$\displaystyle
I_{1}(P)+I_{2}(P)=-\int\frac{dk_{3}}{2\pi}\left(\frac{dn^{+}(E_{1})}{dE_{1}}\frac{p_{3}}{2(p_{0}-p_{3})}-\frac{dn^{-}(E_{1})}{dE_{1}}\frac{p_{3}}{2(p_{0}+p_{3})}\right).$
(A.77)

Now putting ($I_{1}(P)+I_{2}(P)$) in (A.70) we get the form factor $a(P)$ as

$\displaystyle a(p_{0},p)$ $\displaystyle=$
$\displaystyle\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi}\frac{p_{0}p_{3}}{\sqrt{\bar{n}^{2}}\sqrt{{\bar{u}^{2}}}\tilde{P}^{2}}\int\frac{dk_{3}}{2\pi}\left(\frac{dn^{+}(E_{1})}{dE_{1}}~{}\frac{p_{3}}{2(p_{0}-p_{3})}-\frac{dn^{-}(E_{1})}{dE_{1}}~{}\frac{p_{3}}{2(p_{0}+p_{3})}\right)$
(A.78) $\displaystyle=$ $\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi
T}\frac{p_{0}p_{3}}{\sqrt{\bar{n}^{2}}\sqrt{{\bar{u}^{2}}}\tilde{P}^{2}}\int\frac{dk_{3}}{2\pi}\left(n^{+}(E_{1})(1-n^{+}(E_{1}))~{}\frac{p_{3}}{2(p_{0}-p_{3})}\right.$
$\displaystyle\left.-n^{-}(E_{1})(1-n^{-}(E_{1}))~{}\frac{p_{3}}{2(p_{0}+p_{3})}\right),$

which vanishes in the static limit ($p_{0}=0$).

##  Appendix B Calculation of the form factor $b(p_{0},p)$

We can write from Eq. (28) using HTL approximation

$\displaystyle b(p_{0},p)$ $\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi\bar{u}^{2}}T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\frac{\left[k_{0}^{2}+k_{3}^{2}+m_{f}^{2}\right]}{(K^{2}_{\parallel}-m^{2}_{f})(Q^{2}_{\parallel}-m_{f}^{2})},$
(B.79) $\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi\bar{u}^{2}}T\sum_{k_{0}}\int\frac{dk_{3}}{2\pi}\left(\frac{1}{(K^{2}_{\parallel}-m^{2}_{f})}+\frac{2(k^{2}_{3}+m^{2}_{f})}{(K^{2}_{\parallel}-m^{2}_{f})(Q^{2}_{\parallel}-m_{f}^{2})}\right),$
$\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi\bar{u}^{2}}(I_{1}(P)+I_{2}(P)),$

putting the value of $(I_{1}(P)+I_{2}(P))$ from eq. (A.77) in (B.79), the form
factor $b(P)$ can be written as

$\displaystyle
b(p_{0},p)=\sum_{f}\frac{g^{2}|q_{f}B|}{2\pi\bar{u}^{2}}\int\frac{dk_{3}}{2\pi}\left(\frac{dn^{+}(E_{1})}{dE_{1}}~{}\frac{p_{3}}{2(p_{0}-p_{3})}-\frac{dn^{-}(E_{1})}{dE_{1}}~{}\frac{p_{3}}{2(p_{0}+p_{3})}\right).$
(B.80)

The real part of the form factor $b(P)$ is given by

$\displaystyle{\rm Re}~{}b(p_{0},p)$ $\displaystyle=$
$\displaystyle-\sum_{f}\frac{g^{2}|q_{f}B|}{4\pi^{2}\bar{u}^{2}T}\int{dk_{3}}~{}\bigg{\\{}n^{+}(E_{1})(1-n^{+}(E_{1}))~{}\frac{p_{3}}{2(p_{0}-p_{3})}$
(B.81)
$\displaystyle-n^{-}(E_{1})(1-n^{-}(E_{1}))~{}\frac{p_{3}}{2(p_{0}+p_{3})}\bigg{\\}}.$

We have calculated the imaginary part of the form factor $b(p_{0},p)$ using
the identity

$\displaystyle{\rm Im}~{}b(p_{0},p)=\frac{1}{2i}\lim_{\eta\rightarrow
0}\left[b(p_{0}+i\epsilon,p)-b(p_{0}-i\epsilon,p)\right],$ (B.82)

along with the formula which gives the discontinuity across the real axis as

$\displaystyle\frac{1}{2i}\left(\frac{1}{p_{0}+\sum_{j}E_{j}+i\epsilon}-\frac{1}{p_{0}+\sum_{j}E_{j}-i\epsilon}\right)=-\pi\delta(p_{0}+\sum_{j}E_{j}).$
(B.83)

Thus using the above identities Eq. (B.82) and Eq. (B.83), the imaginary-part
of $b(p_{0},p)$ is found to be

$\displaystyle{\rm Im}~{}b(p_{0},p)$ $\displaystyle=$ $\displaystyle
g^{2}\frac{|q_{f}B|m_{f}^{2}p_{0}}{16\pi
T(\frac{p_{3}^{2}}{4}+m_{f}^{2})}\bigg{\\{}n^{+}(\Omega)(1-n^{+}(\Omega))+n^{-}(\Omega)(1-n^{-}(\Omega))\bigg{\\}},$
(B.84)

where $\Omega=\sqrt{\frac{p_{3}^{2}}{4}+m_{f}^{2}}$.

## References

  * [1] I. Arsene et al., BRAHMS Collaboration, Nucl. Phys. A 757,1 (2005).
  * [2] K. Adcox et al., PHENIX Collaboration, Nucl. Phys. A 757, 184 (2005).
  * [3] F. Carminati et al., ALICE Collaboration, J. Phys. G: Nucl. Part. Phys. 30, 1517 (2004).
  * [4] B. Alessandro et al., ALICE Collaboration, J. Phys. G: Nucl. Part. Phys. 32, 1295 (2006).
  * [5] P. Senger, Cent. Eur. J. Phys. 10, 1289 (2012)
  * [6] V. Skokov, A. Illarionov, V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009).
  * [7] V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski, S .A. Voloshin, Phys. Rev. C 83, 054911 (2011).
  * [8] K. Tuchin, Adv. High Energy Phys. 2013, 490495 (2013).
  * [9] L. McLerran, V. Skokov, Nucl. Phys. A 929, 184 (2014).
  * [10] S. Rath and B. K. Patra, Phys. Rev. D 100, 016009 (2019).
  * [11] A. Bandyopadhyay, B. Karmakar, N. Haque, M.G. Mustafa, Phys. Rev. D 100, 034031 (2019).
  * [12] S. Rath, B. K. Patra, JHEP 1712, 098 (2017).
  * [13] S. Rath, B. K. Patra, Eur. Phys. J. A 55, 220 (2019).
  * [14] B. Karmakar, R. Ghosh, A. Bandyopadhyay, N. Haque, M.G.Mustafa, Phys. Rev. D 99, 094002 (2019)
  * [15] Kenji Fukushima, Dmitri E. Kharzeev and Harmen J. Warringa, Phys. Rev. D 78, 074033 (2008).
  * [16] V. Braguta, M. N. Chernodub, V. A. Goy, K. Landsteiner, A. V. Molochkov and M. I. Polikarpov, Phys. Rev. D 89, 074510 (2014).
  * [17] Dmitri E. Kharzeev and Dam T. Son, Phys. Rev. Lett. 106, 062301 (2011).
  * [18] V. P. Gusynin, V. A. Miransky and I. A. Shovkovy, Phys. Rev. Lett. 73, 3499 (1994).
  * [19] P. Braun-Munzinger, J. Stachel, J. Phys. G 28, 1971 (2002)10.
  * [20] J. Cleymans, J. Phys. G 35, 044017 (2008)11.
  * [21] A. Andronic et al., Nucl. Phys. A 837, 65 (2010).
  * [22] K. Fukushima, Y. Hidaka, Phys. Rev. Lett. 117, 102301 (2016).
  * [23] M. Hasan, B. K. Patra, B. Chatterjee, P. Bagchi, Nucl. Phys. A 995, 121688 (2020).
  * [24] G.T. Bodwin, E. Braaten, G.P. Lepage, Phys. Rev. D 51, 1125 (1995).
  * [25] N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Nucl. Phys. B 566, 275 (2000).
  * [26] W. M. Alberico, A. Beraudo, A. De Pace and A. Molinari, Phys. Rev. D 77, 017502 (2008).
  * [27] M. Laine, O. Philipsen, and M. Tassler, P. Romatschke, J. High Energy Phys. 03, 054 (2007).
  * [28] T. Matsui, H. Satz, Phys. Lett. B 178, 416 (1986).
  * [29] A. Beraudo, J. P. Blaizot, C. Ratti, Nucl. Phys. A 806, 312 (2008).
  * [30] N. Brambilla, M. A. Escobedo, J. Ghiglieri, A. Vairo, JHEP 1305, 130 (2013).
  * [31] L. Thakur, U. Kakade, B. K. Patra, Phys. Rev. D 89, 094020 (2014).
  * [32] L. Thakur, N. Haque, U. Kakade, B. K. Patra, Phys. Rev. D 88, 054022 (2013).
  * [33] A. Dumitru, Y. Guo, M. Strickland, Phys. Rev. D 79, 114003 (2009).
  * [34] B. K. Patra, H. Khanchandani, L. Thakur, Phys. Rev. D 92, 085034 (2015).
  * [35] B. K. Patra, H. Khanchandani, Phys. Rev. D 91, 066008 (2015).
  * [36] A. Rothkopf, T. Hatsuda, S. Sasaki, Phys. Rev. Lett. 108, 162001 (2012).
  * [37] D. Lafferty and A. Rothkopf Phys. Rev. D 101, 056010 (2020).
  * [38] Y. Guo, L. Dong, J. Pan, and M. R. Moldes, Phys. Rev. D 100, 036011 (2019).
  * [39] J. Alford and M. Strickland, Phys. Rev. D 88, 105017 (2013).
  * [40] C. Bonati, M. D’Elia, and A. Rucci, Phys. Rev. D 92, 054014 (2015).
  * [41] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro, A. Rucci, and F. Sanfilippo, Phys. Rev. D 94, 094007 (2016).
  * [42] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro, A. Rucci, and F. Sanfilippo, Phys. Rev. D 95, 074515 (2017).
  * [43] M. Hasan, B. Chatterjee, B. K. Patra, Eur. Phys. J. C 77, 767 (2017).
  * [44] B. Singh, L. Thakur, H. Mishra, Phys. Rev. D 97, 096011 (2018).
  * [45] S. A. Khan, B. K. Patra and M. Hasan, arXiv:2004.08868[hep-ph].
  * [46] M. Hasan, B. K. Patra Phys. Rev. D 102, 036020 (2020).
  * [47] M. Hasan, B. K. Patra Int. J. Mod. Phys. A 36 (2021) 06, 2150045.
  * [48] B. Singh, S. Mazumder and H. Mishra, JHEP 2005, 068 (2020). .
  * [49] B. Singh, M. Kurian, S. Mazumder, H. Mishra, V. Chandra and S. K. Das, arXiv:2004.11092[hep-ph].
  * [50] B. Karmakar, A. Bandyopadhyay, N. Haque and M. G. Mustafa, Eur. Phys. J. C. 79, 658, (2019).
  * [51] T. Chyi et al., Phys. Rev. D 62, 105014 (2000).
  * [52] A. Ayala, C. A. Dominguez, S. Hernandez-Ortiz, L. A.Hernandez, M. Loewe, D. Manreza Paret, and R. Zamora, Phys. Rev. D 98, 031501 (2018).
  * [53] H. A. Weldon, Phys. Rev. D 26, 1394 (1982).
  * [54] R. D. Pisarski, Phys. Rev. Lett. 63, 1129 (1989).
  * [55] Uttam Kakade and Binoy Krishna Patra, Phys. Rev. C 92, 024901 (2015).
  * [56] Adrian Dumitru, Yun Guo, and Michael Strickland Phys. Rev. D 79, 114003 (2009).
  * [57] H. A. Weldon, Phys. Rev. D 42, 2384 (1990).
  * [58] A. Mocsy and P. Petreczky, Phys. Rev. Lett. 99, 211602 (2007).


# Pedaling, Fast and Slow
The Race Towards an Optimized Power Strategy

Steven DiSilvio, Anthony Ozerov, Leon Zhou

(February 21, 2022)

###### Abstract

With the advent of power-meters allowing cyclists to precisely track their
power outputs throughout the duration of a race, devising optimal power output
strategies for races has become increasingly important in competitive cycling.
To do so, the track, weather, and individual cyclist’s abilities must all be
considered. We propose differential equation models of fatigue and kinematics
to simulate the performance of such strategies, and an innovative optimization
algorithm to find the optimal strategy.

Our model for fatigue translates a cyclist’s power curve (obtained by fitting
the Omni-Power Duration Model to power curve data) into a differential
equation to capture which power output strategies are feasible. Our kinematics
model calculates the forces on the rider, and with power output models the
cyclist’s velocity and position via a system of differential equations. Using
fine-grained track data, including the slope of the track and velocity of the
wind, the model accurately computes race times given a power output strategy
on the exact track being raced.

To make power strategy optimization computationally tractable, we split the
track into segments based on changes in slope and discretize the power output
levels. As the space of possible strategies is large, we vectorize the
differential equation model for efficient numerical integration of many
simulations at once and develop a parallelized Tree Exploration with Monte-
Carlo Evaluation algorithm. The algorithm is efficient, running in
$O(ab\sqrt{n})$ time and $O(n)$ space where $n$ is the number of simulations
done for each choice, $a$ is the number of segments, and $b$ is the number of
discrete power output levels.

We present results of this optimization for several different tracks and
athletes. As an example, the model’s time for Filippo Ganna in Tokyo 2020
differs from his real time by just 18%, supporting our model’s efficacy.

###### Contents

  1. 1 Introduction
  2. 2 Problem Statement and Assumptions
  3. 3 Modeling Methodology
    1. 3.1 Track Data
    2. 3.2 Kinematic Differential Equation Model
    3. 3.3 Omni-Power Duration Model for Cyclists
    4. 3.4 Fatigue Differential Equation Model
    5. 3.5 Optimization: Tree Exploration with Monte-Carlo Evaluation
  4. 4 Model Results
    1. 4.1 Sample Solutions: Tokyo 2020 Time Trial
    2. 4.2 Sensitivity to Weather
    3. 4.3 Sensitivity to Strategy
    4. 4.4 Results for other courses
  5. 5 Model Evaluation
    1. 5.1 Overall Model
    2. 5.2 Differential Equation Model
    3. 5.3 Fatigue Differential Equation Model
    4. 5.4 Tree Exploration with Monte-Carlo Evaluation
  6. 6 Discussion of Multi-Rider Team Trials
  7. 7 Conclusion
  8. Letter to Directeur Sportif

##  1 Introduction

With the advent of power-meters allowing cyclists to precisely track their
power outputs throughout the duration of a race, devising optimal power output
distributions for races has become increasingly important in the world of
competitive cycling [19]. However, doing so involves considering a wide
variety of factors, ranging from the biological limitations of the particular
competitor to the topographical layout of the track the race will take place
on. In particular, as competitions have gotten more advanced riders have
become specialized for certain disciplines, such as sprinting and ultra-
marathons, leading to drastic differences in rider capabilities which all need
to be accounted for. Additionally, confounding factors such as inevitable
deviations in power output by cyclists from specific power routines and
weather variations on the day of the competition must also be considered when
designing any routine. An ideal power output distribution would minimize the
time necessary to complete a race subject to the constraint that the rider can
feasibly carry out such a routine to its completion.

##  2 Problem Statement and Assumptions

Given the power curve of a cyclist, which for any length of time gives the
maximum amount of constant power output they can sustain, as well as track-
specific data, we are tasked with devising a model to output an optimal power
distribution for completing a race in minimal time. In addition, we analyze
the effect of deviation from the optimal strategy on race time, and allow for
varying weather conditions to be accounted for in the model. In doing so, we
make the following assumptions:

  * •

A rider’s power curve is monotonically decreasing, and in particular is a
bijection from $\mathbb{R}$ to a closed subinterval of $\mathbb{R}^{+}$. This
is reasonable, as it would not be sensible if, for example, someone’s maximum
sustainable power output for 10 seconds is 1000W, and for 9 seconds is 800W.
This assumption is useful as it means the power curve is invertible, and a
unique duration can be computed given a power output level.

  * •

Cyclists are able to recover energy levels as they ride by cycling with less
power. In particular, each cyclist has a critical power output at which they
neither gain or lose energy, and which they could theoretically maintain
indefinitely. Cycling above this critical power will fatigue a cyclist over
time, and cycling below it will allow them to regain energy. This assumption
is supported by the literature, and has been used in several other studies
[4].

  * •

It is not feasible for a cyclist to have a power distribution plan in which
they must continuously vary their power, as doing so would be far too complex
for a rider to follow. Instead, the race must effectively be broken up into
discrete segments based on the distance a rider has covered, and during such a
segment a rider will be instructed to maintain a specific constant power
level.

  * •

A cyclist can effectively be viewed as traveling in one dimension, having a
velocity at any given instant which translates to a speed along the direction
of the path. Dimensional considerations such as varying elevations lead to the
exertion of forces on the rider, but these forces can be broken into
components parallel to the cyclist’s direction of travel (which will affect
their speed) and perpendicular to it (which will impose restrictions on the
velocity attainable to still remain along the course).

##  3 Modeling Methodology

We obtain track data in the form of GPX traces from La Flamme Rouge [1], a
European provider of cycling route data. This data contains latitude,
longitude, and elevation, and matches remarkably well with officially-
published race profile images.

###  3.1 Track Data

Figure 1: Tokyo 2020 track data visualized with local East (y), North (x), and
Up (Height) coordinates in the tangent plane at the track start point.
Elevation exaggerated.

#### Distances

Using this data, we calculate the cumulative distance $d_{i}$ (on the
horizontal plane) at each point $i$ along the track according to:

$d_{i}=\sum_{j=1}^{i-1}d_{j,j+1},$

where $d_{j,j+1}$ is the distance between points $j$ and $j+1$. We calculate
$d_{j,j+1}$ according to the haversine formula [13], which is used for
computing the lengths of great circles based on the latitude and longitude of
the two endpoints:

$d_{ij}=2r\arcsin\left(\sin^{2}\left(\frac{\lambda_{j}-\lambda_{i}}{2}\right)+\cos(\lambda_{i})\cos(\lambda_{j})\sin^{2}\left(\frac{\mu_{j}-\mu_{i}}{2}\right)\right),$

where $\lambda_{i},\lambda_{j}$ are the latitudes of points $i$ and $j$ and
$\mu_{i},\mu_{j}$ are the longitudes of points $i$ and $j$.

#### Slopes

If we let $a_{i}$ be the elevation at point $i$, also calculate the slope
angle $\theta_{i}$ at each point according to:

$\theta_{i}=\arctan\left(\frac{a_{i+k}-a_{i-k}}{d_{i+k}-d_{i-k}}\right)$

where $k$ is a smoothing factor, which we set to 5. This smoothing is
necessary as the elevation data is at one-meter resolution.

#### Curvature

One key track statistic is its curvature. As is explained in the section
below, the curvature metric required is the radius of the circle that the
curve is following. To determine this, we fit a circle to a small set of $k$
points around every point, as using three points to get the unique circle
would be less robust. We do this by solving the following least-squares
problem:

$\min_{x_{c},y_{c},r}\sum_{i=1}^{k}\left|r-\sqrt{(x_{i}-x_{c})^{2}+(y_{i}-y_{c})^{2}}\right|.\;\;\cite[cite]{[\@@bibref{}{circles}{}{}]}$

Effectively, given the set of points with coordinates
$(x_{1},y_{1}),\ldots,(x_{k},y_{k})$, we find the coordinates $(x_{c},y_{c})$
and radius $r$ that minimize the distance from the points to the circle. Note
that to get a radius in meters, the latitudes and longitudes of the points are
converted to coordinates in meters in the local tangent plane at the start
point of the track.

Figure 2: Radii of fitted circles on the Tokyo 2020 Men’s Road time trial
route.

#### Heading

To account for the effect of wind, we must know the velocity of the wind in
the direction of the cyclist. This depends solely on the windspeed and the
angle of the wind relative to their direction of motion. Thus, the direction
of motion at every point along the route must be calculated. With
$\lambda_{i},\lambda_{i+1}$ and $\mu_{i},\mu_{i+1}$ we may use the following
equation to calculate the heading $\gamma_{i}$ at every point, using its
coordinates and those of the next point:

$x=\sin(\mu_{i+1}-\mu_{i})\cos(\lambda_{i+1})$
$y=\cos(\lambda_{i})\sin(\lambda_{i+1})-\sin(\lambda_{i})\cos(\lambda_{i+1})\cos(\mu_{i+1}-\mu_{i})$
$\gamma_{i}=\text{atan2}(x,y).\;\;\cite[cite]{[\@@bibref{}{heading}{}{}]}$

#### Interpolation and gridding

The above computations give track data at every point in the GPX trace.
However, the models below require track data at arbitrary points. To account
for this, we pre-calculate interpolated track data on a 10m grid of points for
efficient access during the simulations. At each gridpoint $x$, a data point
$z(x)$ is calculated according to a linear interpolation between the two
nearest points $i$ and $i+1$:

$w(x)=\frac{d_{i+1}-x}{d_{i+1}-d_{i}}$ $z(x)=w(x)z_{i}+(1-w(x))z_{i+1}$

###  3.2 Kinematic Differential Equation Model

#### Motivation

Differential equation models lend themselves very naturally to the study of
kinematics. In particular, knowing a cyclist’s power output over the entire
duration of a race, as well as the topography of the track itself, it is
possible to set up differential equations modeling incremental changes in
their speed and position, allowing for the time they take to complete a race
to be numerically solved.

#### Preliminary assumptions

In order to simplify our model to be computationally feasible, the following
reasonable assumptions were made:

  * •

Riders can effectively be viewed as traveling in one dimension (i.e. along the
direction of the track), with all imposed forces able to be broken up into
components parallel and perpendicular to this direction of motion

  * •

Turning corners does not affect a rider’s speed along their direction of
motion, but instead imposes a limitation on the maximum speed attainable in
order to stay along the curve (this is in line with the fact that centripetal
force is always perpendicular to the velocity of a moving object, and hence
imposes no work)

  * •

Power outputs from power curves are the effective power that one is able to
produce, i.e. the power driving them forward, and hence internal power loss
(due to gear and chain friction) are not considered

Table 1: Differential equation model variables and constants. Variable | Description |  Units
---|---|---
$x$ | Distance covered along the track |  m
$x_{h}$ | Horizontal distance covered |  m
$v$ | Velocity |  m/s
$a$ | Acceleration |  m/s2
$\theta$ | Angle of incline |  radians
$P$ | Power output |  W
$\gamma$ | Cyclist heading |  radians
$v_{w}$ | Windspeed |  radians
$\gamma_{a}$ | Wind heading |  radians
$v_{a}$ | Effective airspeed against cyclist’s motion |  m/s
$F$ | Net force parallel to cyclist’s motion |  N
Constant |  Description | Value | Units |  Source/Rationale
$A$ |  Frontal area of cyclist | $.4$ | m2 |  [9]
$m$ |  Mass of cyclist | [various] | kg |  [2]
$g$ |  Gravitational acceleration | 9.81 | m/s2 |
$C_{d}$ |  Coefficient of drag | 0.6 | - |  [8]
$\mu_{s}$ |  Static friction coefficient | [various] | - |  [14]
$\mu_{r}$ |  Rolling resistance | .0025 | - |  [5]

#### Horizontal and Track Distance

At any given instance, where $x(t),x_{h}(t)$ are the distances traversed along
the track and horizontally respectively as a function of time $t$, we have
that:

$\frac{dx_{h}}{dt}=\frac{dx}{dt}\cos\big{(}\theta\left(x_{h}\right)\big{)}$

This makes intuitive sense, as assuming a constant angle $\theta$ of
inclination for some stretch of track, by elementary trigonometry clearly
$\Delta x_{h}=\Delta x\cos(\theta)$. Instantaneously at a given time $t$, the
angle $\theta(x_{h}(t))$ is computable via the method described in the
previous section.

#### Air speed

Given a cyclist’s velocity $v$ and heading $\gamma$, as well as the speed
$v_{w}$ and heading $\gamma_{w}$ of the wind, we calculate the effective speed
of the air moving against their motion, $v_{a}$ , as follows:

$v_{a}=v-v_{w}\cos(\gamma-\gamma_{w})$

In particular, under conditions with no wind (i.e. $v_{w}=0$) we have that
$v_{a}=v$, as the speed of the air moving towards the rider is simply the
speed with which they are traveling through the medium.

#### Force parallel to cyclist’s motion

Figure 3: Free body diagram of forces parallel to cyclist’s direction of
motion. Adapted from [12]

Forces acting parallel to the cyclist’s motion under consideration are those
due to gravity, rolling resistance, drag and their own power output (i.e.
forward propulsion) as seen in figure 3. Where $P$ is the power function
taking distance along the track as input which we wish to optimize for, the
forces are given as follows [4]:

$\displaystyle F_{a}$ $\displaystyle=\frac{1}{2}C_{d}Av_{a}^{2}$
$\displaystyle F_{g}$ $\displaystyle=mg\sin\big{(}\theta(x_{h})\big{)}$
$\displaystyle F_{r}$ $\displaystyle=\mu_{r}mg\cos\big{(}\theta(x_{h})\big{)}$
$\displaystyle F_{p}$ $\displaystyle=\frac{P(x)}{v}$

Thus, where each of $F_{p},F_{a},F_{g}$ and $F_{r}$ are positive magnitudes,
we have that the total force acting on the cyclist is given by:

$F=F_{p}-F_{a}-F_{g}-F_{r}$

where the positive direction is the direction of the cyclist’s travel.

#### Forward acceleration

Applying Newton’s Second Law of Motion and substituting the forces as above,
we immediately get that:

$\frac{dv}{dt}=\frac{F}{m}$

#### Maximum velocity around curves

While traveling around a banked curve, forces due to static friction and
gravity act perpendicular to the motion of the cyclist to change their
direction to keep them along the curve while not affecting the magnitude of
their velocity. However, given a portion of track with bank angle $\beta$ and
curvature radius $r$, where $r$ corresponds to the radius of the circle fitted
via the least squares method described above, the maximum attainable velocity
is given by:

$v_{max}=\sqrt{\frac{rg(\sin(\beta)+\mu_{s}\cos(\beta))}{\cos(\beta)-\mu_{s}\sin(\beta)}}\cite[cite]{[\@@bibref{}{bank}{}{}]}$

Traveling above $v_{max}$ would cause the rider to fly out of the curve, and
hence we cap their velocity at $v_{max}$ at any given point along the track as
we assume a cyclist would always break to assure they stay on the path.

###  3.3 Omni-Power Duration Model for Cyclists

#### Model

We calculate power curves of various athletes using the Omni-PD model [15] on
their public data. The model has been previously shown to be accurate in
modeling both high-intensity and endurance cycling. These athletes ranged from
amateurs to Olympians that have competed at the 2020 Tokyo time trials. Data
was pulled from Strava, and we took their maximal mean power levels at various
set intervals. This is the maximum energy they outputted within any interval
divided by the duration of that interval. Usually riders are able to generate
enormous power for seconds, but the power levels quickly taper out due to the
body building up too much lactic acid and will switch from anaerobic to
aerobic mechanisms. At longer intervals of exercise, cyclists enter into a
steady state, which is a relatively flat energy level.

All of this is accounted for by the model. Using non-linear least squares, we
can fit the curve to those points as seen in figure 4.

#### Deriving Critical Power and Work Capacity

By tuning these parameters to fit the power levels at certain intervals, we
can also deduce constants for each individual, such as their Critical Power
threshold $(P_{c})$, and their Anaerobic Work Capacity $(W^{\prime})$.
Critical power is the maximum sustainable energy over an extended period of
time. Usually for an athlete this will be around 30 minutes though it varies
by individuals [15]. Any power level above Critical power will drain from the
work capacity, and cause fatigue when the entire energy capacity is depleted
[6]. On the other hand, going below that threshold will allow riders to
recover their energy slightly.

Table 2: Power Curve variables and constants. Variable | Description |  Units
---|---|---
$P_{max}$ | Max Power |  W
$P_{C}$ | Critical Power |  W
$W^{\prime}$ | Work above $P_{C}$ (Anaerobic Work Capacity) |  W
$t$ | Time |  s
$T_{cpmax}$ | Time sustained at $P_{C}$ |  s
Constant |  Description |  |
$\beta$ |  Linear Constant | - |

$f(t)=\begin{cases}\frac{W^{\prime}}{t}*(1-e^{t*\frac{P_{max}-P_{C}}{W^{\prime}}})+P_{C}&t\leq
T_{cpmax}\\\ \\\
\frac{W^{\prime}}{t}*(1-e^{t*\frac{P_{max}-P_{c}}{W^{\prime}}})+P_{C}-\beta*ln(\frac{t}{T_{cpmax}})&t\geq
T_{cpmax}\end{cases}$ (1)

For each power level, we can take the inverse of the function to determine the
maximum time each athlete can sustain a certain power level for. These power
levels will become the choices that a rider can make at each segment of the
race.

#### Professional Data

All our simulations are based off the power curves we computed for three
professional cyclists shown in table 3. Their MMP were derived from power
recordings they achieved in either races or high-intensity practices.

Table 3: Rider Profile and Information [2] Name | Racer Type | Gender | Weight(kg)
---|---|---|---
Filippo Ganna | Time Trial Specialist | M | 82
Mathieu van der Poel |  Puncheur | M | 75
Chloe Dygert |  Time Trial Specialist | W | 67
Figure 4: Comparison of two different cyclist profiles, as fitted by the Omni-
PD model. The Puncheur, a cyclist who specializes in short steep hills, is
modelled off the data for Mathieu van der Poel. The time trial specialist is
modelled from Filippo Ganna’s rides. [17][18]

###  3.4 Fatigue Differential Equation Model

A cyclist cannot maintain a high level of power output indefinitely, and hence
developed a model for fatigue which tracks their energy level, constrained to
the interval $[0,1]$ and initialized to $1$, throughout the duration of the
race. As a cyclist’s power curve represents how long they can maintain a given
level of power output, we present the following simple model for fatigue,
where $f^{-1}$ is the inverse of a given cyclist’s power curve:

$\frac{dE}{dt}=\begin{cases}\frac{-1}{f^{-1}(P(x_{h}))}&\text{if
}P(x_{h})>P_{C}\\\ \frac{1}{7200P_{C}}(P_{C}-P(x_{h}))&\text{if }P(x_{h})\leq
P_{C}\end{cases}$

This model accounts for the fact that cyclist’s can theoretically sustain
their critical power indefinitely, while levels above this will diminish their
energy and levels below this will allow them to recover energy. In particular,
given a power level above $P_{C}$, the less time a cyclist can maintain it
based on their power curve the greater the rate of decrease in their energy
levels which outputting it would cause. The constant $\frac{1}{7200P_{C}}$ was
chosen such that after 2 hours (7,200 seconds) of not riding (i.e. maintaining
a power level of 0), a cyclist completely fatigued would have their full
energy restored.

In regards to the race simulation, a cyclist can never allow their energy
level to drop below 0, and hence this fatigue model provides a natural
constraint on power output. As a cyclist nears 0 energy, they are incentivized
to output lower power levels to either conserve or regain energy.
Additionally, a cyclist’s energy is capped at 1, and hence cycling below
$P_{C}$ can not indefinitely increase a cyclist’s energy.

###  3.5 Optimization: Tree Exploration with Monte-Carlo Evaluation

With models for fatigue and cycling kinematics, we may now optimize the power
strategy to minimize the time taken to complete the time trial.

#### Track segmentation

The slope of the track is a dominant factor in the cycling kinematics
equation. Indeed, studies of cycling power strategies often focus on the
course profile [19]. For this reason, we segment the track based on changes in
slope. This is done by taking local minima and maxima of elevation (which are
necessarily locations where the slope changes from negative to positive or
positive to negative, respectively) and local minima and maxima of slope—thus
we identify both changes in slope direction and slope extremes. Figure 5
displays the segmentation of the Tokyo route.

Figure 5: Segmentation of Tokyo 2020 Men’s time trial track, done using slope
and altitude data

#### Optimization problem

There are three factors which determine the most effective optimization method
for our formulation of the problem:

  1. 1.

Very large search space. With even just 20 segments and 5 power levels, there
are $5^{20}$ possible power strategies.

  2. 2.

The differential equation model can be vectorized to integrate an arbitrary
number $n$ of power strategies simultaneously in $O(\sqrt{n})$.

  3. 3.

Due to momentum and fatigue levels, the effect on finishing time of choosing a
power level for one segment depends on the choices for all of the other
segments.

#### Optimization method

Given these properties, we design a stochastic optimization method which is
linear in the number of segments and the number of power levels and converges
to the optimal solution with sufficiently large $n$. Starting from the first
segment, we repeatedly choose the power level for the next segment which
coincides with that in the optimal solution (i.e. that which minimizes time).
Since it is computationally intractable to test every possible power strategy,
we may instead choose $n$ random power strategies for the rest of the segments
and get their time results. We then use the minimum result of the random power
strategies to select the power level for the next segment. As $n\to\infty$,
the minimum result will be that of the optimal solution. Figure 6 illustrates
this process. In practice, rather than evaluating choices using the minimum
time, we use the $k$th lowest time. In the limit, this will still find the
optimal solution, but with smaller $n$ this makes the optimization less
sensitive to the stochasticity of the process.

400W⋮450W…500W400W450W……95.4…60.560.261.561.7…80.3 Figure 6: Diagram of the
optimization method. Having already chosen some sequence of power levels
($\ldots$, 400W, 500W) for a sequence of segments, we wish to decide between
400W, 450W, and some other power levels for the next segment. For each of
these choices, we randomly pick $n$ assignments of power levels for the
remaining segments, and calculate the time taken to cycle the route using each
of the random assignments. For 400W, the best time was 60.2 minutes, whereas
for 450W, it was 61.5 minutes. Thus we choose 400W as the power level for the
next segment. We repeat until we have chosen a power level for each segment.

##  4 Model Results

###  4.1 Sample Solutions: Tokyo 2020 Time Trial

We use an Euler approximation to numerically solve the Differential Equation
Model for the three cyclists studied on the Tokyo 2020 Time Trial route.
Figure 7 displays solutions for the Tokyo 2020 Time Trial for three different
athletes.

Figure 7: The optimal strategies for the three cyclists examined when
optimized for the Tokyo 2020 time trial route. Note that the women’s time
trial route is half the distance. We see that all of the strategies tend to
apply more power when going uphill and much less when going downhill. We also
see that the strategies of Filippo Ganna and Van Der Poel do differ slightly.
Ganna reserves bursts of power for major hills. Van Der Poel, a hill climber,
applies more power over less-steep hills as well, e.g. he applies a high power
level over the whole slope from 28,000m to 33,000m where Filippo does not.

###  4.2 Sensitivity to Weather

To determine the optimal strategy’s sensitivity to deviations in weather, we
may simulate the finishing times of the optimal strategy during different
weather conditions than those optimized for.

Figure 8: The left plot displays finishing times for a strategy for Filippo
Ganna in Tokyo 2020 optimized for a northward wind with speed 10m/s (indicated
by the red point) under a range of other windspeeds and directions. The
strategy still performs well under winds to the West and North. The strategy
performs significantly worse under winds blowing East. On the right, having
optimized instead for an eastward wind with speed 10m/s, we see that the
optimization does indeed take into account wind conditions, as much better
times are achieved with eastward winds than when optimizing for 10m/s North.

#### Wind

Ganna’s finishing time with the optimal strategy for a northward 10m/s wind is
4140 seconds. We run 10,000 simulations of this strategy with wind direction
chosen uniformly randomly and windspeed chosen uniformly randomly in the
interval $[0,10]$. Figure 8 displays the finishing times of these simulations.
We find that differences in windspeed and direction do affect the finishing
times, but strategies (whose optimization is not based solely on windspeed),
can still perform as well or better under a range of weather conditions.

Table 4: Filippo Ganna’s time results (in seconds) under the optimal strategy in the Tokyo 2020 time trial under different weather conditions. Note that, when the weather conditions match those optimized for, the times are roughly the same. When the weather conditions do not match, times worsen significantly. Optimized weather | No rain in simulation | Rain in simulation
---|---|---
No rain | 4105 | 4315
Rain | 4416 | 4099

#### Rain

Previous results have shown that the rolling resistance increases with the
thickness of water on the road. This increase is approximately $30\%$ with
light rain conditions of $0.3$mm of water on the road surface[11]. Table 4
displays the time results of optimizations for different rain conditions, and
the results when the actual rain conditions do not match those optimized for.
Rain, by increasing rolling resistance (making it difficult to go fast on
straightaways) and decreasing static friction (forcing cyclists to slow down
more around corners) has a universally negative effect under our model. Yet we
see that, by optimizing with consideration of the weather conditions, its
effects can be mitigated to the point that the stochastic optimization
produces a slightly better finishing time for rain than for no rain.

Figure 9: Filippo Ganna’s split times for the split from to 25,753m to 32,360m
in 1000 simulations with deviations from the optimal power strategy. In every
segment, the deviation from the optimal power level $P$ is calculated under a
Normal distribution with a mean of zero and a standard deviation of $P/50$.
Under the optimal strategy, the split time is 858.5s (indicated with the red
dashed line). Note that most split times fall around the optimal time.
However, some are significantly slower. These are likely choices for power
levels which make Ganna exhausted by the time this split is reached.

###  4.3 Sensitivity to Strategy

Given a power strategy, the expected split times and deviations from them may
be computed. One split of interest in the Tokyo 2020 men’s time trial race is
the time it takes to climb the slope that runs from about 26,000m to 32,000m.
Figure 9 displays Filippo Ganna’s split time for this split when deviations
from the optimal strategy are made. We see that the model is indeed sensitive
to how closely the cyclist follows the strategy, as the optimization naturally
creates strategies which leave energy at 0 by the end of the time trial. If a
cyclist exceeds the target power too much, they will be fatigued and must bike
at their baseline speed.

Figure 10: Filipo Ganna’s power curve in Flanders, Belgium. This track has him
mostly riding at the same power, and this is due to the track being mostly
flat and a straight course without many curves. This lends itself to an
optimal strategy where constant power (though with a burst of power at the
beginning) leads to a better time.

Figure 11: Filipo Ganna’s power curve in a custom designed map of NYC Central
Park. This map is designed to go around the bike path across all of the park,
and will have sharp curves as the reach the edges. Along the path there are
various short hills where he is seen changing his power output to match the
slope.

###  4.4 Results for other courses

Since these strategies are tailored to the rider and we are able to test out
multiple riders on the same Tokyo track, the next benchmark is being able to
place the same rider on different tracks and have them optimize their power
for each track respectively. The results below in Figure 10 show Filippo
Ganna’s optimal racing plan for the UCI 2021 World Championships track in
Flanders. The drastic difference in racing plan is proof that the rider is
also optimizing their plan to tackle the unique aspects of each individual
race, but also that this model is general enough to accommodate different
riders and different courses.

##  5 Model Evaluation

###  5.1 Overall Model

#### Strengths

  * •

The model gives results that roughly coincide with real-world results. For
instance, Ganna finished in 56 minutes in the Tokyo 2020 time trial [3], which
is reasonably close to his finishing times when weather matches that expected
in Table 4 (around 68 minutes). This goes to show that our models for power
curves, fatigue, and kinematics are roughly correct.

  * •

Using real-world data for routes and power curves makes the model directly
applicable to specific races and athletes.

#### Limitations

  * •

For weather and rain, it assumes that these effects are constant throughout
the race. In reality, wind can change direction and rain will either increase
or decrease over time. Since our model returns a power based on the position,
it would be hard to constantly update weather especially when the rider might
deviate from split times.

  * •

The MMP at certain intervals are from race and practice data, which might not
accurately reflect a cyclist’s power curve as well as laboratory tests would
show.

###  5.2 Differential Equation Model

#### Strengths

  * •

The kinematics model is very well-behaved. For instance, in Figure 8, the
finishing times change smoothly as windspeed and heading change.

  * •

The model runs in $O(\sqrt{n})$, where $n$ is the number of simulations done.

#### Limitations

  * •

As the forwards force from the cyclist’s power is determined by $P(x)/v$, the
model is numerically unstable when $v$ is small. This is solved by integrating
with sufficiently small timesteps.

###  5.3 Fatigue Differential Equation Model

#### Strengths

  * •

The model is simple and intuitively matches what a power curve represents.

  * •

The rate at which energy decreases is directly determined by the power curve,
and not through parameters that roughly summarize a power curve.

#### Limitations

  * •

We assume a linear rate of recovery in relation to power for power levels
below critical power, which is most likely not entirely accurate

###  5.4 Tree Exploration with Monte-Carlo Evaluation

#### Strengths

  * •

The optimization algorithm is able to search a very large search space
efficiently. Notably, it is linear in the number of segments the track is
split into, as well as the number of power levels evaluated for each segment.
It thus perfectly matches the optimization problem at hand.

  * •

The algorithm takes advantage of the differential equation model’s ability,
when vectorized, to efficiently solve for the finishing times of many power
strategies at once.

  * •

The algorithm quickly converges to an optimal solution. For example, changing
$n$ from 100 to 1,000 only improves the solution by 23 seconds.

#### Limitations

  * •

As optimization is only done over discrete power levels, the model may miss
better power strategies which take on intermediate power levels.

  * •

Similarly, as optimization is only done over discrete segments of the track,
the model may miss better power strategies which continuously vary power
output over the whole track.

  * •

As noted above, the optimal strategy runs cyclists’ energy down to zero by the
end of a race, which is obvious as, if any energy is remaining by the end, the
cyclist could have outputted more power and gone faster earlier. Any
deviations from the plan could cause the cyclist to be completely fatigued by
the end of the race. As the optimization algorithm takes the $k$th best time
for a choice as the time to evaluate a node with, it cannot naturally take
into account randomness in execution like traditional Monte Carlo Tree Search.
In practice, during a race, a cyclist could try to err on the side of caution,
and leave some energy to the end, and indeed we often see cyclists using up
extra energy in a sprint at the end of a time trial.

##  6 Discussion of Multi-Rider Team Trials

Another component of cycling competitions are the team time trials, where 6
riders compete in a team and are judged by the finishing positions of the
first 4 in the group. This leads to alternative strategies where riders
optimize for the team rather than individual results, and can create room for
specialized cyclists such as the Domestique [16], who leads the team at the
beginning and allows them to conserve energy by taking on all of the wind
resistance.

Optimal results can be simulated by including the power curves of each rider
in the team, and optimizing for the finish speed of any 4 riders. With the
assumption that the riders will be in a linked formation, only the front rider
will experience the effects of wind resistance, and the rest can cycle in the
slipstream generated. One example would be allowing the Domestiques to exceed
their power at the start and run out of energy, if it allowed the other four
riders to break out and sprint towards the finish line.

Different strategies that can be simulated include rotating the lead position
between the riders, which shares the burden of wind resistance equally across
all riders. Since the model optimizes for minimizing time, by including a
choice to rotate positions at each segment, the model will also optimize the
rotation between riders.

By running the optimization over a large set of riders, it could also find the
optimal 6 riders to fit within a team, based on how their combination
optimizes the minimum time.

##  7 Conclusion

Using the innovative Tree Exploration with Monte-Carlo Evaluation algorithm,
we were able to provide a race plan personalized for the rider and the
specific track they are competing on. At any point on the course, the rider
will know how much power to output, as well as their expected times for
reaching different checkpoints and finishing the race overall.

This model incorporates information about the rider’s power levels, weight,
and individual stamina calculated from the Omni-PD model. Data for the tracks
is precise, and allows the model to consider terrain conditions at each
specific point in the race, rather than merely generalizing from similar
tracks. Using the differential equation model, these factors at instantaneous
positions along the track are able to be translated into instantaneous rates
of change in state conditions, namely cyclist energy levels, distance covered
along the track, and velocity. By modeling these incremental changes over
time, we are able to get a very accurate idea of the performance of individual
power strategies, providing a framework for the optimization process.

Also included as inputs are race-day wind direction and weather conditions,
allowing the Directeur Sportif to prepare for all possibilities during the
competition. The plan comes with expected split times at each segment, and a
range of deviations depending on how closely the plan was followed. This is
useful, as riders realistically will differ slightly from the optimal plan,
but will still be able to get an estimation of their lap time. The resulting
race plan is promising, due to its incorporation of real terrain and
physiological factors, robustness to deviations, and ease of access to
cyclists.

## References

  * [1] _La flamme rouge_ , https://www.la-flamme-rouge.eu.
  * [2] _Pro cycling stats_ , https://www.procyclingstats.com.
  * [3] BBC Sport (2021).
  * [4] Faraz Ashtiani, Vijay Sarthy, M Sreedhara, Ardalan Vahidi, Randolph Hutchison, and Gregory Mocko, _Experimental modeling of cyclists fatigue and recovery dynamics enabling optimal pacing in a time trial_ , Jul 2020.
  * [5] Jarno Bierman, _Continental grand prix 5000 s tr rolling resistance review_ , Oct 2021.
  * [6] Matthew I. Black, Andrew M. Jones, Jamie R. Blackwell, Stephen J. Bailey, Lee J. Wylie, Sinead T. J. McDonagh, Christopher Thompson, James Kelly, Paul Sumners, Katya N. Mileva, and et al., _Muscle metabolic and neuromuscular determinants of fatigue during cycling in different exercise intensity domains_ , Journal of Applied Physiology 122 (2017), no. 3, 446–459.
  * [7] I. D. Coope, _Circle fitting by linear and nonlinear least squares_ , Journal of Optimization Theory and Applications 76 (1993), no. 2, 381–388.
  * [8] Timothy N. Crouch, David Burton, Zach A. LaBry, and Kim B. Blair, _Riding against the wind: a review of competition cycling aerodynamics_ , Sports Engineering 20 (2017), no. 2, 81–110.
  * [9] P. Debraux, W. Bertucci, A. Manolova, S. Rogier, and A. Lodini, _New method to estimate the cycling frontal area_ , International Journal of Sports Medicine 30 (2009), no. 04, 266–272.
  * [10] Andrew Duffy, _More circular motion_ , http://physics.bu.edu/~duffy/py105/MoreCircular.html, 2006.
  * [11] Jerzy Ejsmont, Leif Sjögren, Beata Świeczko Żurek, and Grzegorz Ronowski, _Influence of road wetness on tire-pavement rolling resistance_ , Journal of Civil Engineering and Architecture 9 (2015), no. 11.
  * [12] Mathis Fenre and Alex Klein-Paste, _Rolling resistance measurements on cycleways using an instrumented bicycle_ , Journal of Cold Regions Engineering 35 (2021).
  * [13] Kenneth Gade, _A non-singular horizontal position representation_ , Journal of Navigation 63 (2010), no. 3, 395–417.
  * [14] Matthew Michaels, _Coefficients of friction for concrete - the physics factbook_ , 2006.
  * [15] Michael J. Puchowicz, Jonathan Baker, and David C. Clarke, _Development and field validation of an omni-domain power-duration model_ , Journal of Sports Sciences 38 (2020), no. 7, 801–813.
  * [16] Joe Robinson, _The invisible man: The life of a domestique_ , Mar 2019.
  * [17] Strava, _Filippo power curve_ , https://www.strava.com/activities/5324258172/power-curve, May 2021.
  * [18] , _Ronde van vlaanderen - mathieu van der poel’s bike ride_ , https://www.strava.com/activities/5071352854/power-curve, Apr 2021.
  * [19] David Sundström, Peter Carlsson, and Mats Tinnsten, _On optimization of pacing strategy in road cycling_ , Procedia Engineering 60 (2013), 118–123.
  * [20] Elliott White and Jonathan Eisenhamer, _Php 5 in practice_ , Sams, 2007.

## Letter to Directeur Sportif

Dear Directeur Sportif,

As you may very well know, on race day mere seconds differentiate the cyclist
who comes out on top from those who don’t even place. As athletes near the
limits of human potential, utilizing data regarding individual cyclist
capabilities as well as the track on which they’ll compete will become
increasingly necessary in order to devise personalized strategies to gain a
competitive edge. Metrics for cyclist performance are getting evermore
precise, and hence it’s only natural that the way it is utilized to optimize
race strategy keeps up.

The model takes in the specific characteristics of a cyclist and the track to
create a racing plan tailored specifically for them. This way, their
individual strengths are utilized to the utmost, and the rider can be more
confident in the plan.

This is done by creating a representation of the cyclist, and repeatedly
simulating them racing the course thousands of times. At each key part of the
race, they are able to decide how much power they will cycle with for the next
segment. Simulations of the possible finishing times for each choice are done,
and the power level that leads to the best time is selected. The fastest time
for the course is then saved, along with the power the cyclist needs to expend
at each part of the race to achieve that time. The result is a racing plan
that is effective and simple to follow.

Attached to this letter is the map of next year’s 2022 UCI World Championship
Men’s Individual Time Trial event at Wollongong, Australia. Next to it is a
race plan we have designed specifically for Filippo Ganna, the time trial
specialist, on this race track.

Figure 12: UCI 2022 Cycling track and Filippo Ganna’s optimal strategy. Under
this strategy, he can finish the race in 51 minutes

You can see that in the middle chart that Filippo’s power output increases for
climbs up hills and decreases on descents, which matches common intuition. The
strength of this model lies in the precision of how much power should be used
at each segment. On the top is also a chart of his energy over time, which
decreases accordingly to how hard he cycles in each part. The power output is
segmented by key parts of the race, so that it is easier for the rider to
follow the plan, and at any point he can quickly glance at exactly how much
power he needs.

We also account for the maximum speed he can take at turns, and the
acceleration he gets going downhill from gravity, to ensure his safety.

Figure 13: Distance of the track plotted against the time he will take to
reach each split.

As part of the model, we have split times for each of the key segments. The
total projected time he will take to finish the course is 51 minutes. These
numbers can be used to track his progress along the race and to see if he is
above or below target. One example is knowing that he will reach the peak of
the first large hill 9:45 mins into the race, and that at that point he should
be putting in more than 500 watts of power. The model is able to compute
deviations from this, in case he misses the target, and can give an expected
range of when he will reach each segment.

As the event gets closer, you get a better sense of the weather conditions on
race-day. These can be factored into the calculations by giving the model
expected wind speed, wind direction, and rain. All of this will ensure that
you will get the most realistic model possible and the best strategy for
Filippo.

Best Regards,
The Bicycling Optimization Interest Society